# Parameterization and optimizability of pulse-level VQEs

Kyle M. Sherbert  Department of Chemistry, Virginia Tech, Blacksburg, VA 24061
Department of Physics, Virginia Tech, Blacksburg, VA 24061  Virginia Tech
Center for Quantum Information Science and Engineering, Blacksburg, VA 24061,
USA  Hisham Amer  Sophia E. Economou  Edwin Barnes  Department of Physics,
Virginia Tech, Blacksburg, VA 24061  Virginia Tech Center for Quantum
Information Science and Engineering, Blacksburg, VA 24061, USA  Nicholas J.
Mayhall<EMAIL_ADDRESS>Department of Chemistry, Virginia Tech, Blacksburg,
VA 24061  Virginia Tech Center for Quantum Information Science and
Engineering, Blacksburg, VA 24061, USA

###### Abstract

In conventional variational quantum eigensolvers (VQEs), trial states are
prepared by applying series of parameterized gates to a reference state, with
the gate parameters being varied to minimize the energy of the target system.
Recognizing that the gates are intermediates which are ultimately compiled
into a set of control pulses to be applied to each qubit in the lab, the
recently proposed ctrl-VQE algorithm takes the amplitudes, frequencies, and
phases of the pulse as the variational parameters used to minimize the
molecular energy. In this work, we explore how all three degrees of freedom
interrelate with one another. To this end, we consider several distinct
strategies to parameterize the control pulses, assessing each one through
numerical simulations of a transmon-like device. For each parameterization, we
contrast the pulse duration required to prepare a good ansatz, and the
difficulty to optimize that ansatz from a well-defined initial state. We
deduce several guiding heuristics to implement practical ctrl-VQE in hardware,
which we anticipate will generalize for generic device architectures.

††preprint: APS/123-QED

##  I Introduction

Variational quantum eigensolvers (VQEs) are among the most promising
candidates for achieving useful computations in chemistry on near-term quantum
computers [1, 2, 3, 4, 5, 6]. At their core, they are predict-and-test
methods, where a quantum state, determined by a set of classical parameters as
specified by an ansatz, is prepared on the quantum computer, and its energy
measured. Then a new quantum state is prepared with a new set of parameters,
selected by any classical optimization algorithm, and the process is repeated
until the energy is minimized. This minimal energy is an upper-bound to the
ground-state energy of the system, and, if the optimization is carried out
perfectly, the final state is the best approximation to the ground state
attainable with the chosen ansatz.

Clearly, a key decision in any VQE is the choice of ansatz, typically
expressed as a series of parameterized unitary gates in a quantum circuit.
While the original formulation of VQE utilized an ansatz inspired by the
underlying chemistry [7], many novel ansätze have been proposed which balance
chemical intuition with information about which parametric and entangling
operations are most easily implemented on the device [8, 9, 10, 11, 12, 13,
14]. In either case, for experimental realization, the quantum circuit must be
compiled into a sequence of electromagnetic pulses that are applied to the
device. For transmon-type devices (those considered in this work), the control
pulses are microwave signals which drive the states of each physical qubit as
prescribed by each unitary gate in the circuit. Unfortunately, if each gate is
to be implemented with high fidelity in modern hardware, even the most compact
ansätze for moderately-sized systems typically compile into control pulses
with a duration far exceeding the coherence time of the device [15].

Three of the authors previously proposed the algorithm ctrl-VQE [15, 16],
which takes the idea of a hardware-efficient ansatz to the extreme by
parameterizing the actual physical control pulses used in the lab, bypassing
the use of gates entirely. Designing the ansatz at the pulse level allows
drastically shorter evolution times, even approaching the quantum speed limits
imposed by the hardware [17], and hypothetically enabling the VQE to study
much larger or more complex systems. The methods and mathematical theory
behind ctrl-VQE are built upon the closely related field of quantum control,
and we refer the reader to Refs. [18, 19] for introductions to this topic.
Recently, several works have built on the original ctrl-VQE [20, 21, 22, 23,
24, 25, 26, 27], including a handful with proof-of-concept hardware
demonstrations on simple systems.

Most of these efforts have focused on parametric modifications of established
one- and two-qubit gates, while our original work has emphasized de novo
pulses which act more like a multi-qubit gate. In fact, ctrl-VQE represents
the most flexible variational ansatz possible for any given hardware. While
this approach is very powerful, it presents the challenge that there is an
overwhelming number of accessible degrees of freedom. There is typically one
or more control pulses simultaneously applied to each qubit, and each control
pulse is parameterized by an amplitude ($A$), a phase ($\phi$), and an
(angular) frequency ($\nu$), each of which may be varied in time, resulting in
an ostensibly infinite number of parameters (before hardware bandwidth limits
are considered). An important open question is which parameters are most
important for quickly preparing a target state while avoiding optimization
issues.

In this paper, we investigate the impact that each of the pulse parameters
tends to have and establish useful heuristics for practical pulse
parameterizations. For example, we present empirical evidence that
representing the amplitude $A$ and phase $\phi$ together as a complex
amplitude with real and imaginary components enables notably more efficient
optimizations than when varying each pulse parameter independently.
Furthermore, we find that varying the phase $\phi$ is likely to be far more
effective than varying drive frequency $\nu$ in a practical ctrl-VQE
experiment, since varying the latter tends to significantly impede the
efficiency of optimization, without providing any commensurate improvement in
controllability or minimal evolution time. We mainly focus on transmon-type
architectures in this work, but we anticipate our conclusions can guide pulse-
level VQE experiments on a wide array of devices.

The paper is organized as follows: In Section II, we review the mathematical
basis of the ctrl-VQE algorithm, and we introduce the various choices of
parameterization for our ansatz. In Section III, we briefly describe the
methods and parameters used to generate the data presented in this paper. In
Section IV, we present a series of experiments contrasting the optimizability
of each parameterization. Section V summarizes our results and condenses our
findings into a set of guidelines for experimental realizations of ctrl-VQE.

##  II Theory

We take as our starting point a quantum observable $\hat{O}$ (e.g., a
molecular Hamiltonian), and a reference state $\ket{\psi_{\mathrm{0}}}$ (e.g.,
the Hartree-Fock state), both mapped onto qubits via a suitable transformation
(e.g., Jordan-Wigner). Our variational ansatz $\ket{\psi({\bf\it\theta})}$
will take the form:

$\ket{\psi({\bf\it\theta})}=U({\bf\it\theta})\ket{\psi_{\mathrm{0}}},$ (1)

where $U({\bf\it\theta})$ is the unitary state-preparation operator. The
expectation value of $\hat{O}$ is the energy $E({\bf\it\theta})$:

$E({\bf\it\theta})\equiv\expectationvalue{\hat{O}}{\psi({\bf\it\theta})},$ (2)

which is the function that the optimization routine in ctrl-VQE will attempt
to minimize.

In ctrl-VQE, the unitary $U({\bf\it\theta})$ is the physical evolution under a
device Hamiltonian (as distinct from traditional VQEs, where it is typically a
sequence of logical operations on the qubit space). The device Hamiltonian may
be separated into a static component $\hat{H}_{0}$ and a drive component
$\hat{V}(t)$. Variational parameters are typically associated with the time-
dependent drive $\hat{V}(t)$. Frequently, $\hat{V}(t)$ is comprised of several
independent drives; in this work, we consider an independent drive
$\hat{V}_{q}(t)$ applied to each qubit $q$. As such, our choice of ansatz in
ctrl-VQE is essentially the unitary time evolution operator:

$U({\bf\it\theta})=\mathcal{\hat{T}}\exp(-i\int_{0}^{T}\mathop{}\\!\mathrm{d}t~{}\quantity[\hat{H}_{0}+\sum_{q}\hat{V}_{q}({\bf\it\theta},t)]),$
(3)

where $\mathcal{T}$ is the time-ordering operator and $T$ is the total pulse
duration. While the underlying physical mechanism varies for different quantum
computing architectures, the drive term in a qubit Hamiltonian can typically
be modeled as an interaction between a dipole moment (the qubit) and a
classical electric field (the control field):

$\hat{V}_{q}(t)\propto D(t)(\hat{a}_{q}+\hat{a}^{\dagger}_{q}),$ (4)

where $\hat{a}_{q}$ is the bosonic lowering operator acting on the transmon
qubit state. The electric field arises from an applied voltage that oscillates
sinusoidally at microwave frequencies:

$D(t)=A\cos(\nu t+\phi).$ (5)

The parameters we have control over are therefore the amplitude $A$, the
frequency $\nu$, and the phase $\phi$ of the field, each of which may
themselves be functions of time.

When the amplitude is relatively weak with respect to the drive frequency, and
when the difference between the drive frequency and the qubit frequency is
negligible compared to either of these quantities, the rotating wave
approximation can be used to convert the sinusoidal form of Eq. (4) into a
complex phase:

$\hat{V}_{q}(t)\xrightarrow{\mathrm{RWA}}\Omega e^{i\nu
t}\hat{a}_{q}+\mathrm{h.c.},$ (6)

where the complex amplitude $\Omega$ is defined in polar notation by the
amplitude (i.e., modulus) and phase:

$\Omega\equiv Ae^{i\phi}/2.$ (7)

Alternatively, we could parameterize our control field equivalently in
“Cartesian” notation, defining the real-valued amplitudes $\alpha$ and $\beta$
as the real and imaginary parts of $\Omega$:

$\Omega\equiv\alpha+i\beta.$ (8)

This Cartesian parameterization is often more convenient in software
implementations.

In the absence of any control field, a single qubit will evolve under its
static Hamiltonian $H_{0}$, which in transmons can be modeled as an anharmonic
oscillator with a characteristic frequency $\omega$ describing the energy gap
between the lowest two states $\ket{0}$ and $\ket{1}$. In such cases, it is
common to make a coordinate transformation from the lab frame into the
rotating frame of the device.

$U^{(t)}_{\rm RF\leftarrow LAB}=e^{iH_{0}t}$ (9)

The Hamiltonian in this frame is written as:

$V_{R}(t)=\Omega e^{i\Delta t}\hat{a}+\mathrm{h.c.},$ (10)

where the detuning $\Delta$ is simply the difference between the drive
frequency $\nu$ and the qubit frequency $\omega$:

$\Delta\equiv\nu-\omega,$ (11)

so that $\Delta=0$ corresponds to driving on resonance. Because the drive
frequency $\nu$ and the detuning $\Delta$ differ by a constant, they will
behave identically under a gradient-based optimization. In what follows, we
choose to use $\Delta$ as our variational parameter. In multi-qubit systems,
static coupling between qubits leads to a dressing of the bare qubit
parameters, but the detuning $\Delta$ of any given pulse from the frequency
$\omega$ of the qubit it targets remains a useful parameter.

In our previous work [15, 16], we took the amplitude $\Omega$ to be real,
essentially restricting the phase $\phi$ to either 0 (when $\Omega\geq 0$) or
$\pi$ (when $\Omega<0$). In the present work, we instead consider various
parameter sets, including complex amplitudes in some cases; all parameter sets
are listed in Table 1. Each parameter can be varied in an optimization loop,
between control pulses. Furthermore, each parameter could in principle be
varied as an arbitrary function in time, throughout the pulse duration. As a
practical matter, drive frequencies are more difficult to experimentally vary
smoothly, so in this paper as in our previous papers we will take $\Delta$ as
time-independent for any given control pulse. For simplicity, we will limit
ourselves to “windowed” pulses where each window has a constant amplitude and
phase; arbitrary pulse shapes can be approximated by taking windows of
arbitrarily short duration.

Label | $\Omega$ | $\nu$ | Figs
---|---|---|---
$\quantity{\alpha\beta}$ | 2 param. per window, Cartesian | Resonant. |  1, 2, 4, 5
$\quantity{\alpha\beta}_{2}$ | Same as $\quantity{\alpha\beta}$, twice as many windows. | 2
$\quantity{\alpha\beta}_{\infty}$ | Same as $\quantity{\alpha\beta}$, continuous windows. | 2
$\quantity{A\phi}$ | 2/window, polar | Resonant. | 4
$\quantity{\alpha}$ | 1/window, real only | Resonant. | 5
$\quantity{\alpha\Delta}$ | 1/window, real only | 1/pulse | 5
$\quantity{\alpha\beta\Delta}$ | 2/window, Cartesian | 1/pulse | 5
Table 1: Summary of pulse parameterizations considered in this paper. Except
when otherwise noted, optimizations begin with all complex amplitudes $\Omega$
initialized to zero, and drive frequencies $\nu$ initialized on resonance
(i.e., $\Delta=0)$.

##  III Methods

The results presented in this paper focus on the lithium hydride molecule
(LiH) mapped onto four qubits as a case study, usually with a nuclear
separation of 3.0Å. We use the pyscf package [28] to calculate electronic
integrals for each molecular geometry, and we use the qiskit package [29] to
construct the second-quantized molecular Hamiltonian, map it onto qubits with
the parity mapping, and apply two-qubit reduction. [30] We use the Hartree-
Fock singlet state as our reference, and we represent the Hamiltonian as a
matrix in the Hartree-Fock basis so that the reference state is a
computational basis state.

We use our own Julia code [31, 32] to simulate unitary time evolution under
the rotating wave approximation with arbitrary control pulses, using
Trotterized time evolution with twenty time steps per nanosecond. Before the
pulse, the qubits are prepared into the Hartree-Fock reference state, and
after the pulse, the resulting statevector (in the rotating frame) is used to
measure the energy, i.e., the expectation value of our molecular Hamiltonian.
We measure in the rotating frame so that when all pulse amplitudes are set to
zero, the cost function (molecular energy) does not change in time.

The static Hamiltonian $H_{0}$ is modeled as a set of coupled harmonic
oscillators, with each qubit truncated to the first two energy levels:

$\hat{H}_{0}=\sum_{q}\omega_{q}\hat{a}^{\dagger}_{q}\hat{a}_{q}+\sum_{\expectationvalue{p,q}}g_{pq}(\hat{a}^{\dagger}_{p}\hat{a}_{q}+\hat{a}^{\dagger}_{q}\hat{a}_{p}).$
(12)

Note that this is a common model for transmon devices, [33, 34] excepting an
additional anharmonicity term
$\hat{a}^{\dagger}_{q}\hat{a}^{\dagger}_{q}\hat{a}_{q}\hat{a}_{q}$ which is
dropped in the two-level approximation. In this paper, coupling strengths are
fixed to $g_{pq}=0.02\cdot 2\pi\mathrm{~{}GHz}$ if $q=p+1$ (i.e., linear
nearest-neighbor coupling), and qubit frequencies are equally spaced such that
$\omega_{q}=(4.80+0.02q)\cdot 2\pi\mathrm{~{}GHz}$. These parameters roughly
match those found in real IBMQ devices, while being systematically scalable to
larger systems.

We use the Optim.jl implementation of BFGS (a quasi-Newton optimization
algorithm) [35] with analytical gradients (see Appendix A and [36]) to
iteratively minimize the molecular energy with respect to a specified drive
parameterization. Hardware bounds on the maximum amplitude are enforced by
including a smooth penalty term in the cost-function on any time step for
which the pulse amplitude exceeds $\Omega_{\rm max}=0.02\cdot
2\pi\mathrm{~{}GHz}$. Optimization terminates when the maximum of the
analytical gradient has converged to within 1e-6.

##  IV Results

Previous works on ctrl-VQE have emphasized the role of the amplitude $A$ as a
function of time. In this work, we will assess the relative significance of
the phase $\phi$ and the detuning $\Delta$, as measured by the pulse duration
necessary to achieve high-accuracy estimates of a ground state energy. We will
make extensive use of “error-vs-time” plots and “iterations-vs-time” plots,
which show how the accuracy and iteration count of a ctrl-VQE optimization
change as a function of pulse duration.

(a)  (b)  (c)  (d)

Figure 1:  (1a) Ground state energy error vs. pulse duration for a ctrl-VQE
simulation of LiH with nuclear separation 3.0Å, using the
$\quantity{\alpha\beta}$ parameterization. Each pulse is divided into as many
uniformly spaced windows as possible such that window length satisfies $s\geq
3.0\mathrm{~{}ns}$. For the black curve, optimization is carried out with all
parameters starting at zero. Gray dashes mark optimized runs when parameters
are initialized with random values. (1b) The number of BFGS iterations needed
to obtain the energies in (1a). (1c) Optimized pulse amplitudes $A(t)$ for the
ctrl-VQE run at $T=21\mathrm{~{}ns}$. Each color corresponds to the drive on a
different qubit. The dashed line marks the bound above which penalties are
imposed, ostensibly the maximum drive allowed by the device. (1d) Optimized
pulse phases $\phi(t)$ for the ctrl-VQE run at $T=21\mathrm{~{}ns}$. The
dashed lines occur at intervals of $\pi/2$, separating quadrants of the unit
circle.

Each data point in these plots is obtained by simulating ctrl-VQE to optimize
parameters in a pulse with fixed duration $T$. At the start of optimization,
pulse parameters are initialized to zero, meaning the pulse shape starts flat
and the drive frequencies start on resonance. The pulse is divided into
windows with uniform spacing $s$, with longer pulses divided into more windows
(details given in Section IV.1). In the “error-vs-time” plots, we plot the
difference between the final optimized energy and the value obtained from
exact diagonalization of the Hamiltonian matrix. In the “iterations-vs-time”
plots, we plot the number of iterations the BFGS optimizer requires to
converge. This number is unique to our choice of optimizer, but it allows a
heuristic comparison of the expected number of circuit evaluations required to
obtain the accuracies presented in the corresponding “error-vs-time” plot.

Typical behavior in an “error-vs-time” plot is for the error to start off
large, as seen in Fig. 1a: When the pulse duration is very small, the time
evolution is insufficient to drive the system far from the reference state. As
pulse duration is increased, the optimized energy error tends to decrease
exponentially (appearing linear on the semi-log plot), but slowly. However, at
some critical time, there is a sharp transition to highly accurate energies,
where the optimized energy error is on the order of numerical precision;
increasing pulse duration beyond this point has negligible impact on the
accuracy. The pulse duration for which this transition occurs is closely
related to the minimal evolution time (MET) from optimal control theory, which
is the absolute shortest time required to prepare one state from another,
subject to a given set of pulse constraints [17, 37]. According to optimal
control theory, amplitudes will tend to saturate their bounds as the pulse
duration approaches the MET. For real-valued amplitudes (the
$\quantity{\alpha}$ and $\quantity{\alpha\Delta}$ parameterizations), the MET
is achieved with bang-bang pulse shapes which are characterized as pulses
where the amplitude switches between $\pm\Omega_{\rm max}$, as observed in our
previous work [16] and in a control-theoretic analysis of the quantum
approximate optimization algorithm [38]. For complex-valued amplitudes,
optimal control theory imposes no a priori constraints on the phases $\phi$,
but the amplitude parameters $A$ will tend to saturate the bounds
$\pm\Omega_{\rm max}$. One can see this effect in Figs. 1c and 1d, which show
the amplitudes and phases optimized from zero at $T=21\mathrm{~{}ns}$, the
shortest pulse duration in Fig. 1a which results in a ctrl-VQE run obtaining
energy errors less than $10^{-8}\mathrm{~{}Ha}$ when all parameters are
initialized to zero.

Converging to the true MET is extremely computationally expensive, as it
requires, in principle, an exhaustive search over all parameter space. The
approach we have taken previously [16] simply samples over multiple
optimizations using randomly initialized pulses. For example, the gray ticks
we plot in Fig. 1a are the optimized energy errors from 100 runs when the
initial pulse parameters are uniformly sampled from the allowed range. Several
runs with $T=20\mathrm{~{}ns}$ are successful, and it may well be that a
perfect global optimization could locate successful pulses with even lower
duration. However, we choose to use the more computationally tractable zero-
pulse initialization as a heuristic for the remainder of this work, and we
will refer to the evolution time in an “error-vs-time” plot as the “effective”
minimal evolution time (eMET) $T_{0}$.

Interestingly, an “iterations-vs-time” plot (Fig. 1b) exhibits an iteration
count that is relatively low for most pulse durations, except for a broad and
large peak near the eMET $T_{0}$. This is due to the saturation of the
amplitude bounds observed above, straining the optimizer as it balances
reductions to the molecular energy and to the penalty term enforcing the
amplitude bounds. Therefore, we anticipate real ctrl-VQE experiments will want
to target pulse durations sufficiently beyond the eMET, to achieve the best
trade-off between the number of circuit evaluations and decoherence error.

###  IV.1 Characterizing eMET for different window lengths

(a)  (b)

Figure 2:  (2a) Ground state energy error vs. pulse duration of ctrl-VQE
applied to LiH with nuclear separation 3.0Å, using the
$\quantity{\alpha\beta}$ parameterization with three different numbers of time
windows. Each pulse is divided into as many uniformly spaced windows as
possible such that the window length $s$ is constrained as labeled. The
$\quantity{\alpha\beta}_{\infty}$ simulations use as many windows as there are
time-steps in the evolution, such that $s=0.05\mathrm{~{}ns}$. All parameters
are initialized to zero in the optimization steps. (2b) The number of BFGS
iterations needed to obtain the energies in (2a).

Throughout this paper, pulses are divided into a number of windows chosen such
that the duration of each window is as close as possible to, but never below,
some minimum duration $s_{\min}$. For example, if $s_{\min}=3\mathrm{~{}ns}$,
a $T=24\mathrm{~{}ns}$ pulse would be divided into eight windows of
$s=3\mathrm{~{}ns}$ each, while a $T=23.8\mathrm{~{}ns}$ pulse would be
divided into seven windows of $s=3.4\mathrm{~{}ns}$ each. This choice ensures
that longer pulses have more degrees of freedom while simultaneously dampening
high-frequency components in the control signals that are difficult to
implement experimentally. By choosing $s_{\min}=3\mathrm{~{}ns}$, we roughly
approximate an effective bandwidth of roughly $1/2s_{\min}\sim
167\mathrm{~{}MHz}$ in the microwave pulse generator.

Figure 2 shows the “error-vs-time” and “iterations-vs-time” plots comparing
different choices of $s_{\min}$ for a LiH molecule with a nuclear separation
of 3.0Å, using the parameter set $\quantity{\alpha\beta}$ (complex amplitudes
parameterized with Cartesian notation, and the drive frequency for each qubit
is fixed on resonance). The “continuous” curve
$\quantity{\alpha\beta}_{\infty}$ assigns independent $\alpha$ and $\beta$
parameters for each point in the Trotterized time evolution, corresponding to
a $s$ on the order of the Trotter step, $\sim 0.05\mathrm{~{}ns}$. Despite the
vastly larger number of parameters in this case compared to the other two
cases shown in the figure, the eMET $T_{0}$ changes by only about
$3\mathrm{~{}ns}$ or so, emphasizing the importance of evolution time over
degrees of freedom. As can be seen in the “iterations-vs-time” plot, the
larger number of parameters does have a significant effect on computational
cost, so we will adopt the reasonable value of $s_{\min}=3.0\mathrm{~{}ns}$
for most results in this paper.

We note that the somewhat higher asymptote observed for energy errors in Fig.
2a is an artifact of our optimization routine using the $L_{\infty}$ norm
rather than, say, the $L_{2}$ norm as the convergence criterion. Thus,
gradient vectors with a very large number of small but non-zero values will be
somewhat prematurely flagged as converged. We anticipate that a
correspondingly stricter choice of convergence criterion in the runs with more
parameters would result in a taller, broader, and more symmetric “iterations-
vs-time”curve than the one observed in Fig. 2b, though we have not confirmed
this numerically.

We can gain an intuitive understanding for why the number of parameters does
not have a significant impact on the eMET $T_{0}$ by studying the cartoon in
Fig. 3a: When the complex amplitude $\Omega$ varies continuously in time, the
optimized pulses generate the shortest possible path through Hilbert space
between the reference state and the ground state. The time it takes to
traverse this path is determined by the Quantum Speed Limit, which dictates
that the rate at which a system can move through Hilbert space is bounded by
the norm of the drive Hamiltonian, determined in our case by the drive
amplitudes $A$ [17]. Thus, hardware constraints on the maximal drive amplitude
$A$ for each pulse impose a fundamental MET $T_{0}$ required to prepare our
ground state. When we impose finite windows on the pulse shape, the optimal
path traced through Hilbert space is perturbed from the optimum. Longer
windows result in larger perturbations, and therefore longer evolution times.
However, as long as the window duration $s$ is relatively small with respect
to the pulse duration $T$, the dominant contribution to the evolution time
remains the original fundamental MET $T_{0}$.

Fig. 3b provides a particularly striking visualization of this perturbation.
First, we take two optimized pulses which successfully prepare the target
state, one with $s=0.05\mathrm{~{}ns}$ (“continuous”’) and
$T=17\mathrm{~{}ns}$, and one with $s=1.5\mathrm{~{}ns}$ and
$T=18\mathrm{~{}ns}$. We note that optimization of the $s=1.5\mathrm{~{}ns}$
case required randomly initializing pulse parameters to find a successful
pulse with $T=18.0\mathrm{~{}ns}$, as opposed to the data presented in Fig. 2,
which reports results when initializing all pulses from zero. Next, we
calculate the trajectories $\ket{\psi(t)}$, $\ket{\phi(t)}$ in the rotating
frame, when applying each pulse. Finally, for every pair of times $t$,
$t^{\prime}$, we evaluate the fidelity
$|\innerproduct{\psi(t)}{\phi(t^{\prime})}|^{2}$. The white spot in the bottom
left of Fig. 3b indicates that both pulses are applied to the same initial
state, while the white spot in the upper right indicates that both pulses
prepare the same final state after the full pulse duration. The bright band
roughly along the diagonal of the majority of the plot indicates that both
trajectories largely overlapped (except in the beginning), corroborating the
qualitative picture in Fig. 3a.

(a)  (b)

Figure 3:  (3a) A cartoon depiction of how finite window lengths affect
evolution time. When pulse parameters are allowed to vary continuously in
time, denoted by the black curve, ctrl-VQE can in principle find the shortest
possible path between the reference state and the ground state. When
parameters are constrained to be constant over finite windows, ctrl-VQE can
still approximate this ideal path. The total path length, and therefore the
total evolution time, will be marginally longer, but the main contribution to
evolution time is the length of the ideal path. (3b) The fidelity between
quantum states as a function of time, following trajectories solving LiH with
nuclear separation 3.0Å, from two optimized pulses: a $17.0\mathrm{~{}ns}$
pulse where parameters were varied quasi-continuously, and an
$18.0\mathrm{~{}ns}$ pulse where parameters were restricted to window lengths
of $s=1.5\mathrm{~{}ns}$ each. Bright colors indicate high overlap at the
corresponding points in the trajectory, while black indicates orthogonal
states.

###  IV.2 Polar vs Cartesian Parameterization

(a)  (b)

Figure 4:  (4a) Ground state energy error vs. pulse duration for ctrl-VQE
applied to LiH with nuclear separation 3.0Å, parameterizing each pulse’s
complex amplitude with polar $\quantity{A\phi}$ or Cartesian
$\quantity{\alpha\beta}$ notations. Each pulse is fixed on resonance and
divided into as many uniformly spaced windows as possible such that the window
length $s\geq 3.0\mathrm{~{}ns}$. All parameters are initialized to zero in
the optimization. (4b) The number of BFGS iterations needed to obtain the
energies in (4a).

In this section, we contrast ctrl-VQE using Cartesian notation (Eq. 8) and
polar notation (Eq. 7). Either choice leads to the exact same dynamics for any
single pulse, but a finite step in the different parameter sets will have
different resulting pulses, meaning a gradient-based optimizer like BFGS may
well proceed along different optimization trajectories. In Fig. 4, we show
“error-vs-time” and “iterations-vs-time” plots comparing the two different
notations applied to LiH with a bond separation of 3.0Å, holding frequencies
fixed on resonance. The energy plots and the eMET $T_{0}$ are essentially
identical, indicating that both runs identified pulses preparing essentially
the same estimate of the ground state for each pulse duration $T$. Note that
the actual pulse shapes do vary, but more than one pulse generates the same
state. However, the optimization trajectories for the polar notation are
notably more expensive, so we recommend adopting Cartesian notation. We assume
this arises due to larger off-diagonal matrix elements in the pulse Hessian of
the polar representation, but this is a topic we will explore in future work.

###  IV.3 Importance of $\phi$ and $\Delta$

(a)  (b)

Figure 5:  (5a) Ground state energy error vs. pulse duration for ctrl-VQE
applied to LiH with nuclear separation 3.0Å, for various parameterizations.
Each pulse for the $\quantity{\alpha\beta}$ and $\quantity{\alpha\beta\Delta}$
curves is divided into as many uniformly spaced windows as possible such that
the window length $s\geq 3.0\mathrm{~{}ns}$, while those for
$\quantity{\alpha}$ and $\quantity{\alpha\Delta}$ use $s\geq
1.5\mathrm{~{}ns}$, so that the total number of parameters scale similarly.
All parameters are initialized to zero in the optimization. (5b) The number of
BFGS iterations needed to obtain the energies in (5a).

Previous works on ctrl-VQE [15, 16] used the $\quantity{\alpha\Delta}$
parameterization, restricting drive phases to $\phi\in\quantity{0,\pi}$ and
allowing drive frequencies to vary off resonance. In this section, we will
explore what happens when the phase degree of freedom is included
($\quantity{\alpha\beta\Delta}$, in Cartesian notation), when the frequency
degree of freedom is removed ($\quantity{\alpha})$, and when both changes are
effected ($\quantity{\alpha\beta}$). Fig. 5 gives “error-vs-time” and
“iterations-vs-time” plots for all four parameterizations, applied to LiH with
a bond separation of 3.0Å. For the sake of fair comparison, pulses with real-
valued amplitudes are divided into twice as many windows, so that the total
number of parameters is similar for each value of $T$ (parameterizations
varying the frequency include one more parameter for each pulse, for a total
of four more parameters).

Unsurprisingly, the longest eMET $T_{0}$ is observed with the most restrictive
parameterization $\quantity{\alpha}$, and the second longest is the second-
most restrictive, $\quantity{\alpha\Delta}$. (Recall that varying $\Delta$ is
more restrictive than varying $\beta$, since each pulse has only one $\Delta$
parameter, irrespective of the number of windows.) Perhaps surprisingly, the
remaining two parameterizations $\quantity{\alpha\beta}$ and
$\quantity{\alpha\beta\Delta}$ are practically indistinguishable, except that
including the frequency degree of freedom puts a heavy strain on the optimizer
(this effect was also observed in Ref. [25]). Apparently, once the pulse
amplitude is allowed to become complex, there is no further variational
flexibility afforded by the frequency detuning.

We must emphasize that this is not related to the dimensionality of the
optimization problem. The number of degrees of freedom is similar (within a
constant) for all four parameterizations. We verify this numerically by
calculating the effective quantum dimension [39] of our trial state at several
steps of the optimization. This quantity is defined as the rank of the quantum
Fisher information [40], and it describes the number of directions in which
infinitesimal variations in our parameters can drive the system. When the
pulse has a sufficient number of parameters, we anticipate the effective
quantum dimension to saturate at the number of independent degrees of freedom
in the Hilbert space (30, for a four-qubit system). Computational basis states
are singular, isolated points with reduced dimension. Thus, we expect (and
find) initial pulses, with $\Omega=0$, to have an effective quantum dimension
of 16, irrespective of the number of parameters in the pulse. However, random
perturbations from a singular point are almost sure to result in a state with
full effective quantum dimension [39, 41], and indeed we find that the
effective quantum dimension is saturated at 30 for all subsequent trial states
when $T\geq 12\mathrm{~{}ns}$.

We now address why the frequency degree of freedom has no apparent impact on
the eMET, so long as we include the phase degree of freedom (i.e., contrasting
the $\quantity{\alpha\beta}$ and $\quantity{\alpha\beta\Delta}$
parameterizations). This can be partly explained by noting that the drive
phase $\phi$ and drive frequency $\nu$ are intrinsically related: the argument
to the trigonometric function in Eq. (5) is the sum $\nu t+\phi$. In the limit
where the phase $\phi$ is allowed to vary continuously, one could implement
perturbations to the drive frequency $\nu\rightarrow\nu+\Delta$ without
actually changing the drive frequency directly, by selecting
$\phi\rightarrow\phi+\Delta t$.

However, we observe essentially identical eMETs with and without a frequency
degree of freedom even for finite window lengths, where the phase $\phi$
cannot have a linear dependence on time. To obtain a better intuition for this
behavior, we will consider the analytical solution for a single qubit driven
by a pulse with a constant amplitude $A$, detuning $\Delta$, and phase $\phi$.
We will describe our quantum state using the Bloch sphere parameterization
$\quantity(\theta,\varphi)$:

$\ket{\psi}=\cos{\frac{\theta}{2}}\ket{0}+e^{i\varphi}\sin{\frac{\theta}{2}}\ket{1}.$
(13)

The objective in ctrl-VQE is to prepare a particular state from a reference
state; in this example, let us take the north pole $\quantity(0,0)$ as our
reference state ($\ket{\psi_{\mathrm{0}}}=\ket{0}$), and the target state
corresponds to the point $\quantity(\theta_{0},\varphi_{0})$.

Under the rotating wave approximation, the system evolves under the
Hamiltonian in Eq. (10) according to the time-dependent Schrödinger equation,
which can be solved exactly. This is the Rabi problem, which can be found in
standard quantum mechanics textbooks, e.g. see Ref. [42]. The Cartesian
coordinates of the state on a Bloch sphere in the rotating frame are:

$\displaystyle\expectationvalue{X}(t)$
$\displaystyle=x(t)\cdot\cos\phi+y(t)\cdot\sin\phi,$ (14)
$\displaystyle\expectationvalue{Y}(t)$
$\displaystyle=y(t)\cdot\cos\phi-x(t)\cdot\sin\phi,$ (15)
$\displaystyle\expectationvalue{Z}(t)$ $\displaystyle=\cos^{2}(\lambda
t/2)+\sin^{2}(\lambda t/2)\cdot\cos\theta_{L},$ (16)

where

$\theta_{L}=2\arctan(A/\Delta)$ (17)

is the largest attainable $\theta$ on the Bloch sphere, using the following
definitions:

$\displaystyle x(t)\equiv$ $\displaystyle-\sin(\theta_{L}/2)\cdot\sin(\Delta
t)\cdot\sin(\lambda t)$ $\displaystyle+\sin(\theta_{L})\cdot\cos(\Delta
t)\cdot\sin^{2}(\lambda t/2),$ (18) $\displaystyle y(t)\equiv$
$\displaystyle-\sin(\theta_{L}/2)\cdot\cos(\Delta t)\cdot\sin(\lambda t)$
$\displaystyle-\sin(\theta_{L})\cdot\sin(\Delta t)\cdot\sin^{2}(\lambda t/2),$
(19) $\displaystyle\lambda\equiv$ $\displaystyle\sqrt{A^{2}+\Delta^{2}}.$ (20)

The trajectory of the qubit state rotates around an axis oriented along a
polar angle $\theta_{L}/2$ and an azimuthal angle that begins along $\phi$ but
precesses slowly in time at a rate $\Delta$. In addition to setting
$\theta_{L}$, the amplitude $A$ also determines the speed at which this path
is traversed, so we will always want it as large as our hardware allows in
order to minimize the necessary evolution time.

(a)  (b)  (c)

Figure 6: Cartoon depictions of the path on a Bloch sphere traced by a single-
qubit pulse with constant $A$, $\phi$, and $\Delta$, when starting from the
north pole ($\ket{0}$). The path rotates around an axis with azimuthal
orientation starting at $\phi$ and precessing at a rate $\Delta$, and a
vertical tilt given by the ratio $A/\Delta$. When $\Delta=0$ (6a and 6b), the
axis is in the $x$-$y$ plane, so the path is a meridian with longitude given
by $\phi$.

When the detuning vanishes, $\Delta=0$, the axis is in the $x$-$y$ plane;
since we start from the north pole, the trajectory is a meridian, with the
longitude fixed by the pulse phase $\phi$. Therefore, the azimuthal coordinate
of our state $\varphi$ is constant (up to $\pm\pi$). The most direct state
preparation strategy is to choose $\phi=\varphi_{0}$ and drive on resonance
($\Delta=0$). Since the shortest path (i.e., the geodesic) between a pole and
any point on a sphere is along the meridian, this choice guarantees the
shortest possible evolution time, subject to constraints on the amplitude $A$.

If we restrict ourselves to real-valued pulses as in the previous section
(i.e., $\phi\in\quantity{0,\pi}$), the only azimuthal coordinates we can
access with resonant pulses are $\varphi\in\quantity{0,\pi}$ (see Fig. 6a). If
our target state has any other azimuthal coordinate $\varphi_{0}$, we must use
an off-resonant pulse (i.e., $\Delta\neq 0$). Doing so tilts the angle of the
axis about which our trajectory rotates; the trajectory is no longer a
geodesic. Thus, not only are large polar angles $\theta>\theta_{L}$ entirely
inaccessible, but target states with small $\theta<\theta_{L}$ will require a
longer path and therefore a longer pulse duration.

To summarize, reaching any arbitrary state requires varying either $\phi$ or
$\Delta$, but varying $\phi$ allows us to select an initial direction to move
in Hilbert space, while varying $\Delta$ allows us to introduce a curve in the
trajectory. In multi-qubit systems, the ideal trajectory is harder to
characterize, since the optimal trajectory using the control fields accessible
to us experimentally do not necessarily trace a geodesic in the full Hilbert
space. However, the freedom to adjust phases will tend to allow us to select a
more direct path, while adjusting frequency will tend to incur a longer path
length which takes longer to traverse.

##  V Conclusions

In this paper, we have explored the consequences of several different methods
to parameterize a pulse-level ansatz for VQE experiments. For each
parameterization, we analyzed how energy accuracy and optimization difficulty
vary over pulse duration. Here, we review our most salient observations and
summarize our main conclusions.

In every case, it is clear that ctrl-VQE requires a minimal pulse duration to
have any chance of successfully preparing the desired eigenstate. While
operating at this minimal time would be ideal for minimizing qubit
decoherence, and it would seem to minimize runtime, we find that the
optimization loop takes much longer to converge at and near this regime.
Furthermore, even locating the MET (or eMET) would be intractable in practice.
Consequently, a challenge arises for realistic simulations: how does one know
a priori that the specified evolution time in fact exceeds the MET, and thus
provides the possibility of preparing the target state? While answering this
definitively would lie beyond the scope of this paper, we note that the
quality of the final solution can be estimated by metrics such as variance
[43] or Hamiltonian-reconstruction distance [44] to determine whether the
optimization has located an eigenstate.

In Section IV.1, we explored the consequences of restricting pulse parameters
from continuous functions in time to discrete windows. So long as individual
windows are not too long, we found that ctrl-VQE provides essentially
consistent results, indicating that digitization of the analogue signals
should not have a significant impact on the performance. The minimal evolution
time required to obtain high-accuracy ansätze increases marginally as the
number of free parameters is reduced, but the far more significant difference
is that optimization is harder when the parameterization better approximates
continuously variable pulses. In some sense, this is a trivial statement:
optimization with fewer parameters is easier. Nevertheless, it has important
practical significance, since the runtime of ctrl-VQE is ultimately bounded by
the number of optimization iterations.

Similarly, in Section IV.2, we tested parameterizations for a complex
amplitude represented in both polar ($\Omega=Ae^{i\phi}$) and Cartesian
($\Omega=\alpha+i\beta$) coordinates. These two choices are physically
identical; the only difference between them is the trajectory through
parameter space which the optimizer selects to find optimal pulses. This is
analogous to the situation in classical molecular geometry optimization, where
different coordinate representations ($z$-matrix, redundant internals,
Cartesians, etc.) all specify the same molecule, but lead to widely varying
iteration counts for optimization. [45] We find that the polar
parameterization turns out to be significantly more expensive to converge,
especially in the region near the minimal evolution time. Therefore, we
recommend using Cartesian coordinates for experimental implementations of
ctrl-VQE.

Our most important results are found in Section IV.3, where we contrast
parameterizations which allow or disallow the phase and the frequency to vary
in the optimization. We find that the minimal evolution time is longest when
only amplitudes are varied, that it is improved when amplitudes and
frequencies are varied, and that it is improved even more when amplitudes and
phases are varied. Counter-intuitively, we find no measurable improvement to
the minimal evolution time when all three parameters are varied, although the
optimization is disproportionately harder. These results suggest that, in a
transmon-based ctrl-VQE, it is preferable to fix the drive frequencies of each
pulse (perhaps at the associated resonance frequency) and optimize over
amplitudes and phases.

Since single-qubit gates are often implemented with resonant pulses, it may be
surprising that resonant pulses are capable of building the entanglement
needed to represent the ground state of a molecular system. However, in an
architecture with static couplings, pulses applied locally in the “qubit”
basis will nevertheless have a global effect in the eigenbasis of the device.
In some sense, entanglement builds spontaneously, irrespective of pulse
frequency; the pulses serve only to direct it in the appropriate direction and
drive it along faster. Our findings imply that the phase of control fields is
the more effective parameter to provide that direction. Modulating this
parameter in time is also more natural in fixed-frequency, fixed-coupling
devices, such as those presently used by IBMQ [46]. An interesting direction
for future research is to investigate whether similar or complementary
heuristics can be found for alternate architectures, such as those with
tunable couplers, or those where the drive field is applied globally to all
qubits.

###### Acknowledgements.

The authors thank Chenxu Liu, Ayush Asthana, and Christopher Long for edifying
conversations. This research was supported by the U.S. Department of Energy
(DOE) Award No. DE-SC0019199, by the National Science Foundation (Grant No.
2137776). S.E.E. acknowledges support from the DOE Office of Science, National
Quantum Information Science Research Centers, Co-design Center for Quantum
Advantage (C2QA) under contract number DE-SC0012704.

##  VI Supplementary Information

Our latest implementation of ctrl-VQE can be found at [32].

## References

  * Bauer  _et al._ [2020] B. Bauer, S. Bravyi, M. Motta, and G. K.-L. Chan, Quantum Algorithms for Quantum Chemistry and Quantum Materials Science, Chemical Reviews 120, 12685 (2020), publisher: American Chemical Society.
  * McArdle  _et al._ [2020] S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan, Quantum computational chemistry, Reviews of Modern Physics 92, 015003 (2020).
  * Cao  _et al._ [2019] Y. Cao, J. Romero, J. P. Olson, M. Degroote, P. D. Johnson, M. Kieferová, I. D. Kivlichan, T. Menke, B. Peropadre, N. P. Sawaya, _et al._ , Quantum chemistry in the age of quantum computing, Chemical reviews 119, 10856 (2019).
  * Fedorov  _et al._ [2022] D. A. Fedorov, B. Peng, N. Govind, and Y. Alexeev, Vqe method: A short survey and recent developments, Materials Theory 6, 1 (2022).
  * Tilly  _et al._ [2022] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, The Variational Quantum Eigensolver: A review of methods and best practices, Physics Reports The Variational Quantum Eigensolver: a review of methods and best practices, 986, 1 (2022).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Variational quantum algorithms, Nature Reviews Physics , 1 (2021).
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, A variational eigenvalue solver on a photonic quantum processor, Nature communications 5, 1 (2014).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242 (2017).
  * Lee  _et al._ [2018] J. Lee, W. J. Huggins, M. Head-Gordon, and K. B. Whaley, Generalized unitary coupled cluster wave functions for quantum computation, Journal of chemical theory and computation 15, 311 (2018).
  * Ryabinkin  _et al._ [2018] I. G. Ryabinkin, T.-C. Yen, S. N. Genin, and A. F. Izmaylov, Qubit coupled cluster method: a systematic approach to quantum chemistry on a quantum computer, Journal of chemical theory and computation 14, 6317 (2018).
  * Gard  _et al._ [2020] B. T. Gard, L. Zhu, G. S. Barron, N. J. Mayhall, S. E. Economou, and E. Barnes, Efficient symmetry-preserving state preparation circuits for the variational quantum eigensolver algorithm, npj Quantum Information 6, 1 (2020).
  * Grimsley  _et al._ [2019] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature communications 10, 1 (2019).
  * Yordanov  _et al._ [2021] Y. S. Yordanov, V. Armaos, C. H. W. Barnes, and D. R. M. Arvidsson-Shukur, Qubit-excitation-based adaptive variational quantum eigensolver, Communications Physics 4, 1 (2021), publisher: Nature Publishing Group.
  * Tang  _et al._ [2021] H. L. Tang, V. Shkolnikov, G. S. Barron, H. R. Grimsley, N. J. Mayhall, E. Barnes, and S. E. Economou, qubit-adapt-vqe: An adaptive algorithm for constructing hardware-efficient ansätze on a quantum processor, PRX Quantum 2, 020310 (2021).
  * Meitei  _et al._ [2021] O. R. Meitei, B. T. Gard, G. S. Barron, D. P. Pappas, S. E. Economou, E. Barnes, and N. J. Mayhall, Gate-free state preparation for fast variational quantum eigensolver simulations, npj Quantum Information 7, 155 (2021).
  * Asthana  _et al._ [2023] A. Asthana, C. Liu, O. R. Meitei, S. E. Economou, E. Barnes, and N. J. Mayhall, Leakage Reduces Device Coherence Demands for Pulse-Level Molecular Simulations, Physical Review Applied 19, 064071 (2023), publisher: American Physical Society.
  * Deffner and Campbell [2017] S. Deffner and S. Campbell, Quantum speed limits: from heisenberg’s uncertainty principle to optimal quantum control, Journal of Physics A: Mathematical and Theoretical 50, 453001 (2017).
  * Liberzon [2011] D. Liberzon,  _Calculus of variations and optimal control theory_ (Princeton university press, 2011).
  * Ansel  _et al._ [2024] Q. Ansel, E. Dionis, F. Arrouas, B. Peaudecerf, S. Guérin, D. Guéry-Odelin, and D. Sugny, Introduction to theoretical and experimental aspects of quantum optimal control (2024), arXiv:2403.00532 [quant-ph].
  * Magann  _et al._ [2021] A. B. Magann, C. Arenz, M. D. Grace, T.-S. Ho, R. L. Kosut, J. R. McClean, H. A. Rabitz, and M. Sarovar, From pulses to circuits and back again: A quantum optimal control perspective on variational quantum algorithms, PRX Quantum 2, 010101 (2021).
  * Leng  _et al._ [2022] J. Leng, Y. Peng, Y.-L. Qiao, M. Lin, and X. Wu, Differentiable Analog Quantum Computing for Optimization and Control, in  _Advances in Neural Information Processing Systems_, Vol. 35, edited by S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (Curran Associates, Inc., 2022) pp. 4707–4721, arXiv: 2210.15812.
  * Meirom and Frankel [2023] D. Meirom and S. H. Frankel, PANSATZ: pulse-based ansatz for variational quantum algorithms, Frontiers in Quantum Science and Technology 2, 10.3389/frqst.2023.1273581 (2023), publisher: Frontiers.
  * Liang  _et al._ [2023] Z. Liang, J. Cheng, Z. Song, H. Ren, R. Yang, H. Wang, K. Liu, P. Kogge, T. Li, Y. Ding, and Y. Shi, Towards Advantages of Parameterized Quantum Pulses (2023), arXiv:2304.09253 [quant-ph].
  * Egger  _et al._ [2023] D. J. Egger, C. Capecci, B. Pokharel, P. K. Barkoutsos, L. E. Fischer, L. Guidoni, and I. Tavernelli, Pulse variational quantum eigensolver on cross-resonance-based hardware, Physical Review Research 5, 033159 (2023), publisher: American Physical Society.
  * Kottmann and Killoran [2023] K. Kottmann and N. Killoran, Evaluating analytic gradients of pulse programs on quantum computers (2023), arXiv:2309.16756 [quant-ph].
  * Keijzer  _et al._ [2023] R. d. Keijzer, O. Tse, and S. Kokkelmans, Pulse based Variational Quantum Optimal Control for hybrid quantum computing, Quantum 7, 908 (2023), publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften.
  * Entin  _et al._ [2024] N. Entin, M. M. Roses, R. Cohen, N. Katz, and A. Makmal, Molecular groundstate determination via short pulses on superconducting qubits (2024), arXiv:2403.17789 [quant-ph].
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, _et al._ , Pyscf: the python-based simulations of chemistry framework, Wiley Interdisciplinary Reviews: Computational Molecular Science 8, e1340 (2018).
  * Qiskit contributors [2023] Qiskit contributors, Qiskit: An open-source framework for quantum computing (2023).
  * Bravyi  _et al._ [2017] S. Bravyi, J. M. Gambetta, A. Mezzacapo, and K. Temme, Tapering off qubits to simulate fermionic hamiltonians (2017), arXiv:1701.08213 [quant-ph] .
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Review 59, 65 (2017), publisher: Society for Industrial and Applied Mathematics.
  * Sherbert and Mayhall [2023] K. M. Sherbert and N. J. Mayhall, Ctrlvqe.jl, https://github.com/mayhallgroup/CtrlVQE.jl (2023).
  * Krantz  _et al._ [2019] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A quantum engineer’s guide to superconducting qubits, Applied Physics Reviews 6, 10.1063/1.5089550 (2019).
  * Roth  _et al._ [2023] T. E. Roth, R. Ma, and W. C. Chew, The transmon qubit for electromagnetics engineers: An introduction, IEEE Antennas and Propagation Magazine 65, 8–20 (2023), conference Name: IEEE Antennas and Propagation Magazine.
  * Mogensen and Riseth [2018] P. K. Mogensen and A. N. Riseth, Optim: A mathematical optimization package for Julia, Journal of Open Source Software 3, 615 (2018).
  * Khaneja  _et al._ [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms, Journal of Magnetic Resonance 172, 296 (2005).
  * Wilhelm  _et al._ [2020] F. K. Wilhelm, S. Kirchhoff, S. Machnes, N. Wittler, and D. Sugny, An introduction into optimal control for quantum technologies (2020), arXiv:2003.10132 [quant-ph] .
  * Yang  _et al._ [2017] Z.-C. Yang, A. Rahmani, A. Shabani, H. Neven, and C. Chamon, Optimizing variational quantum algorithms using pontryagin’s minimum principle, Physical Review X 7, 021027 (2017).
  * Haug  _et al._ [2021] T. Haug, K. Bharti, and M. Kim, Capacity and Quantum Geometry of Parametrized Quantum Circuits, PRX Quantum 2, 040309 (2021), publisher: American Physical Society.
  * Stokes  _et al._ [2020] J. Stokes, J. Izaac, N. Killoran, and G. Carleo, Quantum Natural Gradient, Quantum 4, 269 (2020), publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften.
  * Funcke  _et al._ [2021] L. Funcke, T. Hartung, K. Jansen, S. Kühn, and P. Stornati, Dimensional expressivity analysis of parametric quantum circuits, Quantum 5, 422 (2021), publisher: Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften.
  * Townsend [2000] J. Townsend,  _A Modern Approach to Quantum Mechanics_, International series in pure and applied physics (University Science Books, 2000).
  * Zhang  _et al._ [2022] D.-B. Zhang, B.-L. Chen, Z.-H. Yuan, and T. Yin, Variational quantum eigensolvers by variance minimization, Chinese Physics B 31, 120301 (2022), publisher: Chinese Physical Society and IOP Publishing Ltd.
  * Moon  _et al._ [2024] L. J. I. Moon, M. M. Sohoni, M. A. Shimizu, P. Viswanathan, K. Zhang, E.-A. Kim, and P. L. McMahon, Hamiltonian-reconstruction distance as a success metric for the Variational Quantum Eigensolver (2024), arXiv:2403.11995 [quant-ph].
  * Baker  _et al._ [1996] J. Baker, A. Kessi, and B. Delley, The generation and use of delocalized internal coordinates in geometry optimization, The Journal of Chemical Physics 105, 192 (1996), https://pubs.aip.org/aip/jcp/article-pdf/105/1/192/19076673/192_1_online.pdf .
  * McKay  _et al._ [2017] D. C. McKay, C. J. Wood, S. Sheldon, J. M. Chow, and J. M. Gambetta, Efficient $Z$ gates for quantum computing, Physical Review A 96, 022330 (2017), publisher: American Physical Society.

##  Appendix A Evaluating Analytical Gradients with Complex Amplitudes

In this appendix, we derive an expression for the analytic gradient of the
expectation value $\expectationvalue{\hat{O}}{\psi}$ with respect to our
control parameters ${\bf\it\theta}$, where
$\ket{\psi}=\ket{\psi({\bf\it\theta})}$ is our fully-evolved state and
$\hat{O}$ is a Hermitian observable. This treatment is an adaptation of GRAPE
[36], tailored to our particular problem.

Formally, the objective function takes the form:

$E({\bf\it\theta})=\expectationvalue{\hat{O}}{\psi_{R}({\bf\it\theta})}$ (21)

By using a subscript $R$ in $\ket{\psi_{R}({\bf\it\theta})}$, we denote that
$\hat{O}$ is measured in the rotating frame.

The variational ansatz takes the form:

$\ket{\psi_{R}({\bf\it\theta})}=U_{R}({\bf\it\theta})\ket{\psi_{\mathrm{0}}}$
(22)

where $U_{R}({\bf\it\theta})$ is the unitary operator describing evolution
under the device Hamiltonian:

$U_{R}({\bf\it\theta})=\mathcal{\hat{T}}\exp(-i\int_{0}^{T}\mathop{}\\!\mathrm{d}t~{}\hat{V}_{R}({\bf\it\theta},t))$
(23)

The device Hamiltonian itself (in the rotating frame) takes the following
form:

$\hat{V}_{R}({\bf\it\theta},t)=e^{it\hat{H}_{0}}\hat{V}({\bf\it\theta},t)e^{-it\hat{H}_{0}}$
(24)

The static component $\hat{H}_{0}$ does not depend on our control parameters,
so it need not be specified further in this derivation. The drive component
$\hat{V}$ has the form:

$\hat{V}({\bf\it\theta},t)=\sum_{q}\hat{V}_{q}({\bf\it\theta},t)=\sum_{q}\quantity(\Omega_{q}({\bf\it\theta},t)e^{i\nu_{q}({\bf\it\theta},t)t}\hat{a}_{q}+\mathrm{h.c.})$
(25)

Recall that $\Omega_{q}$ is a complex number, encapsulating both the amplitude
and phase of the drive pulse. It will be helpful to recast the drive component
into a sum of Hermitian time-independent operators with real scalar time-
dependent coefficients:

$\hat{V}({\bf\it\theta},t)=\sum_{q}\quantity(x_{q}({\bf\it\theta},t)\hat{Q}_{q}+y_{q}({\bf\it\theta},t)\hat{P}_{q})$
(26)

Here $\hat{Q}_{q}\equiv\hat{a}_{q}+\hat{a}^{\dagger}_{q}$ and
$\hat{P}_{q}\equiv i(\hat{a}_{q}-\hat{a}^{\dagger}_{q})$ are the canonical
coordinate and momentum operators, equivalent to Pauli-$X$ and $Y$ spin
operators in the two-level truncation (except that our definition of $\hat{P}$
absorbs an extra negative sign for the sake of notational convenience). The
scalar functions $x_{q}$ and $y_{q}$ are the real and imaginary parts of
$\Omega_{q}\exp(i\nu_{q}t)$.

The simplest approach to characterizing time-dependent functions is to
discretize time. Therefore, let us define a time-interval $\tau\equiv T/r$,
where $r$ should be understood as arbitrarily large. We will abbreviate all
parametric time-dependent $f({\bf\it\theta},t)$ with the shorthand:

$f_{i}\equiv f({\bf\it\theta},i\tau)$ (27)

Now we choose to rewrite the integral in the exponential of $U_{R}$ as a
trapezoidal sum.

$\int_{0}^{T}\mathop{}\\!\mathrm{d}t~{}\hat{V}_{R}({\bf\it\theta},t)\approx\tau\sum_{j=1}^{r}\frac{\hat{V}_{R,j}+\hat{V}_{R,j-1}}{2}$
(28)

Furthermore, since $V_{R}$ is a conjugation of the drive component $\hat{V}$
with respect to a unitary frame rotation $\exp(-it\hat{H}_{0})$, this factor
can be brought outside the exponentiation $\exp(-it\hat{V}_{R})$:

$e^{-i\tau\hat{V}_{R,i}}=e^{it\hat{H}_{0}}e^{-i\tau\hat{V}_{i}}e^{-it\hat{H}_{0}}$
(29)

We may now expand $U_{R}$ as a time-ordered product. At the same time, we will
introduce left and right unitaries $U_{Ri}^{(l)}$, $U_{Ri}^{(r)}$
symmetrically about a single time point $t=i\tau$:

$\displaystyle U_{R}$ $\displaystyle\equiv
U_{Ri}^{(l)}U_{Ri}^{(r)}\text{~{}~{}~{}~{}~{}for any $i\in 0..r$}$ (30)
$\displaystyle U_{Ri}^{(l)}$ $\displaystyle\approx
e^{iT\hat{H}_{0}}\mathop{\overleftarrow{\prod}}_{j=i+1}^{r}\quantity(e^{-i\tau\hat{V}_{j}/2}e^{-i\tau\hat{H}_{0}}e^{-i\tau\hat{V}_{j-1}/2})$
(31) $\displaystyle U_{Ri}^{(r)}$
$\displaystyle\approx\mathop{\overleftarrow{\prod}}_{j=1}^{i}\quantity(e^{-i\tau\hat{V}_{j}/2}e^{-i\tau\hat{H}_{0}}e^{-i\tau\hat{V}_{j-1}/2})$
(32)

In order to take the gradient with respect to control parameters
${\bf\it\theta}$, let us first consider partial derivatives with respect to
$x_{qi}$ and $y_{qi}$. We will begin by writing
$\hat{V}_{i}=\hat{\Sigma}+x_{qi}\hat{Q}_{q}$, where $\hat{\Sigma}$ represents
all the terms in $\hat{V}_{i}$ except $x_{qi}\hat{Q}_{q}$. The analytical
expression for the derivative of the exponential of $\hat{V}_{i}$ is [36]:

$\partial
x_{qi}e^{-i\tau\quantity(\hat{\Sigma}+x_{qi}\hat{Q}_{q})}=e^{-i\tau\hat{\Sigma}}\int_{0}^{1}\mathop{}\\!\mathrm{d}s~{}e^{-i\tau
s\hat{\Sigma}}\quantity(-i\tau\hat{Q}_{q})e^{i\tau s\hat{\Sigma}}$ (33)

Because $\hat{Q}_{q}$ and $\hat{P}_{q}$ do not commute (in the two-level
truncation, the commutator is not even scalar), this has no closed-form
solution. However, so long as $\tau$ is very small, corrections due to the
commutator are negligible. Therefore, we choose to adopt the following much
simpler, approximate expressions:

$\displaystyle\partial_{x_{qi}}U_{Ri}^{(l)}$ $\displaystyle\approx
U_{Ri}^{(l)}\quantity(-i\tau\hat{Q}_{q}/2)$ (34)
$\displaystyle\partial_{x_{qi}}U_{Ri}^{(r)}$
$\displaystyle\approx\quantity(-i\tau\hat{Q}_{q}/2)U_{Ri}^{(r)}$ (35)
$\displaystyle\partial_{x_{qi}}U_{R}$ $\displaystyle\approx-i\tau
U_{Ri}^{(l)}\hat{Q}_{q}U_{Ri}^{(r)}$ (36)

Similarly,

$\partial_{y_{qi}}U_{R}\approx-i\tau U_{Ri}^{(l)}\hat{P}_{q}U_{Ri}^{(r)}$ (37)

Note that the symmetric form of our trapezoidal time integration serves to
further enhance the accuracy of the approximation.

We now define the “gradient signals” $\phi_{qi}^{(x)}$, $\phi_{qi}^{(y)}$:

$\displaystyle\phi_{qi}^{(x)}$
$\displaystyle=\imaginary\expectationvalue{U_{R}^{\dagger}\hat{O}U_{Ri}^{(l)}\hat{Q}_{q}U_{Ri}^{(r)}}{\psi_{\mathrm{0}}}$
(38) $\displaystyle\phi_{qi}^{(y)}$
$\displaystyle=\imaginary\expectationvalue{U_{R}^{\dagger}\hat{O}U_{Ri}^{(l)}\hat{P}_{q}U_{Ri}^{(r)}}{\psi_{\mathrm{0}}}$
(39)

These can be computed numerically in a time comparable to that of simulating
time-evolution itself. It is straightforward to show that the gradient of the
objective function $E$ with respect to $x$ and $y$ is:

$\displaystyle\partial_{x_{qi}}E$ $\displaystyle=2\tau\phi_{qi}^{(x)}$ (40)
$\displaystyle\partial_{y_{qi}}E$ $\displaystyle=2\tau\phi_{qi}^{(y)}$ (41)

Finally, a simple chain rule gives us the expression for the energy gradient
with respect to our control parameters ${\bf\it\theta}$:

$\partial_{\theta_{k}}E=2\tau\sum_{q,i}\quantity[\quantity(\partial_{\theta_{k}}x_{qi})\cdot\phi_{qi}^{(x)}+\quantity(\partial_{\theta_{k}}y_{qi})\cdot\phi_{qi}^{(y)}]$
(42)

We note that the above formula is in effect a discretized integral over time,
and accuracy can be tuned by adopting any quadrature rule (e.g., the
trapezoidal rule used in Eq. (28)). We conclude by explicitly writing the
integral formulation for continuous time:

$\displaystyle\partial_{\theta_{k}}E=2\sum_{q}\int_{0}^{T}\mathop{}\\!\mathrm{d}t~{}$
$\displaystyle\left[\frac{\partial
x_{q}({\bf\it\theta},t)}{\partial\theta_{k}}\cdot\phi_{q}^{(x)}(t)\right.$
$\displaystyle\left.+\frac{\partial
y_{q}({\bf\it\theta},t)}{\partial\theta_{k}}\cdot\phi_{q}^{(y)}(t)\right]$
(43)