††∗<EMAIL_ADDRESS>
# Approximations in transmon simulation

Tyler Jones1,2,3,∗ Kaiah Steven1 Xavier Poncini1 Matthew Rose1 Arkady
Fedorov2,3 1 Max Kelsen, Spring Hill, Queensland 4000, Australia
2 ARC Centre of Excellence for Engineered Quantum Systems, St Lucia,
Queensland 4072, Australia
3 School of Mathematics and Physics, University of Queensland, St Lucia,
Queensland 4072, Australia

###### Abstract

Classical simulations of time-dependent quantum systems are widely used in
quantum control research. In particular, these simulations are commonly used
to host iterative optimal control algorithms. This is convenient for
algorithms which are too onerous to run in the loop with current-day quantum
hardware, as well as for researchers without consistent access to hardware.
However, if the model used to represent the system is not selected carefully,
an optimised control protocol may be rendered futile when applied to hardware.
We present a series of models, ordered in a hierarchy of progressive
approximation, which appear in quantum control literature. The validity of
each model is characterised experimentally by designing and benchmarking
control protocols for an IBMQ cloud quantum device. This result demonstrates
error amplification induced by the application of a first-order perturbative
approximation. Furthermore, the emergence of errors which cannot be corrected
by simple amplitude scaling of control pulses is demonstrated in simulation,
due to an underlying mistreatment of non-computational dynamics. Finally, an
evaluation of simulated control dynamics reveals that despite the substantial
variance in numerical predictions across the proposed models, the complexity
of discovering local optimal control protocols appears invariant in the simple
control scheme setting.

###### pacs:

Valid PACS appear here

##  I Introduction

Since the turn of the century, superconducting qubits have emerged as a
leading candidate for the realisation of a universal fault-tolerant quantum
computer. To harness the potential of these systems, there remains an
important focus on improving the speed and accuracy of superconducting qubit
gate operations. This focus is particularly illustrated by recent work in the
field of quantum optimal control Allen  _et al._ (2017); Spiteri  _et al._
(2018); Werninghaus  _et al._ (2021); Kirchhoff  _et al._ (2018); Wu  _et al._
(2020); García-Ripoll  _et al._ (2020); Riaz  _et al._ (2019), an area
concerned with determining control protocols that minimise gate error.
Notably, state-of-the-art developments in machine learning have been adapted
to the quantum control domain An and Zhou (2019); Zhang  _et al._ (2019); Xu
_et al._ (2019); Wauters  _et al._ (2020); Zahedinejad  _et al._ (2016); Niu
_et al._ (2019); Bukov  _et al._ (2018); Wu  _et al._ (2019); Abdelhafez  _et
al._ (2020), leading to the discovery of control schemes that result in order-
of-magnitude reductions in gate infidelity Niu  _et al._ (2019). These
techniques are largely deployed in classical simulations rather than on
quantum hardware due to experimental constraints, including iteration speed
and difficulty in extracting the quantum state. Whilst the move to simulation-
based optimisation does carry significant benefits in accessibility and
iteration time, the models used have to rely on abstracted representations of
the relevant quantum system derived through a series of approximations. There
should be due caution taken when reporting optimal control results that
simplify system dynamics in such a way.

In this paper, we consider the system of a transmon qubit Koch  _et al._
(2007) coupled to a readout resonator in the dispersive regime. The simplicity
of this design, along with its charge-noise insensitive nature, has led to its
prevalence in research laboratories. The circuit architecture of the transmon
is nominally identical to the Cooper pair box architecture, with the
additional requirement of operation in a regime of high $E_{J}/E_{C}$
(Josephson energy and charging energy respectively). While this introduces
charge noise insensitivity, it also results in a weakened anharmonicity,
reducing isolation of the computational subspace from higher energy levels.
This enhanced risk of leakage intrinsically reduces gate operation speed, an
undesirable consequence considering the limited coherence times available in
modern devices. The dichotomy of suppressing leakage error and accelerating
operation speed has led to the discovery of several control schemes designed
to exploit nuances in system dynamics, most notably the derivative removal by
adiabatic gate (DRAG) pulse Motzoi  _et al._ (2009).

The problem of finding optimal quantum control schemes is often addressed
algorithmically, whether through standard gradient-based approaches like
Krotov, GRAPE or GOAT Reich  _et al._ (2012); Machnes  _et al._ (2018);
Khaneja  _et al._ (2004), or through the application of neural network-based
algorithms Niu  _et al._ (2019); Bukov  _et al._ (2018); Wu  _et al._ (2019).
A combination of the derivative computations and magnitude of iterations
required for these algorithmic approaches frequently limits them to operation
within simulations, rather than closed-loop optimisation. Given this, it is
imperative that relevant dynamics of the quantum system are not erased through
model abstraction.

This research seeks to examine the validity of a series of approximations made
to transmon-resonator systems, especially in the context of single qubit
control. Simulation environments will be developed for a defined hierarchy of
Hamiltonian models as seen in Fig. 1. These environments will be configurable
to match the energy parameters of any coupled transmon-resonator system,
permitting the simulation of IBMQ’s Qiskit Pulse compatible backends. This
will allow for the contrast of system properties at each level of
approximation, exposing the ramifications of each simplification. A set of
control parameters will be optimised to design a given gate using each
environment, and validated through deployment to the IBMQ device ibmq_armonk
to illustrate the enhanced fidelity induced by using the more faithful model.
Noting that the particular error type demonstrated on the IBMQ device can be
mitigated by parameter calibration in experiment, model disparities in the
frequency domain are presented in order to highlight a more pathological set
of deviations that are not able to be resolved through simple experimental
calibration. Finally, simulated control landscapes of each model will be
contrasted and traversed in order to quantify the ease of finding optimal
controls for each system. It is anticipated that these results will emphasise
the importance optimal control researchers should place on faithful system
modelling.

Figure 1: A hierarchy of system Hamiltonians, ordered by level of
approximation, is proposed. $H_{\mathrm{CPB}}$ represents the model which most
accurately describes a coupled transmon-resonator system. Following an arrow
indicates the application of an approximation which obfuscates the system
dynamics in some manner.

##  II  Transmon Theory

The design of the transmon qubit Koch  _et al._ (2007) is similar to the
Cooper pair box (CPB) qubit, containing a superconducting island consisting of
one or two111Using one Josephson junction results in a fixed frequency qubit,
while the usage of two in a dc superconducting quantum interference device
setup allows for the flux tuning of $E_{J}$. Both cases result in the same
Hamiltonian, provided junctions are identical. Josephson junctions. The
transmon qubit is distinguished by a further large capacitance shunting of the
Josephson junction, which alters the relative strengths of energy parameters
in the system. The system dynamics can be reduced to the isolated qubit CPB
Hamiltonian Koch  _et al._ (2007)

$\hat{H}=4E_{C}(\hat{n}-n_{g})^{2}-E_{J}\cos(\hat{\varphi}).$ (1)

Here, $E_{J}$ and $E_{C}$ represent the Josephson and charging energies
respectively, with $E_{C}=\frac{e^{2}}{2C_{\Sigma}}$ where $C_{\Sigma}$ is the
effective capacitance $C_{J}+C_{g}+C_{s}$ (see Fig. 2). Any effective offset
charge bias in the device is also included as $n_{g}$. This accounts for
unwanted charge noise arising from the environment, or for an external gate
voltage $V_{g}$, used for control purposes. The variables $\hat{n}$ and
$\hat{\varphi}$ are the canonical conjugate variables that represent the
Cooper pair number and gauge-invariant phase difference operators
respectively.

The choice of the ratio $E_{J}/E_{C}$ distinguishes different classes of
superconducting qubits; in the transmon the Josephson energy dominates the
charging energy, resulting in small charge dispersion. This in turn reduces
sensitivity to gate charge exponentially with $E_{J}/E_{C}$ Koch  _et al._
(2007), giving rise to longer coherence times. This is done at the cost of
reduced anharmonicity, although it has a much weaker dependence given as
$~{}(E_{J}/E_{C})^{-1/2}$.

Figure 2: Effective circuit diagram of the Cooper pair box. A Josephson
junction is placed in parallel to a capacitance $C_{s}$, and capacitively
coupled to a voltage source. In the transmon architecture, this Josephson
junction is further shunted by a large secondary capacitance.

For the effective operation of the transmon as a qubit, it is conventional to
couple it to a cavity resonator for readout purposes Koch  _et al._ (2007).
This can be achieved via a gate capacitance. In the event that the gate
capacitance is small compared to the transmon and resonator capacitances
($C_{g}\ll C_{r},C_{\Sigma}$), this system is well described by circuit QED,
with the effective Hamiltonian ($\hbar$ = 1) written as

$\displaystyle\hat{H}_{\mathrm{CPB}}=$ $\displaystyle
4E_{C}(\hat{n}-n_{g})^{2}-E_{J}\cos(\hat{\varphi})$ (2)
$\displaystyle+\omega_{r}\hat{a}^{{\dagger}}\hat{a}-2e\beta_{C}V_{rms}^{0}\hat{n}(\hat{a}+\hat{a}^{{\dagger}}),$

where $\hat{a}$ and $\hat{a}^{\dagger}$ are the canonical lowering and raising
operators for the resonator, $\beta_{C}=C_{g}/C_{\Sigma}$ is the voltage
divider ratio and $V^{0}_{rms}=\sqrt{\omega_{r}/2C_{r}}$.

In our discussion we take Eq. (2) to be the most faithful Hamiltonian, and use
it as the baseline model to make further comparisons. While it is possible to
solve Eq. (1) analytically in terms of Mathieu functions, an analytic solution
to the cavity coupled Hamiltonian Eq. (2) does not exist; instead, it is
solved through numerical diagonalisation. Given the reliance on this numerical
approach, it is useful to consider approximations of the underlying model to
reduce computational overhead of simulation. In this investigation, we place
our attention on three approximations of the underlying model present in the
literature; in increasing simplicity, these models will be referred to as the
sextic Duffing oscillator ($\mathrm{DO_{3}}$), the generalised quantum Rabi
(GR) model, and the two-level quantum Rabi (R) model.

The $\mathrm{DO_{3}}$ model constitutes the first level of approximation and
is derived by treating the transmon as a perturbed harmonic oscillator. This
approach is justified in the limit of $E_{J}/E_{C}\gg 1$, resulting in a weak
anharmonicity and conversely, eigenstates with highly localised phase,
$\Delta\hat{\varphi}\ll 1$. Consequently the cosine term can be treated as a
series expansion truncated at $K^{\mathrm{th}}$ order,

$\cos(\hat{\varphi})\approx\sum_{k=0}^{K}\frac{(-1)^{k}\hat{\varphi}^{2k}}{(2k)!},$
(3)

where at any finite $K$ the periodic boundary conditions are neglected. In the
limit of $E_{J}/E_{C}\gg 1$, terms beyond quadratic are treated as corrections
to the harmonic oscillator. It follows that the number and phase operator can
be expressed as

$\displaystyle\hat{n}$
$\displaystyle=\frac{i}{2\sqrt{\eta}}(\hat{b}^{\dagger}-\hat{b}),$ (4)
$\displaystyle\hat{\varphi}$
$\displaystyle=\sqrt{\eta}(\hat{b}^{\dagger}+\hat{b}),$ (5)

with $\eta=\sqrt{2E_{C}/E_{J}}$. For $K>1$, the uncoupled transmon system can
be written in the form of a Duffing oscillator,

$\displaystyle\hat{H}_{\mathrm{DO}_{K}}$
$\displaystyle=\sqrt{8E_{C}E_{J}}(\hat{b}^{\dagger}\hat{b}+\frac{1}{2})-E_{J}$
$\displaystyle-$ $\displaystyle
E_{J}\sum_{k=2}^{K}\left(\frac{2E_{C}}{E_{J}}\right)^{k/2}\frac{(-1)^{k}(\hat{b}^{{\dagger}}+\hat{b})^{2k}}{(2k)!},$
(6)

defining a family of Hamiltonians parametrised by the truncation order $K$.

With this approximation, resonator-induced coupling between transmon states
can be justifiably limited to nearest neighbours (see Appendix A). The coupled
system Hamiltonian is thus expressed in fundamental transmon literature Koch
_et al._ (2007) as the canonical Duffing oscillator model,

$\displaystyle\hat{H}_{\mathrm{DO_{2}}}=$
$\displaystyle\sqrt{8E_{C}E_{J}}(\hat{b}^{\dagger}\hat{b}+\frac{1}{2})-E_{J}-\frac{E_{C}}{12}(\hat{b}+\hat{b}^{\dagger})^{4}$
$\displaystyle+\omega_{r}\hat{a}^{{\dagger}}\hat{a}+ig(\hat{a}+\hat{a}^{\dagger})(\hat{b}^{\dagger}-\hat{b}).$
(7)

Note that the prefactors of the asymptotic number operator are absorbed into
$g$.

Preliminary analysis in an experimentally relevant $E_{J}/E_{C}$ regime with
the inclusion of a sextic term ($K=3$) proved to have non-negligible effects
on the system spectrum in comparison to the quartic case ($K=2$), causing us
to introduce the $\mathrm{DO_{3}}$ Hamiltonian,

$\displaystyle\hat{H}_{\mathrm{DO_{3}}}=$
$\displaystyle\sqrt{8E_{C}E_{J}}(\hat{b}^{\dagger}\hat{b}+\frac{1}{2})-E_{J}-\frac{E_{C}}{12}(\hat{b}+\hat{b}^{\dagger})^{4}$
$\displaystyle+\frac{E_{J}}{720}\left(\frac{2E_{C}}{E_{J}}\right)^{3/2}(\hat{b}+\hat{b}^{\dagger})^{6}$
$\displaystyle+\omega_{r}\hat{a}^{{\dagger}}\hat{a}+ig(\hat{a}+\hat{a}^{\dagger})(\hat{b}^{\dagger}-\hat{b}).$
(8)

The GR model can be derived from Eq. (II) with the inclusion of a perturbative
approximation. In the limit $E_{J}/E_{C}\gg 1$, the relative effect of the
sextic term in Eq. (II) becomes negligible and the quartic term can be
approximated by first-order eigenenergy corrections. In this case, the
Hamiltonian takes the form of a harmonic oscillator with modified
eigenenergies for each eigenstate $\ket{m}$,

$E_{m}=-\frac{E_{C}}{12}(6m^{2}+6m+3),$ (9)

which can be incorporated through the operator,

$\hat{E}_{m}=-\frac{E_{C}}{12}(6(\hat{b}^{\dagger}\hat{b})^{2}+6\hat{b}^{\dagger}\hat{b}+3).$
(10)

This correction yields the GR Hamiltonian

$\displaystyle\hat{H}_{\mathrm{GR}}$
$\displaystyle=\sqrt{8E_{C}E_{J}}(\hat{b}^{\dagger}\hat{b}+\frac{1}{2})-E_{J}$
$\displaystyle-\hat{E}_{m}+\omega_{r}\hat{a}^{{\dagger}}\hat{a}+g(\hat{a}+\hat{a}^{\dagger})(\hat{b}^{\dagger}+\hat{b}),$
(11)

where the unitary transformation $U=e^{-\frac{i\pi}{2}b^{\dagger}b}$ has been
applied to express the coupling term in its conventional form. The key
simplification made by this model is the omission of off-diagonal elements in
the transmon Hamiltonian. In this picture, the transition frequency between
the ground and first excited states can be recovered as
$\hbar\omega_{01}\approx\sqrt{8E_{C}E_{J}}-E_{C}$. This model of a transmon
system is used within quantum control literature Zahedinejad  _et al._ (2016);
Spiteri  _et al._ (2018); Allen  _et al._ (2017), often in conjunction with
the rotating wave approximation (RWA) Liebermann and Wilhelm (2016);
Daraeizadeh  _et al._ (2020); Machnes  _et al._ (2018); Werninghaus  _et al._
(2021); Basilewitsch  _et al._ (2019); Kirchhoff  _et al._ (2018); Wu  _et
al._ (2020); Niu  _et al._ (2019). Usage of the RWA induces negligible error
within the parameter regime investigated in this work (and thus is not
reported on), but can introduce significant error of its own in regimes of
strong coupling or strong driving Ballester  _et al._ (2012). The R model can
be obtained directly from the GR representation by neglecting anharmonic
terms, as seen in

$\displaystyle\hat{H}_{\mathrm{R}}$
$\displaystyle=(\sqrt{8E_{C}E_{J}}-E_{C})\frac{\hat{\sigma}_{z}}{2}$
$\displaystyle+\omega_{r}\hat{a}^{{\dagger}}\hat{a}+g\hat{\sigma}_{x}(\hat{a}+\hat{a}^{\dagger}),$
(12)

where the $\sigma_{i}$ denote the Pauli matrices. Omitting higher energy
levels is a negligent approach to transmon simulation. As such, the use of
this model is generally restricted to research that focuses on innovation
within the control algorithm at the expense of modelling a specific qubit
architecture Zhang  _et al._ (2019); An and Zhou (2019); Bao  _et al._ (2018).

###  II.1 Control

Transmon control can be realised through a capacitively coupled time-dependent
driving voltage $V_{g}(t)$. The driving Hamiltonian of interest is

$\hat{H}_{d}(t)=\hat{n}V_{g}(t).$ (13)

For models operating in the anharmonic oscillator basis
($\hat{H}_{\mathrm{DO_{3}}}$, $\hat{H}_{\mathrm{GR}}$, and
$\hat{H}_{\mathrm{R}}$), Eq. (13) can be re-expressed in terms of the operator
introduced in Eq. (4),

$\hat{H}_{d}=\frac{i}{2\sqrt{\eta}}(\hat{b}^{\dagger}-\hat{b})V_{g}(t).$ (14)

We format the drive pulse into a general sum of sinusoidal carrier frequencies
modulated by an envelope function,

$V_{g}(t)=\sum_{i}\Omega_{i}(t)\cos(\omega_{i}^{dr}t+\gamma_{i}),$ (15)

where $\omega_{i}^{dr}$ is the driving frequency, $\gamma_{i}$ is the offset
phase and $\Omega_{i}$ is the envelope function. In experimental situations,
the amplitude of the envelope function should be constrained to zero at $t=0$
and $t=T$, with $T$ representing the final time of the pulse. Additional
constraints are also introduced by the sampling rate of the hardware, which
constrains the minimum time under which an arbitrary envelope function may be
expressed.

Information leakage outside the computational subspace is a predominant issue
in transmon control. A common technique to combat leakage is to apply a two-
component drive pulse expressed as

$V_{g}(t)=\Omega(t)\cos(\omega^{dr}t)-\beta\dot{\Omega}(t)\sin(\omega^{dr}t),$
(16)

with parameter $\beta$ representing a coefficient that maximally suppresses
leakage when tuned. This technique is known as derivative removal by adiabatic
gate (DRAG) Motzoi  _et al._ (2009).

##  III  Methods

We have presented four models of a transmon qubit architecture. Taking the CPB
Hamiltonian $\hat{H}_{\mathrm{CPB}}$ as the baseline description of the
system, $\hat{H}_{\mathrm{DO_{3}}}$, $\hat{H}_{\mathrm{GR}}$, and
$\hat{H}_{\mathrm{R}}$ are approximations with decreasing levels of
sophistication. As detailed in Sec. II, each model is approximately equivalent
in the regime $E_{J}/E_{C}\gg 1$. Despite the prevalence of these
approximations in the literature, a comparative study of the associated models
with experimentally relevant parameters is lacking (although
$\hat{H}_{\mathrm{CPB}}$ has been well examined Willsch  _et al._ (2017)). We
seek to explore this area within the context of quantum control, where
correspondence between results of simulated protocols and protocols run on
hardware is of clear importance.

We initially explore the differences in behaviour evident in the spectrum of
each model, specifically the parameters relevant to quantum control: qubit
frequency and anharmonicity. The extent of spectral deviation observed between
individual models provides preliminary insight into the level of approximation
sufficient for simulation of a real transmon device.

As we are concerned with the validity of these assumptions in the context of
engineering high-fidelity pulse sequences, the dynamics of each model will be
evaluated when probed by identical Gaussian pulses. Each model will then be
used to design an optimal Gaussian $\frac{\pi}{2}$-pulse, by finely sweeping
over amplitudes for Gaussian pulses of constant timespan and observing
$\ket{0}$ and $\ket{1}$ populations. In order to generate pulses for
deployment on an experimentally relevant system, we match our simulation
environment parameters to that of a real experimental device, for which
purpose we select the remotely available IBMQ device ibmq_armonk. Device
parameters are obtained by a series of experiments conducted on the
ibmq_armonk system; details of these experiments are available in Appendix B.
Generated $\frac{\pi}{2}$-pulses will then be deployed on ibmq_armonk and
evaluated by randomised benchmarking procedures to establish gate fidelities
Magesan  _et al._ (2011). This will subsequently allow us to quantify the
losses when transferring to hardware introduced by each level of model
simplification.

In practice, the simplicity of the control parametrisation used in the above
experiment means that conventional amplitude scaling and frequency calibration
in laboratory conditions have the capacity to mitigate error induced by use of
a simplified model. We emphasise that these methods cannot be used in general
to correct approximation-induced errors. To do this, physical origins of the
differences between model Hamiltonians are described and numerical simulations
are performed to demonstrate the importance of the pulse frequency profile.

(a) Qubit frequency and $E_{J}/E_{C}$

(b) Anharmonicity and $E_{J}/E_{C}$

Figure 3: Spectral features of a coupled transmon-resonator system, under a
range of modelling schemes and $E_{J}/E_{C}$ ratios. Panel (a) demonstrates
significant qubit frequency deviations between modelling approaches for an
experimentally relevant regime ($n_{\mathrm{exp}}=1$) with $\omega_{01}$ in
the GR model held constant, while (b) demonstrates further deviations in the
qubit anharmonicity across models with $E_{C}$ (equivalent to $\alpha$ in the
GR model) held constant. In each case, the model spectra are seen to converge
as $E_{J}/E_{C}\rightarrow\infty$. Inset graphs illustrate how $E_{C}$ and
$E_{J}$ were individually varied to increase the $E_{J}/E_{C}$ ratio; in each
case, optimised to maintain the GR spectral feature of interest.

To extend this investigation, we finally seek to quantify the difficulty of
reaching an optimal control solution in each model, rather than contrasting
the solutions themselves. To this end, the Gaussian control landscape
evaluated by first excited level population and parametrised by pulse time and
amplitude will be compared for each model. In order to produce these
landscapes, the evolution of a ground-state transmon-resonator system is
simulated (with each of the four modelling approaches) over a sweep of each
control parameter.

The landscape investigation will then be extended to a DRAG pulse
parametrisation. The gradient optimisation technique GOAT Machnes  _et al._
(2018) will be used to generate a set of control trajectories, iteratively
determining an optimal Gaussian amplitude and tunable DRAG coefficient from a
distribution of starting points. With this tool, we seek to quantify the
complexity of navigation in a given control landscape more stringently,
introducing the metric $R_{\gamma}$ of a trajectory $\gamma$ as Nanduri  _et
al._ (2013)

$\displaystyle R_{\gamma}=\frac{d_{\gamma}}{d_{E_{\gamma}}},$ (17)

where $d_{\gamma}$ and $d_{E_{\gamma}}$ are the path length of $\gamma$ and
the Euclidean distance between the initial and final points of $\gamma$
respectively. As $d_{\gamma}\rightarrow d_{E_{\gamma}}$, the lower bound of
the $R_{\gamma}$-metric saturates and the control trajectory $\gamma$ is
considered optimal. Insight into the structural complexity can be established
by averaging the $R_{\gamma}$ metric over a set of gradient-optimised control
trajectories $\Gamma$. The distribution of starting points for these control
trajectories is chosen to be a multivariate Gaussian distribution centred
around an established optimum.

##  IV  Results

###  IV.1 System characterisation

In order to produce results relevant to current transmon-resonator systems,
system parameters residing within order-of-magnitude bounds of publicly
accessible IBMQ devices are selected, and expressed in Table 1.

$E_{C}/2\pi$ (GHz) |  $E_{J}/2\pi$ (GHz) |  $g/2\pi$ (GHz) |  $\omega_{r}/2\pi$ (GHz)
---|---|---|---
0.348 | 10.158 | 0.02 | 6.99
Table 1: Energy parameters for transmon-resonator simulation.

As is common in experimental transmon literature, the ratio of $E_{J}$/$E_{C}$
is of the $10^{1}-10^{2}$ order of magnitude. The frequent usage of the
assumption $E_{J}/E_{C}\gg 1$ in model derivations motivates investigation
into the effects of this energy scaling on system properties. Figure 3
illustrates the dependence of qubit frequency and anharmonicity on the
magnitude of $E_{J}/E_{C}$ for each model. Given the two degrees of freedom
that one can use to increase the ratio of $E_{J}/E_{C}$, we opt to perform
this operation in two distinct ways in Fig. 3: (a) constant Rabi model
frequency $\omega_{01}=\sqrt{8E_{J}E_{C}}-E_{C}$ and (b) constant Rabi model
anharmonicity $\alpha=-E_{C}$. Taking the previously defined simulation
parameters, $E_{J}^{\mathrm{exp}}/2\pi=10.158$ GHz and
$E_{C}^{\mathrm{exp}}/2\pi=0.348$ GHz, we define the parameter
$n_{\mathrm{exp}}$ to be the multiplier of the standard experimental ratio

$\displaystyle\frac{E_{J}}{E_{C}}=\left(\frac{E_{J}^{\mathrm{exp}}}{E_{C}^{\mathrm{exp}}}\right)n_{\mathrm{exp}}=29.19\>n_{\mathrm{exp}}.$
(18)

Figure 4: The dynamics induced in each model by a Gaussian pulse with period
142.2 ns, over a set of pulse amplitudes from 0 to 75 MHz. The first excited
level population is taken to be representative of the dynamics in each case.
There is a clear divergence in population between the family of models in red
(CPB) and orange ($\mathrm{DO_{3}}$), and the family of models in dark blue
(GR) and light blue (R). This result implies that the efficacy of optimised
control fields may depend markedly on the model used in the optimisation
process.

Figure 3 demonstrates significant spectral variance across models, with qubit
frequency and anharmonicity deviating from CPB to GR by over 15 and 50 MHz,
respectively, at $n_{\mathrm{exp}}=1$. To establish whether these model
inconsistencies affect system dynamics, each simulated system is probed by an
identical Gaussian pulse mixed with a carrier signal at their respective
frequencies. A pulse period of 142.2 ns is used, with amplitudes ranging from
0–75 MHz.

Examining Fig. 4, the population of the first excited level is plotted against
the amplitude of the probing Gaussian pulse for each of the four models in the
regime $n_{\mathrm{exp}}=1$. All are observed to exhibit periodic behaviour
interpreted as Rabi oscillations. As the amplitude is increased, the influence
of approximations on the dynamics of each system becomes apparent.

It can be observed that the first-order perturbative correction, employed for
the GR and R models but absent for the CPB and $\mathrm{DO_{3}}$ models,
induces a pronounced deviation in dynamics at later Rabi periods due to a
mismatch of Rabi frequencies. Appendix C presents an auxiliary result that
verifies that this difference is predominantly attributable to the first-order
perturbative correction rather than the neglection of the sextic term, both of
which occur in the transition from $\mathrm{DO_{3}}$ to GR. Meanwhile, the
omission of effects from outside the computational subspace (comparing GR and
R) leads to a small but perceptible variance in dynamics at high drive
amplitudes, while the basis change, cosine expansion and truncation implicit
in the comparison of CPB and $\mathrm{DO_{3}}$ have no apparent effect.

###  IV.2 Validation on hardware

To verify that the hierarchy of modelling approaches presented corresponds to
decreasing fidelity with a real system, and to examine the degree of these
differences, a validation procedure on hardware (ibmq_armonk) is conducted. To
do this, single qubit gates are designed in simulation for each model and
subsequently evaluated by a random benchmarking protocol on hardware.

In order to characterise the ibmq_armonk device hardware and design the
appropriate gates, a set of standard calibration experiments are performed to
approximate system parameters. These calibrations are detailed in Appendix B.
Energy parameters in simulation are then replaced with the values obtained via
these experiments.

Figure 5: Randomised benchmarking results for a single-qubit Clifford gate
set, transpiled using RX($\frac{\pi}{2}$) pulses. The pulse parameters used in
(a) and (b) were optimised using $\mathrm{DO_{3}}$ and GR models respectively.
The EPC (error per Clifford gate), established by exponential curve fitting,
increases by over a factor of 2 from $\mathrm{DO_{3}}$ to GR. This result
substantiates the earlier inference that the efficacy of an optimised control
field depends markedly on the model used in the optimisation process.

To universally perform any single-qubit gate, only RZ($\theta$) and
RX($\frac{\pi}{2}$) pulses require definition. To see this, consider the
single-qubit universal rotation gate

$\displaystyle U3(\theta,\phi,\lambda)$
$\displaystyle=\begin{pmatrix}\cos{\frac{\theta}{2}}&-e^{i\lambda}\sin{\frac{\theta}{2}}\\\
e^{i\phi}\sin{\frac{\theta}{2}}&e^{i(\phi+\lambda)}\cos{\frac{\theta}{2}}\end{pmatrix}$
$\displaystyle=RZ(\phi)RX(-\frac{\pi}{2})RZ(\theta)RX(\frac{\pi}{2})RZ(\lambda).$
(19)

The Qiskit application programming interface (API) provides optimised pulses
for these two gates, implementing RZ($\theta$) as a virtually executed change
of frame and RX($\frac{\pi}{2}$) as a 142.2 ns pulse. To evaluate the accuracy
of each proposed transmon-resonator model, we seek to replace the default RX
gate with a Gaussian pulse, described by optimal parameters established within
the corresponding simulation environment. In each case, the pulse period is
fixed (142.2 ns) and amplitude is optimised to produce the RX($\frac{\pi}{2}$)
gate of highest fidelity. The Qiskit-optimised default pulse is then replaced
for transpilation by the numerically optimised pulse, and randomised
benchmarking sequences are subsequently transpiled and evaluated with the
single-qubit Clifford gate set {H, X, Y, Z, S, S†}.

The numerical optimisation (via amplitude sweeping) of RX($\frac{\pi}{2}$)
pulses over 142.2 ns leads to amplitudes of 3.38 MHz for the CPB and
$\mathrm{DO_{3}}$ models, whilst amplitudes of 3.25 MHz are obtained for the
GR and R models, as seen in Table 2. Randomised benchmarking sequences
conducted with CPB- and $\mathrm{DO_{3}}$-optimised pulses produce error per
Clifford metrics of $1.98\times 10^{-2}$ in application to hardware.
Corresponding sequences conducted with GR- and R-optimised pulses produce
error per Clifford metrics of $4.11\times 10^{-2}$, an error magnification of
over 100% (as seen in Fig. 5). This substantiates the claim that simulated
control protocols optimised using Rabi models will yield inflated error rates
when directly migrated to hardware. As will be addressed in the upcoming
control landscape results, this divergence is likely to be amplified under a
stronger drive.

| $\hat{H}_{\mathrm{CPB}}$ | $\hat{H}_{\mathrm{DO_{3}}}$ | $\hat{H}_{\mathrm{GR}}$ | $\hat{H}_{\mathrm{R}}$
---|---|---|---|---
Amplitude (MHz) | 3.38 | 3.25
EPC | $1.98\times 10^{-2}$ | $4.11\times 10^{-2}$
Table 2: Optimised amplitudes for a 142.2 ns RX($\frac{\pi}{2}$) Gaussian
pulse, and error per Clifford (EPC) associated with each as evaluated by
randomised benchmarking on ibmq_armonk.

The error per Clifford metric approaches or exceeds $2\times 10^{-2}$ in all
cases. Improvements in this metric are available through two major areas. The
first is a refined model, through modelling noise and dissipation, and/or
obtaining access to precise system parameters, thereby removing reliance on
calibration experiments. A second pressing area is control pulse enhancements;
specifically, the usage of more sophisticated pulse parametrisations, or a
more efficient Clifford gate transpilation. This transpilation would be
achieved by using a basis set of fundamental pulses that can compile Clifford
gates using less gates on average than that of Eq. (19). Despite this scope in
improvement with respect to absolute fidelity, it is clear that the first-
order perturbation applied when moving to the Rabi family of models incurs a
damaging cost in simulation precision even for simple single-qubit controls.

It is important to note that in a laboratory environment, the calibration of
drive amplitudes for single-qubit control occurs with the system in-the-loop.
Indeed, the error type described above can be largely mitigated by a simple
drive amplitude calibration across models. We now seek to characterise a set
of control settings in which errors caused by these approximations cannot be
mitigated in experimental conditions.

###  IV.3 Detuned control

Throughout this section, experiment-like parameter calibrations are applied
within simulation to maximally match DO3 and GR dynamics in the computational
subspace. Specifically, we apply an amplitude scaling factor to each applied
pulse, calibrated such that Rabi oscillations are matched between each model
(mitigating the error seen in Fig. 4). Subsequently, we vary drive frequency
(for each pulse strength used) and observe the periodicity of Rabi
oscillations to approximate the Stark-shifted frequency associated with each
(note that differences in Stark shift are trivial). We perform these steps to
emphasise that the model-induced deviations are in general not resolvable
through simple parameter calibrations, and originate from the omission of off-
diagonal Hamiltonian elements and overall mistreatment of dynamics in the
noncomputational subspace.

The differences between the compositions of $\hat{H}_{\mathrm{DO}_{3}}$ and
$\hat{H}_{\mathrm{GR}}$ are minimal. Energy parameters EJ and EC can be tuned
to match qubit frequency and anharmonicity across transmon Hamiltonians, and
the drive Hamiltonians are by definition equivalent. Under these conditions,
the only material differences between the models are frequencies of higher-
level transitions such as $\ket{2}\rightarrow\ket{3}$ and so on, and the
existence of off-diagonal elements in the Hamiltonian for
$\hat{H}_{\mathrm{DO}_{3}}$ that provide couplings between distant energy
states such as $\ket{0}\rightarrow\ket{2}$ and $\ket{1}\rightarrow\ket{3}$.

Figure 6: Deviations in simulated control between the DO3 and GR models. The
$x$ axis is detuning of a 5 ns drive pulse from the qubit frequency. The $y$
axis is $E_{J}/E_{C}$ values for defining the GR Hamiltonian. $E_{J}$ and
$E_{C}$ values of a comparable ratio were chosen to precisely match the qubit
frequency and anharmonicity for the DO3 Hamiltonian. Amplitudes around 190 MHz
were chosen for the drive, scaled to maximally match dynamics. Both models
undergo evolution via 5 ns square pulses. The $z$ axis represents the maximum
infidelity experienced between models from a subset of times during the
evolution of the 5 ns square pulse, beginning in a computational superposition
state. As can be seen, maximum infidelities lie within the range of
frequencies associated with higher-level transitions, and exist across all
$E_{J}/E_{C}$ values. Infidelity is artificially increased at around
$E_{J}/E_{C}\sim 50$ due to the qubit frequency briefly nearing the resonator
frequency.

The extent to which leakage out of the computational subspace becomes
problematic for different energy regimes is illustrated in Fig. 6. Deviations
in predicted leakage fluctuate with two key properties; decreases in the
$E_{J}/E_{C}$ ratio, which is previously established as a measure of model
divergence, and the frequency profile of the drive pulse, which determines the
activations of energy level transitions and thereby the intrinsic model
differences in noncomputational dynamics. There exist clear infidelity maxima
in the negative drive detuning region. This corresponds to the range of drive
frequencies associated with higher-level and two-photon transitions. To
confirm that it is these transitions that give rise to the infidelities
observed, a cross-section of the plot taken at the indicative white line is
presented in Appendix D, with relevant higher-level transitions and two-photon
transitions characterised.

The inability to mitigate these errors through control parameter calibration
indicates that these are the major deviations of concern in practice.
Algorithmic searches for optimal quantum control schemes, in particular those
facilitated by neural network-based algorithms, yield pulse shapes that have
unintuitive profiles in both the frequency and time domains. It is evident in
Fig. 6 that designing any control pulse with non-negligible frequency
components within a sub-GHz detuning range relative to the qubit frequency
will risk significant and irreversible fidelity loss, without a measured
consideration of the model in use. This danger is most pertinent for
unorthodox tailored pulse shapes, but should also be taken into consideration
in simulation of existing gate implementations that exploit the
noncomputational subspace Chow  _et al._ (2013); Magesan and Gambetta (2020);
Negîrneac  _et al._ (2021).

A set of recommendations are synthesised for control practitioners as follows.

####  IV.3.1 Resonant pulse design

The term resonant is used here to denote pulse shapes tightly distributed
around the target transition frequency (e.g. low amplitude Gaussians). For a
single-qubit or few-qubit experiment operating within this control scheme, it
is likely that the deviations demonstrated in this work will have little
functional importance. The nature of the induced errors in this regime and the
small quantity of control parameters means that calibration in the loop is
feasible, expressly that pulse shapes can be straightforwardly scaled and
frequency-shifted to achieve maximal fidelities.

As systems scale up to many-qubit algorithms, the full characterisation of
individual components will become increasingly prohibitive and simulated
spectral predictions for leakage transitions will become an important factor.
Spectral crowding within neighbouring qubits (directly resulting in crosstalk
error) means that an accurate knowledge of leakage transitions within the
frequency domain becomes crucial to the design of control protocols. It is
clear in Fig. 6 that even with qubit frequency and anharmonicity values in the
GR model scaled to match those in the DO3 model, the dynamical disagreements
around leakage transitions are substantial enough to prohibit its use for such
a task.

We do not recommend great care when selecting a simulation model for resonant
single-qubit control protocols. In contrast, designing a protocol that
successfully navigates a crowded spectrum (e.g. many-qubit algorithms) while
minimising the activation of undesired transitions forbids the use of the
first-order perturbation approximation intrinsic to the general Rabi model, on
account of the infidelities observed in Fig. 6. This is a problem that plagues
this class of models throughout the range of experimentally relevant EJ/EC
ratios.

####  IV.3.2 Advanced pulse design

In the endeavour to increase overall operational fidelity of near-term
superconducting systems, high-dimensional multiparameter optimisation tools
for discovering efficient pulse protocols has blossomed into a popular field
of research Niu  _et al._ (2019); Bukov  _et al._ (2018); Wu  _et al._ (2019).
It is well established that a) the bulk of this optimisation is often
performed in simulation in order to facilitate large iteration volumes and b)
the pulses developed have nonstandard profiles in both the time and frequency
domains, which they are able to exploit to increase gate fidelities. In
addition to algorithmically designed controls, there are existing gate
implementations which make use of higher levels and ‘exotic’ transitions such
as two-photon transitions. It is these two audiences which we aim to reach
principally.

In this case, the recommendations are clear. Algorithmic approaches and exotic
gate implementations rely on the mitigation or exploitation of leakage
transitions, and to maximise the power and applicability of your protocol it
is crucial that the dynamics in these frequency regimes are modelled
accurately. We expect that the use of explicit higher-order nonlinear terms in
the transmon Hamiltonian (such as in $\hat{H}_{\mathrm{CPB}}$ and
$\hat{H}_{\mathrm{DO}_{3}}$) will help to avoid losses as simulated in Fig. 6
at all experimentally relevant EJ/EC ratios.

###  IV.4 Control landscapes

With our attention so far on understanding the regimes within which different
models will provide meaningfully different globally optimal control protocols,
we now seek to quantify an objective difficulty of determining this solution
within each model. In short, we seek to synthesise and demonstrate a method of
measuring ‘ease of control’. A general characterisation of this problem is
provided by quantum control landscapes, which are the subject of active
research Chakrabarti and Rabitz (2007); Moore  _et al._ (2008); Larocca  _et
al._ (2018); Day  _et al._ (2019); Kosut  _et al._ (2019). These structures
are defined as the manifold generated by mapping between control Hamiltonian
parameters and some control performance metric. A thorough analysis of control
landscapes requires an examination of both the topology (local and global
optima) and nontopological features (local structure) of the manifold. Pairing
these considerations will provide a basis from which the difficulty of
discovering optimal control protocols in each regime can be assessed.

(a) Gaussian control landscape ($\hat{H}_{\mathrm{CPB}}$)

(b) Deviations in Gaussian control landscape
($\lvert\hat{H}_{\mathrm{DO_{3}}}-\hat{H}_{\mathrm{GR}}\rvert$)

Figure 7: (a) The control landscape of the CPB model under application of a
Gaussian control pulse. The $z$ axis shows the population of the coupled
system state associated with computational $\ket{1}$, after being driven by a
Gaussian control pulse with parameters defined on the $x$ and $y$ axes. (b) A
visualisation of the deviations in simulated control between the
$\mathrm{DO_{3}}$ and GR models. The axes represent pulse parameters in the
same way as (a). The colour intensity of the figure is then the absolute
difference in first-level populations after application of the relevant
Gaussian control pulse. This figure illustrates the magnitude of model
deviation that can occur for high-energy control fields.

We realise the landscape objective function as the population of the
computational $\ket{1}$ state, to act as a probe of the system dynamics. As a
preliminary investigation, landscapes are constructed in a parameter space
defined by the amplitude and timespan of the Gaussian control pulse. Across
all models similar topological features are present, with landscapes
punctuated by ridges of high population interpreted as Rabi oscillations. This
basic phenomenon can be observed in Fig. 7(a), which presents the control
landscape under implementation of $\hat{H}_{\mathrm{CPB}}$. In the case of the
$\hat{H}_{\mathrm{R}}$ model, these ridges possess a degenerate structure and
consequently globally optimal protocols are accessible to naive gradient-based
optimisation algorithms. As for the remaining models, the first excited level
population is increased by travelling along a ridge in the direction of
decreasing pulse amplitude, conducive to reducing leakage. Contrasting
globally optimal points across the models $\hat{H}_{\mathrm{CPB}}$,
$\hat{H}_{\mathrm{DO_{3}}}$ and $\hat{H}_{\mathrm{GR}}$, we find that they
correspond to different ridges (Rabi oscillations). This implies that not only
will the initial conditions of a gradient-based optimisation determine whether
the process will achieve a globally optimal solution, but the optimal initial
conditions also vary across each model.

In Fig. 7(b), the absolute difference between control landscapes for the
$\mathrm{DO_{3}}$ and GR models are displayed. These models are chosen for
contrast due to the earlier identification that the transition between them
encompassed the most disruptive assumption (first-order perturbation
correction with $E_{J}/E_{C}\gg 1$). This figure, which demonstrates excited
state population deviations in excess of 70% within a realistic pulse
parametrisation space, acts as an extension of Fig. 4 in illustrating the
inaccuracies inherent in the generalised quantum Rabi model, although we again
stress that for this simple control parametrisation the error can be mitigated
by simple amplitude scaling.

Additionally, Fig. 7(b) lends further legitimacy to the randomised
benchmarking results presented in Fig. 5. Due to drive amplitude limits
imposed on the ibmq_armonk device used for benchmarking, the final optimised
gates had strengths of just $\sim$3 MHz. Locating this drive within the
parameter space in Fig. 7(b) sees a model deviation of just a few percent,
indicative of the fact that this extreme low-amplitude regime produces minimal
model deviations in comparison to other regions in parameter space. It can be
extrapolated that the deviations observed in randomised benchmarking results
would be intensified considerably on hardware without these drive amplitude
constraints.

With a simple Gaussian control protocol, we have seen that the modelling
approach impacts both the distribution of global optima and their
corresponding fidelity markedly. Elaborating on this initial analysis, we seek
to examine the control landscape of a DRAG pulse protocol. This
parametrisation removes the triviality of finding local optima in the Gaussian
control space (methodology of minimising amplitude at the cost of pulse time)
in favour of well-defined optima which maximally suppress leakage. Upon
identifying these optima, $10^{3}$ starting points are randomly sampled from a
Gaussian distribution around each. The gradient-based pulse optimisation
algorithm GOAT is then used to generate optimal control trajectories for each
point. In each trajectory, termination is achieved by meeting one of two
conditions: (1) first excited level population exceeding $1\\!-\\!(5\times
10^{-5})$, or (2) trajectory steps exceeding 100. Generally, trajectories are
terminated by condition (1). Terminations that occurred by condition (2), but
had their final point in a position where other trajectories had successfully
converged from, are assumed to converge along the path of other said
trajectories. In Table 3, the distribution of $R_{\gamma}$ metrics [defined in
Eq. (17)] over the set of control trajectories generated by GOAT is presented
for each model. The premise of the $R_{\gamma}$ metrics is detailed in Sec.
III, whereby the metric tends towards unity as trajectories approach
linearity.

$\hat{H}_{\mathrm{CPB}}$ | $\hat{H}_{\mathrm{DO_{3}}}$ | $\hat{H}_{\mathrm{GR}}$ | $\hat{H}_{\mathrm{R}}$
---|---|---|---
(1.185$\pm$0.165) | (1.190$\pm$0.157) | (1.188$\pm$0.178) | (1.167$\pm$0.175)
Table 3: Mean $R_{\gamma}$ metric ($\pm$ one standard deviation) statistics
over $10^{3}$ gradient-optimised DRAG protocols.

These results indicate that the complexity of landscape navigation is
effectively invariant across the models presented, at least under the DRAG
pulse parametrisation. This invariance demonstrates that local landscape
features are consistent, and can be combined with earlier results to determine
that optima location is the key symptom of model deviation. This conclusion
can be validated by observing the control trajectory endpoints produced by the
gradient-control algorithm. These endpoints are structurally identical for
each model, terminating on the surface of an elliptical parameter space that
had a large range in DRAG coefficient $\beta$-space and a small range in
Gaussian amplitude $\Omega$-space. This behaviour is demonstrated in Fig. 8,
omitting trajectories that began with initial points within the termination
parameter space.

Figure 8: The endpoint distribution of control trajectories generated by GOAT.
Initial points in parameter space sampled from a multivariate Gaussian
distribution around a previously known maximum are provided to the GOAT
algorithm for each model. The blue and red represent the trajectory endpoints
of the $\hat{H}_{\mathrm{CPB}}$ and $\hat{H}_{\mathrm{DO_{3}}}$ models
respectively, the distributions of which are near-identical in location and in
shape. Models $\hat{H}_{\mathrm{GR}}$ and $\hat{H}_{\mathrm{R}}$, in gold and
black, respectively, differ significantly in parameter-space location, but
share a similar elliptical distribution. These results demonstrate consistency
in the nature of solutions across models, but considerable differences in the
parameter-space regions at which they occur.

The resemblance in local landscape structure that is shared across models for
standard control techniques invites a future research direction for the
development of a model-transfer mapping function. If control landscapes for
simple conventional control protocols possess similar local curvature
features, it could be reasonably expected that an empirical mapping function
could reliably learn to translate a high-population location in parameter
space for a simplified model, to the equivalent location of a more accurate
model. In light of analysis in Sec. IV.3, within conventional resonant
(tightly distributed around the qubit frequency) control settings, this
proposed transfer function would approach a diagonal matrix scaling the
dimensions of control. Moving beyond resonant control, the mapping function
would acquire off-diagonal contributions. In this event a mapping function
could, for example, be parametrised by a neural network tasked with learning
the nontrivial transfer features across models. This would allow control
optimisation algorithms to operate in low-precision simulations and cast their
results to the same parameter space in more accurate models, aiming to provide
a reduction in computational resource requirements of up to 16 times as
discussed in Appendix F. Such a solution could conceivably be most valuable in
accelerating a many-qubit control simulation while mitigating crosstalk from a
crowded spectrum.

##  V Conclusion

We have presented and characterised four relevant and distinct models of a
coupled transmon-resonator architecture. Examining their respective dynamics
in a single-qubit control scenario, we have concluded that the reduction in
required computational resources of the commonly used generalised quantum Rabi
model $\hat{H}_{\mathrm{GR}}$ and two-level quantum Rabi model
$\hat{H}_{\mathrm{R}}$ comes at the cost of significant accuracy loss in
hardware simulation. We have demonstrated loss in hardware fidelity across
models through randomised benchmarking processes on a cloud quantum device for
a simple qubit control scheme. We have provided insight into the origins of
these deviations, presented simulations of a particularly pathological class
of errors, and produced recommendations to optimal control practitioners based
on these findings. We have also characterised the models in the context of
standard control parameter spaces, presenting both an analysis methodology for
and insight into the difficulty of optimal control for each.

Contemporary superconducting quantum devices operate in regimes where the
validity of $E_{J}/E_{C}\gg 1$ is ambiguous, naturally obscuring the optimal
choice of effective Hamiltonian. This investigation provides guidelines for
the use of a common set of Hamiltonians, illustrates the loss in fidelity
which can be experienced by being liberal with the application of common
assumptions, and provides insight into the nature of deviations between these
models, all within a set of energy parameters which are experimentally
relevant for noisy intermediate-scale quantum devices.

Burgeoning techniques for precise quantum control rely on software which
simulates the unitary evolution of a quantum system as accurately as possible,
with auxiliary channels responsible for noise and other non-unitary processes
Youssry  _et al._ (2020). The deployment of learnings from this investigation
will maximise the overall performance of these techniques. Further natural
extensions of this work arise in more stringently evaluating the validity of
these results for the multiqubit case, and investigating the efficacy of the
proposed model-transfer mapping function.

##  Acknowledgements

We acknowledge the IBM Quantum team for helpful discussions, and access to
devices for the experimental portion of this research. We also acknowledge the
team behind QuTip, a package used extensively for the simulation portion of
this research. We thank the anonymous referees responsible for providing
valuable input on this manuscript. Finally, we thank Dr. Maciej Trzaskowski
and the wider management at Max Kelsen for support and direction.

##  Appendix A Nearest-neighbour coupling

With this approximation, the cavity coupling term can be re-expressed by
calculating matrix elements of the number operator in the basis of the
uncoupled transmon eigenstates $\ket{i}$,

$\hat{n}^{(t)}_{ij}=\bra{j}\hat{n}^{(c)}\ket{i}.$ (20)

The superscripts $t$ and $c$ denote the transmon basis and charge basis
respectively. In the asymptotic limit of large $E_{J}/E_{C}$, perturbation
theory for the matrix elements yields,

$\displaystyle\bra{j+1}\hat{n}\ket{j}$
$\displaystyle\approx\sqrt{j+1}\left(\frac{1}{2\sqrt{\eta}}\right)$ (21)
$\displaystyle\bra{j+k}\hat{n}\ket{j}$ $\displaystyle\approx 0,$ (22)

where $k>1$. From this, we make the approximation that resonator-induced
coupling between transmon states is limited to nearest neighbours.

##  Appendix B Hardware characterisation

The design of optimal Gaussian pulses for deployment on a real device requires
careful characterisation of said device. We utilise Qiskit Pulse access to the
IBMQ backend ibmq_armonk to obtain these parameters. A significant volume of
the code used in this process is adapted from the IBM Qiskit textbook Asfaw
_et al._ (2020).

###  B1 Qubit spectroscopy

IBMQ backends provide regularly updated estimates of qubit frequencies through
their API; nonetheless, the characterisation approach we take includes a qubit
spectroscopy process at the beginning of each calibration to ensure the
integrity of the subsequent steps. To perform this process, Gaussian pulses
are individually mixed with carrier signals at 50 frequencies evenly
distributed in a 50 MHz region around the estimated qubit frequency and
applied to the qubit. A Lorentzian distribution is fitted to the resultant
transmission data to approximate the point of peak transmission (qubit
frequency).

###  B2 $\pi$-pulse calibration

For further characterisation, a $\pi$-pulse that can be used to reliably flip
between computational states is required. For this calibration, a carrier
signal is fixed at the previously established qubit frequency, whilst a
constant-time (of the order 10-1 $\mu$s) Gaussian was swept over a range of
pulse amplitudes (on the order of MHz). A sinusoidal fit to the resultant
transmission data yielded a Rabi oscillation period, half of which can be
adopted as an optimal $\pi$-pulse amplitude.

###  B3 $\ket{e}-\ket{f}$ spectroscopy

In order to find the second excited level transition frequency, one can apply
the same principle as qubit spectroscopy with an initial state of the first
excited level $\ket{e}$. At the time of experimentation, the configuration of
IBM’s Qiskit Pulse package was such that only a single frequency could be
associated with a drive line, which prohibited the intuitive protocol of a
$\pi$-pulse followed by a sweep of carrier frequencies (mixed with Gaussian
pulses) over the frequency region in which the second excited transition is
expected. As such, a method of controlling the frequency of the signal
indirectly through the drive pulse envelope is required.

To do this, the second Gaussian drive pulse (to excite the higher level
transition) is further mixed with the sinusoid $e^{2i\pi\Delta ft}$ to
generate sidebands at $f_{c}\pm\Delta f$ (where $f_{c}$ is the frequency of
the carrier signal). Mixing a carrier signal at the qubit frequency with a
$\pi$-pulse followed by a $\Delta f$ sweep of sideband pulses allows the
second excited transition frequency to be obtained through a Lorentzian fit of
transmission data, as in qubit spectroscopy.

###  B4 Dispersive shift

The dispersive shift of a transmon-resonator system materialises as half the
difference between the frequencies of the $\ket{g}$ and $\ket{e}$ transmission
peaks. As such, a $\frac{\pi}{2}$-pulse is used to prepare the transmon in a
superposition state, and the measurement pulse is swept over a sufficiently
fine range of frequencies in order to resolve the two separate transmission
peaks apparent at $\omega_{qb}-\chi$ and $\omega_{qb}+\chi$ (corresponding to
$\ket{g}$ and $\ket{e}$, respectively). This resolution is achieved for the
transmon system by sweeping the resonator frequency over 50 points within a
500 kHz range around the raw resonator frequency.

###  B5 System parameters

Energy parameters $E_{J}$, $E_{C}$, $g$ and $\omega_{r}$ are required for the
modelling approaches detailed in this article. Resonator frequency
$\omega_{r}$ is provided by the Qiskit API. Methods for system parameter
estimation vary depending on the model in use.

For the Rabi models, the process is simplified by the employment of relations
that can be derived from the Rabi set of assumptions. Transmon energy
parameters $E_{C}$ and $E_{J}$ are estimated from the above measurements using
Rabi relations
$\omega_{\ket{e}\rightarrow\ket{f}}-\omega_{\ket{g}\rightarrow\ket{e}}=E_{C}$
and $\omega_{\ket{g}\rightarrow\ket{e}}=\sqrt{8E_{C}E_{J}}-E_{C}$. The
coupling strength $g$ is estimated by the dispersive regime relation
$\chi=-g^{2}\frac{E_{C}}{\Delta(\Delta-E_{C})}$.

For the CPB and $\mathrm{DO_{3}}$ models, solutions for $E_{J}$, $E_{C}$, and
$g$ that return the measured spectrum and dispersive shift are calculated
numerically, beginning in the neighbourhood of the parameters estimated by
using Rabi relations.

##  Appendix C Perturbation or sextic omission?

Results in Sec. IV.1 highlight deviations in model predictions that manifest
in the transition between $\hat{H}_{\mathrm{DO_{3}}}$ and
$\hat{H}_{\mathrm{GR}}$. This divergence presents itself when comparing the
dynamics of each model under a Gaussian drive of period 142.2 ns. The
procedure to obtain $\hat{H}_{\mathrm{GR}}$ from $\hat{H}_{\mathrm{DO_{3}}}$
relies on neglecting the higher-order sextic term in the latter model, and
subsequently applying a first-order perturbative approximation to the
associated quartic term. To verify that the model deviations are largely
explainable by the perturbative approximation, the model
$\hat{H}_{\mathrm{GR_{3}}}$ is introduced. This Hamiltonian replicates
$\hat{H}_{\mathrm{GR}}$, while additionally including leading-order
corrections for the sextic term from $\hat{H}_{\mathrm{DO_{3}}}$, expressed as

$\displaystyle\hat{H}_{\mathrm{GR_{3}}}$
$\displaystyle=\sqrt{8E_{C}E_{J}}(\hat{b}^{\dagger}\hat{b}+\frac{1}{2})-E_{J}$
$\displaystyle-\hat{E}_{m_{3}}+\omega_{r}\hat{a}^{{\dagger}}\hat{a}+gi(\hat{a}+\hat{a}^{\dagger})(\hat{b}^{\dagger}-\hat{b}),$
(23)

where in contrast to Eq. (10), the eigenenergy corrections introduced by
$\hat{E}_{m_{3}}$ include third-order terms,

$\displaystyle\hat{E}_{m_{3}}$
$\displaystyle=\frac{E_{J}}{720}\left(\frac{2E_{C}}{E_{J}}\right)^{3/2}(20(\hat{b}^{\dagger}\hat{b})^{3}+30(\hat{b}^{\dagger}\hat{b})^{2}+40\hat{b}^{\dagger}\hat{b}+15)$
$\displaystyle-\frac{E_{C}}{12}(6(\hat{b}^{\dagger}\hat{b})^{2}+6\hat{b}^{\dagger}\hat{b}+3).$
(24)

The $\mathrm{DO_{3}}$, $\mathrm{GR_{3}}$ and $\mathrm{GR}$ models are each
probed by an identical Gaussian pulse mixed with a carrier signal at their
respective frequencies, as in Sec. IV.1. The driving pulse uses a period of
142.2 ns, with amplitudes ranging from 0–75 MHz.

To underscore the model deviations between these intermediate models, we
introduce the variable

$\displaystyle\Delta_{\mathrm{M_{2}}}^{\mathrm{M_{1}}}$
$\displaystyle=P_{\mathrm{M_{1}}}-P_{\mathrm{M_{2}}},$ (25)

where $P_{\mathrm{M_{1}}}$ and $P_{\mathrm{M_{2}}}$ denote the first excited
level populations as calculated in the models M1 and M2. The differences
$\Delta_{\mathrm{M_{2}}}^{\mathrm{M_{1}}}$ quantify the influence of the
associated approximations on the system dynamics.

Figure 9: The behaviour of the metric
$\Delta_{\mathrm{M_{2}}}^{\mathrm{M_{1}}}$ under application of a 142.2 ns
Gaussian driving pulse, plotted with respect to the amplitude of said driving
pulse. This metric measures the difference in simulated first energy level
population between models M1 and M2. In blue, substantial differences between
models $\mathrm{GR_{3}}$ and $\mathrm{DO_{3}}$ can be observed, with the
difference between these models being the application of a first-order
perturbative approximation. In comparison, the red points display the
relatively small differences between models $\mathrm{GR_{3}}$ and
$\mathrm{GR}$, attributable to omission of the sextic perturbative term in
simulation. This justifies the statement that the first-order perturbative
correction in deriving $\hat{H}_{\mathrm{GR}}$ from
$\hat{H}_{\mathrm{DO_{3}}}$ is responsible for the magnitude of divergence in
simulated dynamics between these models.

The behaviours of $\Delta_{\mathrm{DO_{3}}}^{\mathrm{GR_{3}}}$ and
$\Delta_{\mathrm{GR}}^{\mathrm{GR_{3}}}$ across the defined range of driving
protocols are demonstrated in Fig. 9. The effects of neglecting the sextic
term are trivial in comparison to applying the first-order perturbative
approximation, illustrated in red and blue respectively. The validity of the
perturbative approximation (predicated on $E_{J}/E_{C}\gg 1$) is therefore
taken to be the dominant assumption when obtaining $\hat{H}_{\mathrm{GR}}$
from $\hat{H}_{\mathrm{DO_{3}}}$.

##  Appendix D Origins of deviations in the frequency domain

The origins of divergences between the DO3 and GR models introduced in Fig. 6
are elucidated in Fig. 10. The parameter calibration in amplitude and
frequency leads to a significant mitigation in error at the
$\ket{0}\rightarrow\ket{1}$ transition, where maximum infidelity is less than
1% over the duration of a 5 ns pulse. In comparison, the two-photon transition
at $\ket{0}\rightarrow\ket{2}$ increases error to around 1.2%. Finally, the
frequency region that comprises a set of noncomputational transitions
routinely sees infidelities in excess of 10%, most notably at
$\ket{1}\rightarrow\ket{3}$ where a significant portion of the original
superposition state is actively driven to the third excited level. This result
acts to demonstrate that appropriate amplitude scaling and Stark shift
corrections are insufficient measures to ensure a viable control solution, for
a control pulse with an arbitrary frequency profile.

Additionally, it is valuable to note that the off-diagonal matrix elements in
$\hat{H}_{\mathrm{DO_{3}}}$ that are responsible for the dynamical deviations
at transitions such as $\ket{0}\rightarrow\ket{2}$ tend towards zero in the
limit of $E_{J}/E_{C}\rightarrow\infty$. This is a key assumption that allows
the first-order perturbation approximation inherent to
$\hat{H}_{\mathrm{GR}}$. The clear discrepancies in these models at
$E_{J}/E_{C}\sim 100$, a relatively large experimental ratio, further
highlights the questionable justification of this approximation for standard
experimental energy regimes.

Figure 10: Deviations in simulated control between frequency-matched and
amplitude-scaled DO3 and GR models at $E_{J}/E_{C}\sim 100$. The $x$ axis is
detuning of a 5 ns square drive pulse from the qubit frequency. Amplitudes
around 190 MHz are chosen for the drive, scaled to maximally match dynamics.
The system begins in a 50-50 superposition between the ground and excited
states. The $y$ axis represents the maximum infidelity experienced between
models from a subset of times during the evolution of the 5 ns square pulse.
Significant infidelity peaks exist at frequencies that correspond to higher-
level transitions and two-photon transitions, as well as some small deviation
at the qubit frequency. In particular, comparison of the models at the two-
photon transition $\ket{1}\rightarrow\ket{3}$ yields an infidelity in excess
of 50%.

##  Appendix E Numerical simulation

All model simulations (propagations of time-dependent Hamiltonians) are
performed using the zvode ordinary differential equation (ODE) solver provided
in Python package Qutip Johansson  _et al._ (2012). The solver is deployed
with an absolute tolerance of $10^{-8}$, relative tolerance of $10^{-6}$, and
$5\times 10^{3}$ time steps.

The accuracy of numerical simulation of quantum systems can be inhibited by
the number of energy levels used. To determine the required Hilbert space to
accurately simulate dynamics for each model, we drive an $N$-level transmon
from $\ket{0}$ to $\ket{\psi_{F}}$ with a standard Gaussian pulse of
experimentally relevant parameters, from $N=2$ to $N=20$ until an error
tolerance criteria is met. This criteria is the mean absolute population
difference in the computational subspace between $n$ and $n-1$ across the
duration of the pulse $T$ (which consists of $f$ timesteps, [$t_{0}=0$,
$t_{f}=T$]). This can be expressed as

$d\left(p_{m}(n,t)\right)=\frac{1}{f}\sum_{j=0}^{f}|p_{m}(n,t=t_{j})-p_{m}(n-1,t=t_{j})|,$
(26)

where $p_{m}(n,t)$ returns a vector of the model m populations with $n$ energy
levels at time step $t$. The model is deemed to converge when this metric
reaches 10-5. For $\hat{H}_{\mathrm{CPB}}$, $\hat{H}_{\mathrm{DO_{3}}}$ and
$\hat{H}_{\mathrm{GR}}$, required Hilbert space dimensionality was found to be
13, 12, and 6, respectively. Model $\hat{H}_{\mathrm{R}}$ requires only two
energy levels by definition.

##  Appendix F Computational requirements

Given the iterative optimisation protocols often used for quantum control, a
balance must be struck between the model fidelity and computational footprint.
As resource requirements will be largely constrained by the model Hilbert
space, the approach detailed in Appendix E is used to define the number of
energy levels required for each model. To examine the computational resource
requirements of each model, each simulated transmon is driven at its
transition frequency by a Gaussian pulse of amplitude 75 MHz, with period
ranging between 0 and 150 ns. Linearly fitting these results provides
characteristic runtime statistics for each model, which in conjunction with
results presented in the body of this paper yields a clear picture of the
balance between accuracy and resource requirement. Figure 11 illustrates these
runtime statistics, using a zvode-centric ODE solver.

Figure 11: Normalised simulation times for Gaussian control pulses of varying
periods. The most faithful model (CPB) is seen in blue. The proposed transmon
basis model of highest precision ($\mathrm{DO_{3}}$) is seen in green, while
the generalised quantum Rabi (GR) and two-level quantum Rabi (R) models offer
significant accelerations in computation time as seen in red and black
respectively. This result communicates that the merit of the Rabi family of
models (GR and R) lies in their relatively low computational demand.

Computational efficiency is a clear advantage of the simplified models (GR and
R), attributable to both the smaller Hilbert space they require and their
diagonal matrix structure. Linear fits reveal that the $\mathrm{DO_{3}}$
modelling approach offers a marginal 6% runtime advantage over the full CPB
model, in comparison to GR and R offering accelerations of 3.6x and 16.3x
respectively. For the set of problems in which precise dynamics is expendable
in favour of simulation runtime, Fig. 11 is presented in order to demonstrate
the value of the simplified quantum Rabi models.

## References

  * Allen  _et al._ (2017) Joseph L. Allen, Robert Kosut, Jaewoo Joo, Peter Leek, and Eran Ginossar, “Optimal control of two qubits via a single cavity drive in circuit quantum electrodynamics,” Phys. Rev. A 95, 042325 (2017).
  * Spiteri  _et al._ (2018) Raymond J Spiteri, Marina Schmidt, Joydip Ghosh, Ehsan Zahedinejad, and Barry C Sanders, “Quantum control for high-fidelity multi-qubit gates,” New J. Phys. 20, 113009 (2018).
  * Werninghaus  _et al._ (2021) M. Werninghaus, D. J. Egger, F. Roy, S. Machnes, F. K. Wilhelm, and S. Filipp, “Leakage reduction in fast superconducting qubit gates via optimal control,” NPJ Quantum Inf. 7, 14 (2021).
  * Kirchhoff  _et al._ (2018) Susanna Kirchhoff, Torsten Keßler, Per J. Liebermann, Elie Assémat, Shai Machnes, Felix Motzoi, and Frank K. Wilhelm, “Optimized cross-resonance gate for coupled transmon systems,” Phys. Rev. A 97, 042348 (2018).
  * Wu  _et al._ (2020) Xian Wu, S. L. Tomarken, N. Anders Petersson, L. A. Martinez, Yaniv J. Rosen, and Jonathan L. DuBois, “High-fidelity software-defined quantum logic on a superconducting qudit,” Phys. Rev. Lett. 125, 170502 (2020).
  * García-Ripoll  _et al._ (2020) J. J. García-Ripoll, A. Ruiz-Chamorro, and E. Torrontegui, “Quantum control of frequency-tunable transmon superconducting qubits,” Phys. Rev. Appl. 14, 044035 (2020).
  * Riaz  _et al._ (2019) Bilal Riaz, Cong Shuang, and Shahid Qamar, “Optimal control methods for quantum gate preparation: a comparative study,” Quantum Inf. Process. 18, 100 (2019).
  * An and Zhou (2019) Zheng An and D. L. Zhou, “Deep reinforcement learning for quantum gate control,” EPL 126, 60002 (2019).
  * Zhang  _et al._ (2019) Xiao-Ming Zhang, Zezhu Wei, Raza Asad, Xu-Chen Yang, and Xin Wang, “When does reinforcement learning stand out in quantum control? a comparative study on state preparation,” NPJ Quantum Inf. 5, 85 (2019).
  * Xu  _et al._ (2019) Han Xu, Junning Li, Liqiang Liu, Yu Wang, Haidong Yuan, and Xin Wang, “Generalizable control for quantum parameter estimation through reinforcement learning,” NPJ Quantum Inf. 5, 82 (2019).
  * Wauters  _et al._ (2020) Matteo M. Wauters, Emanuele Panizon, Glen B. Mbeng, and Giuseppe E. Santoro, “Reinforcement-learning-assisted quantum optimization,” Phys. Rev. Research 2, 033446 (2020).
  * Zahedinejad  _et al._ (2016) Ehsan Zahedinejad, Joydip Ghosh, and Barry C. Sanders, “Designing high-fidelity single-shot three-qubit gates: A machine-learning approach,” Phys. Rev. Appl. 6, 054005 (2016).
  * Niu  _et al._ (2019) Murphy Yeuzhen Niu, Sergio Boixo, Vadim N. Smelyanskiy, and Hartmut Neven, “Universal quantum control through deep reinforcement learning,” NPJ Quantum Inf. 5, 33 (2019).
  * Bukov  _et al._ (2018) Marin Bukov, Alexandre G. R. Day, Dries Sels, Phillip Weinberg, Anatoli Polkovnikov, and Pankaj Mehta, “Reinforcement learning in different phases of quantum control,” Phys. Rev. X 8, 031086 (2018).
  * Wu  _et al._ (2019) Re-Bing Wu, Haijin Ding, Daoyi Dong, and Xiaoting Wang, “Learning robust and high-precision quantum controls,” Phys. Rev. A 99, 042327 (2019).
  * Abdelhafez  _et al._ (2020) Mohamed Abdelhafez, Brian Baker, András Gyenis, Pranav Mundada, Andrew A. Houck, David Schuster, and Jens Koch, “Universal gates for protected superconducting qubits using optimal control,” Phys. Rev. A 101, 022321 (2020).
  * Koch  _et al._ (2007) Jens Koch, Terri M. Yu, Jay Gambetta, A. A. Houck, D. I. Schuster, J. Majer, Alexandre Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, “Charge-insensitive qubit design derived from the cooper pair box,” Phys. Rev. A 76, 042319 (2007).
  * Motzoi  _et al._ (2009) F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm, “Simple pulses for elimination of leakage in weakly nonlinear qubits,” Phys. Rev. Lett. 103, 110501 (2009).
  * Reich  _et al._ (2012) Daniel M. Reich, Mamadou Ndong, and Christiane P. Koch, “Monotonically convergent optimization in quantum control using krotov’s method,” J. Chem. Phys. 136, 104103 (2012).
  * Machnes  _et al._ (2018) Shai Machnes, Elie Assémat, David Tannor, and Frank K. Wilhelm, “Tunable, flexible, and efficient optimization of control pulses for practical qubits,” Phys. Rev. Lett. 120, 150401 (2018).
  * Khaneja  _et al._ (2004) Navin Khaneja, Timo Reiss, Cindie Kehlet, Thomas Schulte-Herbrüggen, and Steffen J. Glaser, “Optimal control of coupled spin dynamics: design of nmr pulse sequences by gradient ascent algorithms,” J. Magn. Reson. 172, 296–305 (2004).
  * Liebermann and Wilhelm (2016) Per J. Liebermann and Frank K. Wilhelm, “Optimal qubit control using single-flux quantum pulses,” Phys. Rev. Appl. 6, 024022 (2016).
  * Daraeizadeh  _et al._ (2020) S. Daraeizadeh, S. P. Premaratne, N. Khammassi, X. Song, M. Perkowski, and A. Y. Matsuura, “Machine-learning-based three-qubit gate design for the toffoli gate and parity check in transmon systems,” Phys. Rev. A 102, 012601 (2020).
  * Basilewitsch  _et al._ (2019) Daniel Basilewitsch, Francesco Cosco, Nicolino Lo Gullo, Mikko Möttönen, Tapio Ala-Nissilä, Christiane P Koch, and Sabrina Maniscalco, “Reservoir engineering using quantum optimal control for qubit reset,” New J. Phys. 21, 093054 (2019).
  * Ballester  _et al._ (2012) D. Ballester, G. Romero, J. J. García-Ripoll, F. Deppe, and E. Solano, “Quantum simulation of the ultrastrong-coupling dynamics in circuit quantum electrodynamics,” Phys. Rev. X 2, 021007 (2012).
  * Bao  _et al._ (2018) Seraph Bao, Silken Kleer, Ruoyu Wang, and Armin Rahmani, “Optimal control of superconducting gmon qubits using pontryagin's minimum principle: Preparing a maximally entangled state with singular bang-bang protocols,” Phys. Rev. A 97, 062343 (2018).
  * Willsch  _et al._ (2017) D. Willsch, M. Nocon, F. Jin, H. De Raedt, and K. Michielsen, “Gate-error analysis in simulations of quantum computers with transmon qubits,” Physical Review A 96 (2017), 10.1103/physreva.96.062302.
  * Magesan  _et al._ (2011) Easwar Magesan, J. M. Gambetta, and Joseph Emerson, “Scalable and robust randomized benchmarking of quantum processes,” Phys. Rev. Lett. 106, 180504 (2011).
  * Nanduri  _et al._ (2013) Arun Nanduri, Ashley Donovan, Tak-San Ho, and Herschel Rabitz, “Exploring quantum control landscape structure,” Phys. Rev. A 88, 033425 (2013).
  * Chow  _et al._ (2013) Jerry M Chow, Jay M Gambetta, Andrew W Cross, Seth T Merkel, Chad Rigetti, and M Steffen, “Microwave-activated conditional-phase gate for superconducting qubits,” New Journal of Physics 15, 115012 (2013).
  * Magesan and Gambetta (2020) Easwar Magesan and Jay M. Gambetta, “Effective hamiltonian models of the cross-resonance gate,” Physical Review A 101 (2020).
  * Negîrneac  _et al._ (2021) V. Negîrneac, H. Ali, N. Muthusubramanian, F. Battistel, R. Sagastizabal, M. S. Moreira, J. F. Marques, W. J. Vlothuizen, M. Beekman, C. Zachariadis, N. Haider, A. Bruno, and L. DiCarlo, “High-fidelity controlled- z gate with maximal intermediate leakage operating at the speed limit in a superconducting quantum processor,” Physical Review Letters 126 (2021).
  * Chakrabarti and Rabitz (2007) Raj Chakrabarti and Herschel Rabitz, “Quantum control landscapes,” Int. Rev. Phys. Chem. 26, 671–735 (2007).
  * Moore  _et al._ (2008) Katharine Moore, Michael Hsieh, and Herschel Rabitz, “On the relationship between quantum control landscape structure and optimization complexity,” J. Chem. Phys. 128, 154117 (2008).
  * Larocca  _et al._ (2018) Martín Larocca, Pablo M. Poggi, and Diego A. Wisniacki, “Quantum control landscape for a two-level system near the quantum speed limit,” J. Phys. A Math. Theor. 51, 385305 (2018).
  * Day  _et al._ (2019) Alexandre G. R. Day, Marin Bukov, Phillip Weinberg, Pankaj Mehta, and Dries Sels, “Glassy phase of optimal quantum control,” Phys. Rev. Lett. 122, 020601 (2019).
  * Kosut  _et al._ (2019) Robert L Kosut, Christian Arenz, and Herschel Rabitz, “Quantum control landscape of bipartite systems,” J. Phys. A Math. Theor. 52, 165305 (2019).
  * Youssry  _et al._ (2020) Akram Youssry, Gerardo A. Paz-Silva, and Christopher Ferrie, “Characterization and control of open quantum systems beyond quantum noise spectroscopy,” NPJ Quantum Inf. 6, 95 (2020).
  * Asfaw  _et al._ (2020) A. Asfaw, L. Bello, Y. Ben-Haim, S. Bravyi, N. Bronn, L. Capelluto, A. Carrera Vazquez, J. Ceroni, R. Chen, A. Frisch, J. Gambetta, S. Garion, L. Gil, S. De La Puente Gonzalez, F. Harkins, T. Imamichi, D. McKay, A. Mezzacapo, Z. Minev, R. Movassagh, G. Nannicni, P. Nation, A. Phan, M. Pistoia, A. Rattew, J. Schaefer, J. Shabani, J. Smolin, K. Temme, M. Tod, S. Wood, and J. Wootton, “Learn quantum computation using qiskit,”  (2020).
  * Johansson  _et al._ (2012) J.R. Johansson, P.D. Nation, and Franco Nori, “Qutip: An open-source python framework for the dynamics of open quantum systems,” Comput. Phys. Commun. 183, 1760–1772 (2012).