Modeling low- and high-frequency noise in transmon qubits with resource-efficient measurement

TL;DR

This paper develops a resource-efficient model of noise in transmon qubits, capturing both low- and high-frequency effects, validated through experiments and fidelity measurements, improving understanding of decoherence mechanisms.

Contribution

The paper introduces a hybrid noise model combining low- and high-frequency effects using the Redfield master equation, validated with IBM quantum devices, and emphasizes the importance of including both noise types for accurate predictions.

Findings

Model accurately predicts fidelity decay under various conditions.

Including both high- and low-frequency noise improves model accuracy.

Experimental validation confirms the model's effectiveness.

Abstract

Transmon qubits experience open system effects that manifest as noise at a broad range of frequencies. We present a model of these effects using the Redfield master equation with a hybrid bath consisting of low and high-frequency components. We use two-level fluctuators to simulate 1/f-like noise behavior, which is a dominant source of decoherence for superconducting qubits. By measuring quantum state fidelity under free evolution with and without dynamical decoupling (DD), we can fit the low- and high-frequency noise parameters in our model. We train and test our model using experiments on quantum devices available through IBM quantum experience. Our model accurately predicts the fidelity decay of random initial states, including the effect of DD pulse sequences. We compare our model with two simpler models and confirm the importance of including both high-frequency and 1/f noise in order to accurately predict transmon behavior.

Tables

Table 1. Table 1: System-bath parameter values extracted using the fitting procedure of Section   IV.2 , and corresponding to the minima indicated by the green circles in Fig.   3 for Quito (top row) and Lima (bottom row).

Params$/\left(2\pi \right)$	Quito	Lima
$g_x$ [MHz]	$5.734$	$4.782$
$g_z$ [MHz]	$4.413$	$9.393$
${\omega }_x^c$ [GHz]	$1.948$	$2.340$
${\omega }_z^c$ [MHz]	$5.690$	$5.979$
$b$ [MHz]	$0.598$	$0.323$
${\gamma }_{\max }$ [GHz]	$0.051$	$0.083$

Equations

$$H_{\mathrm{trans}}=4E_C{\left(\hat{n}-n_g\right)}^2-E_J\cos \hat{\phi }.$$

$$H_{\mathrm{trans}}=$$

$$\displaystyle-\frac{E_{J}}{2}\sum_{-n_{\max}}^{n_{\max}}\left(|n\rangle\langle n+1|+|n+1\rangle\langle n|\right)\$$

$$H_{\mathrm{trans}}^{\mathrm{eigen}}=S{H}_{\mathrm{trans}}{S}^{†}=\sum _{k\ge 0}{\omega }_k|k⟩⟨k|,$$

$$H_{\mathrm{sys}}=H_{\mathrm{trans}}^{\mathrm{eigen}}+\epsilon (t)\cos \left({\omega }_{\mathrm{d}}t+{\varphi }_{\mathrm{d}}\right)\hat{n},$$

$$\hat{n}\approx \sum _{k\ge 0}⟨k|\hat{n}|k+1⟩|k⟩⟨k+1|+h.c.$$

$${\tilde{H}}_{\mathrm{sys}}$$

$$+\frac{\epsilon (t)}{2}\sum _{k\ge 0}{g}_{k,k+1}(|k⟩⟨k+1|+|k+1⟩⟨k|),$$

$${\int }_0^{t_g}\epsilon (t)dt=\pi ,$$

$$\epsilon (t)=\epsilon \left[G(t,t_g,\sigma )-G(0,t_g,\sigma )\right]\left(\Theta (t)-\Theta \left(t-t_g\right)\right),$$

$$G(t,t_g,\sigma )=\exp \left(-\frac{{\left(t-t_g/2\right)}^2}{2{\sigma }^2}\right).$$

$$H_X(t)=\frac{\epsilon (t)}{2}\left(|0⟩⟨1|+|1⟩⟨0|\right)=\frac{\epsilon (t)}{2}{\sigma }^x,$$

$$F_g=\frac{1}{6}\sum _{j=\pm x,\pm y,\pm z}\mathrm{Tr}\left[U_{\mathrm{ideal}}{\rho }_j^{1\mathrm{q}}{U}_{\mathrm{ideal}}^{†}\Pi [\rho (t_g)]\right],$$

$$\tilde{\epsilon }(t)=\epsilon (t)+i\alpha \frac{\dot{\epsilon }(t)}{{\eta }_q},$$

$$H_{\mathrm{SB}}=\sum _{i=x,y,z}{g}_i{A}_i\otimes B_i,$$

$$H_{\mathrm{SB}}=g_x{A}_x\otimes B_x+A_z\otimes \left(g_z{B}_z+\sum _k{b}_k{\chi }_k\left(t\right)I_{\mathrm{B}}\right),$$

$$A_x$$

$$A_z$$

$$C_{ij}(t)=\mathrm{Tr}\{U_{\mathrm{B}}(t)B_i{U}_{\mathrm{B}}^{†}(t)B_j{\rho }_{\mathrm{B}}\},$$

$$\rho_{\mathrm{B}}=e^{-\beta H_{\mathrm{B}}}/\mathrm{Tr}\big{(}e^{-\beta H_{\mathrm{B}}}\big{)}\ ,$$

$$\frac{\partial {\rho }_{\mathrm{S}}}{\partial t}=-i[H_{\mathrm{sys}}+A_z\sum _{k=1}^{10}{b}_k{\chi }_k(t),{\rho }_{\mathrm{S}}]+{\mathcal{L}}_{\mathrm{R}}({\rho }_{\mathrm{S}}),$$

$${\mathcal{L}}_{\mathrm{R}}({\rho }_{\mathrm{S}})=-\sum _{i=x,z}[A_i,{\Lambda }_i(t){\rho }_{\mathrm{S}}(t)]+h.c.,$$

$${\Lambda }_i(t)={\int }_0^t{C}_i(t-\tau )U_{\mathrm{sys}}(t,\tau )A_i{U}_{\mathrm{sys}}^{†}(t,\tau )\mathrm{d}\tau ,$$

$${\gamma }_j(\omega )={\int }_{-\infty }^{\infty }{C}_j(\tau )e^{i\omega \tau }\mathrm{d}\tau .$$

$${\gamma }_j(\omega )=2\pi \eta g_j^2\frac{\omega {\mathrm{e}}^{-|\omega |/{\omega }_j^c}}{1-{\mathrm{e}}^{-\beta \omega }},$$

$$C_{|\psi ⟩,s}=\frac{1}{N}\sqrt{\sum _{i=0}^{N-1}{\left(P_{|\psi ⟩,s}^{\mathrm{Sim}}(t_i)-P_{|\psi ⟩,s}^{\mathrm{Exp}}(t_i)\right)}^2},$$

$$a\equiv i\sum _{k\ge 0}{g}_{k,k+1}|k⟩⟨k+1|=i\sum _{k\ge 0}{\tilde{g}}_k\sqrt{k+1}|k⟩⟨k+1|,$$

$$g_{k,k+1}$$

$$\equiv {\tilde{g}}_k\sqrt{k+1},$$

$${\tilde{g}}_k\approx g\left(1-\frac{k}{2}\frac{{\eta }_q}{{\omega }_q}\right).$$

$$H_{\mathrm{sys}}$$

$$H_{\mathrm{drive}}$$

$${\tilde{H}}_{\mathrm{sys}}$$

Full text

I Introduction

Quantum computing based on superconducting qubits has made significant progress. Starting with the first implementations of superconducting qubits [1, 2], the field has developed several flavors of qubits, broadly classified as charge, flux, and phase qubits [3]. However, the real workhorse behind many of the recent critical developments [4, 5, 6, 7, 8, 9, 10, 11, 12, 13] in gate-based quantum computing is the transmon qubit [14]. Transmons are designed by adding a large shunting capacitor to charge qubits, the result being that they are almost insensitive to charge noise. Transmon-based cloud quantum computers (QCs) are now widely available to the broad research community for proof of principle quantum computing experiments [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25].

Quantum computers in their current form have high error rates. This includes coherent errors (originating from imperfect gates), state preparation and measurement errors (SPAM), and incoherent errors (environment-induced noise). The latter, which results in dephasing and relaxation errors, is a pernicious problem in quantum information processing. Characterizing and modeling these open quantum system effects is crucial for advancing the field and improving the prospects of fault-tolerant quantum computation [26, 27, 28]. Various procedures for modeling decoherence and control noise affecting idealized qubits have been discussed [29, 30, 31]. Still, modeling noise effects from first principles, i.e., starting at the circuit level of transmons and including $1/f$ noise, is relatively unexplored [32].

In this work, we develop a framework to model environment-induced noise effects on a transmon qubit using the master equation formalism. We use a hybrid quantum bath with an Ohmic-like noise spectrum to model dephasing and relaxation processes in multi-level transmons. We also include classical fluctuators and use a hybrid Redfield model to account for both low- ($1/f$) and high-frequency noise. We develop a simple noise learning procedure relying on dynamical decoupling (DD) [33, 34, 35, 36, 37] to obtain the noise parameters (see Ref. [38] for early experimental work in this area). Our procedure relies only on measurements of quantum state fidelity with and without a single type of DD sequence, and so is quite resource-efficient compared to protocols requiring full quantum state tomography or DD-based spectral analysis. We test our noise model via fidelity preservation experiments on IBMQE processors for random initial states and find that the model can correctly capture these experiments. The model is, moreover, capable of reproducing the effects of time-dependent dynamical decoupling pulses on the main qubit. Finally, we compare the predictions based on our model with two simpler models using ideal two-level qubits, excluding the fluctuators and assuming ideal, zero-width DD pulses. In contrast to our complete model, these simpler models fail to capture noise simultaneously in both the low and high-frequency regimes. As a result, whether with or without DD, they underperform in capturing fidelity preservation experiments.

This paper is organized as follows. In Section II, we develop our numerical method focusing on simulating multi-level transmon qubit and single-qubit gates, which form the DD sequences. Next, we discuss our open quantum system in Section III and describe our noise learning method in Section IV. We then test our learned model on Quito using DD experiments with random initial single-qubit states in Section V. We extend our method to Lima, which relies on a different calibration procedure, in Section VI. We conclude in Section VII. The appendix provides additional details and calculations in support of the main text.

II Numerical model of transmons

In this section, we focus on the circuit-level description of the transmon qubit that we use in our model. We start with the transmon Hamiltonian and find an effective Hamiltonian to simulate single-qubit time-dependent microwave gates. We include the Derivative Removal of Adiabatic Gates (DRAG) [39] technique in our numerical model. DRAG is a state-of-the-art technique used to improve the performance of single-qubit gates by suppressing leakage and phase errors. The former refers to non-zero population of non-computational levels at the end of a pulse, while the latter is a type of coherent error that results from the non-zero population in non-computational levels during the pulse: the admixture of such levels leads to a phase shift of the computational levels, resulting in a net phase error at the end of the pulse. By including the DRAG technique (used in the IBMQE devices) and considering the residual errors it is unable to suppress, we more accurately model the transmon behavior.

II.1 Transmon Hamiltonian

The Hamiltonian of a fixed-frequency transmon qubit is [14]:

[TABLE]

We work in units where $\hbar=1$. $E_{C}=e^{2}/(2C)$ is the charging energy ($C$ is the capacitance, and $e$ is the electron charge), $E_{J}=I_{C}/(2e)$ is the potential energy of the Josephson junction ($I_{C}$ is the critical current of the junction) and $n_{g}$ represents the charge offset number which can result in charge noise. In the operating regime of a transmon qubit, i.e., $E_{J}/E_{C}\gg 1$, the lowest few energy levels of the transmon are almost immune to charge noise, in which case $n_{g}$ can be safely ignored. The two operators $\hat{n}$ and $\hat{\varphi}$ are a canonically conjugate pair analogous to momentum and position. They satisfy the commutation relation $\left[\hat{n},\hat{\varphi}\right]=i$; $\hat{n}$ is the number operator for the Cooper pairs transferred between the superconducting islands of the Josephson junction, and $\hat{\varphi}$ is the gauge invariant phase difference across the Josephson junction, i.e., between the islands.

II.2 Time-dependent drives

To numerically simulate the time-dependent drive-pulses or gates, we start with Eq. 1 and write it in the charge basis (the eigenbasis of $\hat{n}$) such that the number of Cooper pairs takes values from $-n_{\max}$ to $n_{\max}$. Eq. (1) thus reduces to

[TABLE]

where we have taken $n_{g}=0$ since we are in the transmon regime. We truncate to $n_{\max}$ (later we set $n_{\max}=50$) and diagonalize the resulting Hamiltonian:

[TABLE]

where $\omega_{k}$ for $k=0,1,...$ represents the energy of the $k^{\rm{th}}$ level in the transmon eigenbasis, and $S$ is the unitary similarity transformation. The eigenfrequencies are $\omega_{ij}\equiv\omega_{i}-\omega_{j}$. The bare qubit frequency is $\omega_{q}\equiv\omega_{10}$ and the anharmonicity is $\eta_{q}\equiv\omega_{10}-\omega_{21}$. Since $\omega_{q}$ and $\eta_{q}$ are the two quantities accessible via experiments, we use these values to obtain the fitting parameters $E_{C}$ and $E_{J}$ in Eq. 1, which is the starting point of our transmon model.111In more detail, we try different values of $E_{J}$ and $E_{C}$ in Eq. 1 by diagonalizing the corresponding $H_{\rm{trans}}$ and comparing the result with the experimental values of $\omega_{q}$ and $\eta_{q}$. Once we find the values of $E_{J}$ and $E_{C}$ yielding the closest match, we proceed to Eq. 3 to find the full spectrum.

Next, we add the coupling to the microwave drive, which couples to the transmon charge operator. The total system Hamiltonian can then be written as

[TABLE]

where $\varepsilon(t)$ is the pulse envelope, $\omega_{\rm d}$ is the drive frequency, and $\phi_{\rm d}$ is the phase of the drive. We can simplify Eq. 4 by first writing the charge operator in the transmon eigenbasis of Eq. 3, i.e., $\hat{n}=\sum_{k,k^{\prime}}\bra{k}\hat{n}\ket{k^{\prime}}|{k}\rangle\!\langle k^{\prime}|$ and considering the charge coupling matrix elements. Only nearest-level couplings $\bra{k}\hat{n}\ket{k\pm 1}$ are found to be non-negligible, allowing us to ignore all higher order terms:

[TABLE]

Transforming into a frame rotating with the drive and employing the rotating wave approximation (RWA), we obtain, for $\phi_{\rm d}=0$, the effective Hamiltonian

[TABLE]

where $g_{k,k+1}\equiv\bra{k}\hat{n}\ket{k+1}$. By tuning $\phi_{\rm d}$, we can implement a rotation about any axis in the $(x,y)$ plane of the qubit subspace (after an additional projection). In particular, taking $\phi_{\rm d}=0$ or $\pi/2$ corresponds to a rotation about the $x$ or $y$ axis, respectively. Appendix A provides a derivation of Eq. 6 from Eq. 4.

The pulse envelope $\varepsilon(t)$ plays a vital role in the final implementation of the gate. Since we are interested mainly in applying $\pi$ pulses, we choose

[TABLE]

where $t_{g}$ is the pulse or gate duration. For our numerical simulations, we choose Gaussian-shaped pulses with envelopes given by

[TABLE]

where

[TABLE]

Here $\varepsilon$ is the maximum drive amplitude during the pulse, $\Theta(t)$ is the step function, and $\sigma$ is the standard deviation of the Gaussian pulse.

An essential aspect of gate design is that population should not leak to higher levels of the transmon, i.e., the drive pulses should be bandwidth-limited (adiabatic). An accurate measure of these off-resonant excitations is the Fourier transform of the pulse envelope at the detuning frequencies [40, 41]. For example, consider a Gaussian pulse with standard deviation $\sigma$. Its Fourier transform has a standard deviation proportional to $1/\sigma$, which means that the drive pulse applied at the qubit frequency $\omega_{q}$ has a frequency spread close to $1/\sigma$ about $\omega_{q}$. If $1/\sigma$ is of the order of the anharmonicity of the transmon, the pulse spectral width will overlap with some of the higher-level transitions. Fig. 1 shows the Gaussian pulse envelope and its Fourier transform, and illustrates how choosing a shorter gate time results in a larger frequency spread and vice versa. The use of DRAG pulses mitigates this leakage, as discussed further below.

In the two-state (qubit) subspace, Eq. 6 reduces to

[TABLE]

where $g_{0,1}=g$ has been absorbed into $\varepsilon(t)$ and $\sigma^{x}$ is the Pauli $X$ matrix. When $H_{X}(t)$ is evolved for a time $t_{g}$ such that Eq. 7 is satisfied, the resulting unitary is an ideal $X_{\pi}$ gate. To include the effect of higher levels, we first use the full gate Hamiltonian from Eq. 6 and then project the result to the qubit subspace.

The gate fidelity averaged over all input states in the qubit Hilbert space can be written as the average over the six polar states (i.e., the six eigenstates of $\sigma^{x},\sigma^{y}$, and $\sigma^{z}$) [42, 39]:

[TABLE]

where $\rho_{j}^{1\mathrm{q}}$ is the single qubit density matrix, and $U_{\rm ideal}$ represents the ideal unitary corresponding to the gate we wish to study. $\Pi[\rho(t_{g})]$ is the projection of the full density matrix onto the single qubit subspace. $F_{g}$ compares the density matrix $\rho$ after application of the gate (i.e., at $t=t_{g}$) with the expected density matrix in the qubit subspace.

To reduce phase errors caused by the presence of additional levels, a commonly used trick to implement single qubit gates such as $X_{\pi}$ is to break the gate into two halves where each half performs a $\pi/2$ rotation, accompanied by some virtual $Z$ rotations [43, 44]. Numerically, we observe that with four levels included in the transmon Hamiltonian and a total gate duration $t_{g}=70\;{\rm ns}$, with $\sigma=t_{g}/6$, the average single-qubit gate error $1-F_{g}$ can be suppressed by around $20\%$ if we use two such pulses instead of a single long pulse. The exact quantitative improvement depended on other model parameters. We observed this error reduction in a closed system setting with no environmental coupling, and so any fidelity improvement may be counteracted by open-system effects. For a detailed numerical study on the error of time-dependent gates with transmon qubits in the open system settings, see Ref. [45].

II.3 DRAG

Derivative Reduction by Adiabatic Gate (DRAG) [39, 46] is a useful technique to reduce both the leakage and the phase errors which accumulate during the operation of single-qubit gates. The most standard DRAG technique is to add the derivative of the pulse envelope to the quadrature component such that the final form of the pulse envelope $\tilde{\varepsilon}(t)$ is given by

[TABLE]

where $\eta_{q}$ is the transmon anharmonicity and $\alpha$ is a constant that can have different values depending on which errors need to be suppressed, e.g., $\alpha=1$ to suppress leakage and $\alpha=1/2$ to suppress coherent phase errors [47]. Eq. 12 implies that if we need to apply a single $X$ gate whose pulse envelope is given by $\varepsilon(t)$, then $\alpha{\dot{\varepsilon}(t)}/{\eta_{q}}$ needs to be applied along the $y$-axis.

IBMQE devices use DRAG, but the exact pulse parameters are not available to users. Instead, the value of $\alpha$ in our simulations can be optimized numerically to match the experimentally reported gate fidelity. This is the approach we take here, with the goal being to model experiments on IBMQE devices with single-qubit gate errors of the order of $10^{-3}$. This is the value reported using randomized benchmarking [48], which equals the average gate infidelity ($1-F_{g}$) when the gate set has gate-independent errors [49, 50], an assumption we make here to justify the use of $10^{-3}$ as our target gate infidelity.

We perform closed-system simulations and find that without DRAG, the gate infidelity (due to leakage and phase errors) is in the range $10^{-2}-10^{-3}$. With DRAG, we find that varying $\alpha$ from $1/2$ to $1$ increases the gate infidelity from $\sim 10^{-6}$ to $10^{-3}$, respectively. We thus choose $\alpha=1$ to match the reported fidelity. This suggests, as described in Ref. [47], that the remaining errors are mostly phase error. Indeed, we find that even without DRAG, single-qubit $X$ and $Y$ gates (both implemented using two $\pi/2$ pulses as explained above) have leakage errors well below $10^{-5}$. This is unsurprising, as these long gates are quite narrowband compared to the transmon anharmonicity. We note that this attributes all errors to coherent closed-system effects rather than decoherence. We expect incoherent errors to be of the order of $t_{g}/T_{2}\approx 5\times 10^{-4}$, suggesting that this choice is defensible (see Table 2). Furthermore, given that both gate fidelity and coherence times drift over hour-long timescales, we focus only on matching the correct order of magnitude for fidelity with our coherent error model.

III Open quantum system simulation

This section describes the noise model and discusses the hybrid Redfield equation used for the open quantum system simulations. For all simulations we truncate to the lowest four levels of the transmon qubit.

III.1 Interaction Hamiltonian

The single-qubit system bath interaction Hamiltonian in the lab frame can be written as

[TABLE]

where ${A}_{i}$ and ${B}_{i}$ represent the dimensionless system and bath coupling operators, respectively, and the coupling strengths $g_{i}$ have dimensions of energy. There are several contributions to decoherence and noise for a multi-level transmon circuit. With fixed-frequency architectures, charge noise and fluctuations in the critical current contribute most to decoherence. In contrast, in the flux-tunable variants of transmon qubits, the largest contribution comes from flux noise [14]. These considerations determine which coupling operators are needed to describe the noise model for a given architecture. In the IBMQE processors used, the transmons are fixed-frequency. We therefore choose appropriate noise operators below.

We consider the following system-bath interaction Hamiltonian:

[TABLE]

where the coupling operators ${A}_{x}$ and ${A}_{z}$ correspond to the charge coupling operator and to the Josephson energy operator and are defined as

[TABLE]

where $c_{1}$ and $c_{2}$ are fixed constants that depend on the charge energy $E_{c}$ and the Josephson energy $E_{J}$ of the transmon qubit, respectively. We expect, based on the discussion in Section II.2 – and observe in our simulations – that ${A}_{x}$ and ${A}_{z}$ act like $\sigma^{x}$ and $\sigma^{z}$ when projected into the qubit subspace. We find numerically that Eq. 14 is an adequate model accounting for the nearly equal decay of the $\ket{+}$ and $\ket{i}$ states, which is why we do not include a separate $\sigma^{y}$ coupling term. Note, however, that a (dependent) $\sigma^{y}$ component appears when we transform Eq. 14 from the lab frame into a frame rotating with the drive.

Previous studies have found that noise in the superconducting circuit can be separated into high and low-frequency components [51]. To account for this observation, we combine two noise models. We choose the bath operators $B_{x}$ and $B_{z}$ in Eq. 14 to be bosonic bath operators, which generally represent the high-frequency component of the noise. However, this is not always the case, as we argue in Section IV.

To account for the low-frequency noise component, which is a dominant noise source for superconducting qubits [52], we include a sum over classical fluctuators in Eq. 14, via the term proportional to the bath identity operator $I_{\rm B}$. This semiclassical contribution, when parameterized properly, can simulate the behavior of $1/f$ noise. We model the fluctuators as having equal coupling strengths, i.e., we set $b_{k}=b$ (with dimensions of energy) for $k=1,\cdots,10$. Each fluctuator can be characterized by a stochastic process $\chi_{k}(t)$ that switches between $\pm 1$ with a frequency $\gamma_{k}$, which is log-uniformly distributed between $\gamma_{\rm min}$ and $\gamma_{\max}$ [53].

III.2 Hybrid Redfield model

To simulate the reduced system dynamics of the interaction Hamiltonian in Eq. 14, we use a hybrid form of the Redfield (or TCL$2$) master equation [55, 56]. We first define the standard bath correlation function

[TABLE]

where $U_{\mathrm{B}}(t)=e^{-iH_{\rm B}t}$ is the unitary evolution operator generated by the bath Hamiltonian $H_{\mathrm{B}}$, and the reference state $\rho_{\mathrm{B}}$ is the Gibbs state of $H_{\mathrm{B}}$:

[TABLE]

where $\beta=1/T$ is the inverse temperature. Assuming the bath operators $B_{x}$ and $B_{z}$ are uncorrelated, i.e., $C_{xz}(t)=C_{zx}(t)=0$, we construct the following hybrid Redfield equation

[TABLE]

where $\mathcal{L}_{\mathrm{R}}$ is the Redfield Liouville superoperator

[TABLE]

and

[TABLE]

where $C_{j}(\tau)\equiv C_{jj}(\tau)$ [from LABEL:{eq:C_ij}] and $U_{\mathrm{sys}}(t)$ is the unitary evolution operator generated by the system Hamiltonian $H_{\mathrm{sys}}$. The reduced system dynamics are obtained by averaging the solution of Eq. 18 over all the realizations of $\chi_{k}(t)$ for $k=1,...,10$.

The bath component correlation functions $C_{j}(\tau)$ are the Fourier transforms of the bath component noise spectra

[TABLE]

We choose the bath to be Ohmic, which means that the component noise spectra have the form

[TABLE]

where $\omega^{c}_{j}$ is the cut-off frequency for bath operator $B_{j}$, and $\eta$ is a positive constant with dimensions of time squared arising in the specification of the Ohmic spectral function.

Lastly, the hybrid Redfield equation (18) can be transformed into a frame rotating with the drive frequency $\omega_{\rm d}$ (see Appendix A for details) by replacing every operator with the interaction-picture one [specifically, the $A_{i}$ operator in Eq. 20 needs to be replaced by $A_{i}(\tau)$]. We simulate the Redfield master equation in this rotating frame in the methodology and results we discuss next.

IV Methodology and Fitting Results

This section discusses our methodology for modeling a transmon qubit’s open quantum system behavior in a multi-qubit processor. We refer to the qubit of interest as the main qubit and all the others as spectator qubits. The goal is to extract the bath parameters in our open quantum system model and then use this model to predict the outcomes of experiments on the main qubit, including dynamical decoupling sequences. We treat qubit 1 (Q1) of the Quito processor as our main qubit. We are interested only in the main qubit’s behavior here; hence, we measure only the main qubit. Appendix B describes the procedure to extract and analyze the experimental data.

IV.1 Free and DD Evolution Experiments

Our procedure involves two types of experiments, as shown in Fig. 2. The first type, which we call a free-evolution experiment, consists of initializing all the qubits in a given state by applying a particular unitary operation $U3(\theta,\phi,\lambda)$ [57] (denoted as $U$ in Fig. 2) to each of the qubits. We then apply a sequence of identity gates on the main qubit, which we vary in number. Simultaneously we also apply the XY4 DD sequence to all the other (spectator) qubits (i.e., $Xf_{\tau}Yf_{\tau}Xf_{\tau}Yf_{\tau}$, where $f$ denotes free-evolution in the absence of pulses for a duration of $\tau$ [37]) for the same total duration as that of the identity gates on the main qubit. As shown in [54], DD sequences applied to spectator qubits suppress unwanted $ZZ$-interactions, i.e., $ZZ$-crosstalk between the main qubit and the spectator qubits. Without crosstalk suppression, we observe oscillations in the probability decay as a function of time [20, 29]; see Appendix C. With the crosstalk suppression scheme, i.e., DD applied to the spectator qubits, these oscillations disappear, and the main qubit is now primarily affected only by environment-induced noise. Finally, we apply the inverse of $U3(\theta,\phi,\lambda)$ and measure in the $Z$-basis. The result is how we compute the initial states’ decay probability.

Everything remains the same in the second type of experiment, which we call DD-evolution, except that we now apply the XY4 sequence to the main qubit and identity gates on the spectator qubits. As discussed in [54], when we apply the XY4 sequence only to the main qubit, we suppress the $ZZ$-crosstalk between the latter and all the spectator qubits and also decouple unwanted interactions between the main qubit and the environment. We note that, in contrast to experiments using DD to perform noise spectroscopy, here we use only a single type of DD sequence and do not vary any of its parameters.

IV.2 Fitting Procedure

We perform the free-evolution and DD-evolution experiments for the six Pauli states as initial states, i.e., we choose $U3(\theta,\phi,\lambda)$ to prepare $\ket{0}$, $\ket{1}$, $\ket{+}$, $\ket{-}$, $\ket{i}$ and $\ket{-i}$, and use the hybrid Redfield equation described in Section III.2 to simulate the dynamics of these experiments. We sweep over different values of the bath parameters and obtain the simulated probability decay as a function of time. To identify the simulation parameters that optimally match the experimental results, we define a cost function $C$ for a given initial state $\ket{\psi}$ as the $l_{2}$ norm distance between the experimental probabilities $P_{\ket{\psi},s}^{\rm Exp}(t_{i})$ and the simulation probabilities $P_{\ket{\psi},s}^{\rm Sim}(t_{i})$ for every instant $t_{i}$:

[TABLE]

where $s\in\{\rm{free},\rm{DD}\}$ and $N$ is the total number of instants. Note that we compensate for state preparation and measurement (SPAM) errors by shifting the experimental results such that in all cases, $P_{\ket{\psi},s}^{\rm Exp}(0)=1$.

We limit the number of free parameters requiring fitting to six: the coupling strengths $g_{x}$, $g_{z}$, and $b_{k}\equiv b$ [Eq. 14], and the cutoff frequencies $\gamma_{\max}$ (for $1/f$ noise), $\omega^{c}_{x}$, and $\omega^{c}_{z}$ [Eq. 22]. We set the bath temperature $T=20\;$mK ($\sim$ the fridge temperature), $\gamma_{\min}=10^{-4}\;$GHz, and $\eta=10^{-4}\;{\rm GHz}^{-2}$; these values are the same as in our previous work [54], which showed strong agreement between open system simulations and experiments using other IBMQE devices, and remain fixed throughout our fitting procedure. This procedure consists of three steps, which we detail next.

IV.2.1 Step I: free-evolution for $\ket{1}$

We first focus on the free-evolution experiment for the initial state $\ket{1}$, the first excited state in the transmon eigenbasis. Since the free-evolution for this state is only affected by charge noise, i.e., noise along the $x$-axis, the only contribution to the decay of $\ket{1}$ should come from the $g_{x}A_{x}\otimes B_{x}$ term in Eq. 14. Thus, we consider only this term in our numerical simulations for this initial state. For a given set of values of the coupling strength $g_{x}$ and bath cutoff frequency $\omega^{c}_{x}$, we compute the cost function $C_{\ket{1},\rm{free}}$ using Eq. 23, and obtain the contour plot shown in Fig. 3(a).

In our simulations, we vary $\omega^{c}_{x}/(2\pi)$ from $0.5$ to $3$ GHz and $g_{x}$ from [math] to $10^{-2}$ GHz, each with $20$ equidistant points so that the contour plot has a total of $400$ data points. We take the position of the global minimum of the cost function on this grid as the optimal set of bath parameter values. To reduce the resulting discretization uncertainty, we interpolate the contour plot and use the Nelder-Mead optimization method to locate the minima. We numerically find the global minimum at $\omega^{c}_{x}/(2\pi)=1.948$ GHz and $g_{x}/(2\pi)=0.573\times 10^{-2}$ GHz, denoted by the green circle in Fig. 3(a). With this, we have two out of the six bath parameters, and we use these learned parameters in the subsequent steps.

IV.2.2 Step II: DD-evolution for all six Pauli states

The second step involves the DD-evolution experiment for all six Pauli states. This requires including the term $g_{z}A_{z}\otimes B_{z}$ in Eq. 14, along with the first term whose bath parameters we already obtained. We do not include the semiclassical term in Eq. 14 consisting of fluctuators since it is expected to be strongly suppressed when DD is applied to the main qubit. We simulate time-dependent gates with DRAG corrections to model the DD pulses, as discussed in Section II. We are again left with just two bath parameters to optimize: $\omega^{c}_{z}$ and $g_{z}$. Fig. 3(b) shows the average of the cost function [Eq. 23] over the six Pauli states. The global minimum is found at $\omega^{c}_{z}/(2\pi)=0.569\times 10^{-2}$ GHz and $g_{z}/(2\pi)=0.441\times 10^{-2}$ GHz.

IV.2.3 Step III: free-evolution for $\ket{+}$

The final step requires optimizing the two remaining free parameters associated with the fluctuators: $\gamma_{\max}$ and $b$. Here we focus on the free-evolution experiment for initial state $\ket{+}$. We now employ the full system-bath Hamiltonian in Eq. 14 with the optimal parameters found in Steps I and II. Fig. 3(c) shows the contour plot for the cost function [Eq. 23], where, as in Step I, we again use $20$ different values of $\gamma_{\max}$ and $b$ each. The global minimum is found at $\gamma_{\max}/(2\pi)=0.051$ GHz and $b/(2\pi)=0.598\times 10^{-3}$ GHz.

IV.3 Methodology wrap-up

Let us briefly summarize our methodology and add a few technical details. As explained above, we extract the bath parameters by performing free-evolution experiments for two initial states ($\ket{1}$ and $\ket{+}$) and DD-evolution experiments for up to six initial states (the Pauli states). Our optimization procedure is iterative and is thus not guaranteed to yield the globally optimal values of all the bath parameters, but this is by design: we choose initial states that allow us to isolate the bath parameters one pair at a time, which renders the optimization problem tractable.

This methodology is quite general and can be used to characterize all the transmon qubits on a given transmon processor or, much more broadly, to characterize single qubits on any quantum information processing platform capable of supporting individual qubit gates and measurements, provided a sufficiently accurate and descriptive model of the qubits and the system-bath interaction is available. Our procedure inherently suppresses the effects of crosstalk due to the neighboring qubits via DD applied either to the spectator qubits (free-evolution experiments) or the main qubit (DD-evolution experiment), which reduces the number of free parameters of the noise model by eliminating the need to model crosstalk.

To obtain the contour plots shown in Fig. 3, we solve the Redfield master equation (Section III.2) for each point (i.e., each set of model parameters), requiring a total of $400$ simulation runs for each optimization. In the final step, including classical fluctuators to obtain Fig. 3(c), we use the trajectory version of the Redfield model introduced in Section III.2 to simulate a total of $600$ trajectories at each point. This is large enough to yield negligible error bars ($<2\times 10^{-2}$). The experimental results are obtained using the standard bootstrap method (see Appendix B). In defining the cost function [Eq. 23], we use the mean value of the experimental fidelity obtained after bootstrapping (see Appendix C) and ignore the associated tiny error bars ($\leq 6\times 10^{-3}$). These error bars are much smaller than the error induced by the discrete nature of our $40\times 40$ parameter grid, and so we can safely ignore them. We confirmed that varying the probabilities to the extremes of the error bars does not affect the values of the bath parameters we have extracted to the least significant digit we report. Table 1 summarizes the extracted values and the parameters we have fixed.

The accuracy of our results depends on the number of points in the contour plots in Fig. 3 (we used a $20\times 20$ grid for each panel). Even though we interpolate the otherwise discrete contour plots and find the minima over the resulting smooth surface, the limited number of points affects the precision of the learned bath parameters. Increasing this precision requires more sophisticated optimization techniques to speed up the process of obtaining the bath parameters. This becomes especially acute when extending the model to learning a multi-qubit system-bath Hamiltonian with correlated noise, as in this case, the number of bath parameters increases significantly. Here, we aim to demonstrate the model and methodology and illustrate both via the example of a single transmon qubit, and so we perform a simple brute-force search of the parameter space. Note that our methods for extracting the bath parameters also work with density matrices (from state tomography) instead of just probabilities. In that case, the $l_{2}$ norm distance in the cost function of Eq. 23 can be replaced by the trace-norm distance between the density matrices obtained from the simulation and the experiment. However, quantum state tomography imposes a much higher cost in terms of the number of required experiments and is thus less practical to scale up with a larger number of qubits. Our protocol requires only fidelity measurements and so is more resource-efficient.

V Model Prediction Results

V.1 Full Model

We now test our model for different initial states of the main qubit of the Quito processor. Since we always apply the DD sequence to the spectator qubits during the free-evolution experiments, the initial states of the latter do not matter due to the suppressed $ZZ$ coupling. This section considers a total of $16$ initial states, consisting of the six Pauli states and ten Haar-random states. We model the experimental results using the bath parameters we extracted in the previous section (Table 1). This serves as a stringent test of the model: we now use the previously fitted model to predict the outcome of experiments not included in Steps I-III of Section IV.2, i.e., the results with different initial states. The data for all the experiments (both fitting and testing) was obtained in one batch. We used only the data for the six Pauli initial states needed for Steps I-III to perform the fitting. We used the data for all $16$ initial states in the testing phase.

We consider the same two kinds of experiments: free-evolution and DD-evolution. Fig. 4 (top row) shows our model’s prediction accuracy for the ten random and six Pauli states. The top left panel [Fig. 4(a)] corresponds to the free-evolution case, where we present the relative error in the prediction of our model as a function of time compared with the experimental results. The relative error is defined as $({\rm mean_{exp}}-{\rm mean_{sim}})/{\rm mean_{exp}}$, where ${\rm mean_{exp}}$ is the bootstrapped average over $8192$ experimental repetition and ${\rm mean_{sim}}$ is the average over $600$ trajectory simulations of the hybrid Redfield model for any given time instant. The box plot contains the spread in the relative error over all $16$ states, showing that the relative error of the model is always below $8\%$ over the total time considered here. The median and the mean over the $16$ states are confined well below $3\%$ for every instant.

The performance of our model for the DD-evolution experiments is shown in the top right panel [Fig. 4(b)]. Here the relative error is always below $2\%$. The median and the mean are below $1\%$. The closer agreement of the model with the DD-evolution experiments is expected, given that in contrast to the free-evolution experiments, DD suppresses the low-frequency noise affecting the main qubit, and the limitations of our fluctuator model of this noise is a likely source of modeling error.

The first and fourth columns of Fig. 5 show the relative error of our model over the $16$ states and all $70$ instants of the free and DD-evolution experiments, respectively. The results of the latter are better, as expected from Fig. 4. In both cases, however, we observe that the model has a relative prediction error of just a few percent.

V.2 Simplified Models

As discussed in Section II, our numerical simulations use the circuit model Hamiltonian of a transmon qubit truncated to the four lowest transmon eigenstates. The gates are applied with time-dependent pulses of non-zero duration. To test the robustness of our detailed model and learning procedure, we compare it with two simpler models, SM1 and SM2, derived from our detailed model. The simpler models use a more straightforward qubit description where we truncate the transmon Hamiltonian to only two levels. The time-dependent gates are replaced with instantaneous (zero-duration), ideal gates. Moreover, we focus only on the Ohmic bath terms in Eq. 14, thus simplifying the noise model by removing the classical fluctuators. To test these simpler models’ predictive power, we follow the same procedure as in Section IV, but using only Steps I and II. The difference between SM1 and SM2 lies in Step II, where SM1 uses the DD-evolution experiments for the six Pauli states, whereas SM2 uses the free-evolution experiments for the same states. In both cases, we extract the model parameters and then use the resulting learned models to predict the outcomes of both the free-evolution and the DD-evolution experiment.

Fig. 5 shows the comparison between our detailed model and the simpler models SM1 and SM2. We observe that SM1 has the largest relative error for the free-evolution experiments, whereas SM2 has the largest relative error for the DD-evolution case. Our full model has the smallest relative error among the three models considered here for both the free and DD evolution experiments. However, the performance of SM1 and SM2 is essentially indistinguishable from the full model results in the DD and free-evolution cases, respectively. This is not unexpected, given that SM1 (SM2) is trained on the DD (free) evolution experiments and predicts these well. In other words, SM1 (SM2) captures the high (low)-frequency noise well, as expected since for SM1, the use of DD suppresses most of the low-frequency noise, while for SM2, the use of free-evolution means that the low-frequency noise remains a dominant source of decoherence. The added value of the detailed model and the use of Step III is that this provides enough information to capture both the low and high-frequency components of the noise, which yields a more complete model with better predictive power. We do note that taking the qubit approximation and treating DD pulses as instantaneous does not seem to appreciably worsen the performance of the simple models in their regime of accuracy, as SM1 (SM2) are roughly as accurate as the full model in the DD (free) evolution case. This suggests that an intermediate model, taking the qubit and instantaneous-pulse approximations but retaining the fluctuators, may be accurate and computationally efficient.

VI Calibration-independent learning

For multi-qubit superconducting processors, calibrating single-qubit drive frequencies is crucial for gate operations. In the presence of $ZZ$ coupling, the state of the spectator qubits modifies the eigenfrequency of the main qubit. This results in different choices of calibration frequencies depending on the spectator qubits’ state [54]. So far, we have focused on one particular device (Quito), which is calibrated while keeping the spectator qubits in the $\ket{+}$ state (see Appendix C and Ref. [54]). The all-$\ket{+}$ or all-$\ket{0}$ are usually the two preferred choices for the spectators’ state while calibrating a given qubit in a multi-qubit processor. When we perform state protection experiments on the main qubit initialized in the $\ket{+}$ state while keeping all the spectator qubits in $\ket{0}$, there exists a frequency mismatch which results in $ZZ$ crosstalk oscillations (see Appendix C); to remove these oscillations, we applied DD to the spectator qubits before starting our noise learning procedure. When device calibration is performed while keeping the spectators in the $\ket{0}$ state, a similar state protection experiment does not result in any oscillations, as evidenced by our Lima results (see Appendix C).

Therefore, we extend our noise learning method to Lima in this section. Following our procedure from Section IV, we again perform free-evolution (no DD is applied to any qubit) and DD-evolution (XY4 is applied just to the main qubit) experiments. The only difference from the Quito case is that, for the reasons explained above, the free-evolution experiment does not require the application of DD to the spectator qubits to suppress crosstalk oscillations. Fig. 3 (bottom row) shows the contour plots for each of the three steps involved in our learning methodology as described in Section IV.2. We find the global minima at the parameter values given in Table 1. Comparing the Quito and Lima parameters in Table 1, we observe that the coupling strength $g_{z}$ of the Ohmic bath along the $z$-axis is roughly double for Lima, whereas the strength of the fluctuators is roughly double for Quito. This indicates that Quito is more prone to low-frequency ($1/f$) noise.

Fig. 4 (bottom row) shows the Lima prediction results for the $16$ different initial states described above, using the learned noise parameters. The bottom left panel [Fig. 4(c)] shows the relative error in the prediction of our model of the free-evolution experiments as a function of time compared to the experimental results for all $16$ states. The bottom right panel [Fig. 4(d)] shows the same for the DD-evolution experiments. The relative error is always below $17\%$ and $4\%$ for free-evolution and DD-evolution, respectively. Similar to the Quito case, the relative error is significantly lower for the DD-evolution experiments. For longer evolution times ($\gtrsim 14~{}\mu$s), the agreement worsens for the free-evolution experiments. As for the Quito results, the closer agreement of the model with the DD-evolution experiments is likely because DD suppresses the low-frequency noise affecting the main qubit, which dominates the free-evolution case’s simulation error.

Fig. 6 shows the time-integrated version of the Lima results of Fig. 4(c,d), where we have combined all $16$ states and $70$ time instants into one box each for the free-evolution and DD-evolution experiments. Except for a few outliers, almost all the data points for the free-evolution case have relative errors below $9\%$. The outliers all correspond to evolution times longer than $11\mu$s, as indicated by the color bar. Similarly, for DD-evolution, all data points have relative errors below $3\%$, except for a few outliers. There are two main reasons for the larger errors at longer evolution times. First, the Redfield equation is based on the weak coupling approximation, and its accuracy degrades as we increase the simulation time (for rigorous error bounds, see Ref. [58]). Second, as time increases, the effect of distant fluctuators is increasingly felt, thus reducing the accuracy of our fluctuator model, which relies on a fixed number of fluctuators. Adding more fluctuators and exploring a distribution of fluctuator strengths, as opposed to our assumption of a fixed fluctuator strength $b$, is expected to improve the predictive power of our model.

Finally, as another test of our learning procedure, we computed $T_{1}$ from the learned models, using the values we report in Table 2. We did this by simulating the fidelity decay of the $\ket{1}$ state for $19.6\mu$s and fitting $\exp(-t/T_{1})$ to estimate $T_{1}$. We find $T_{1}=92.5\mu$s for Quito and $86.5\mu$s for Lima, compared to the reported $T_{1}=98.6\mu$s and $105\mu$s, respectively. Our result gives the correct order of magnitude and is particularly reasonable for Quito. The discrepancy may be in part due to the relatively short simulation time of $19.6\mu$s (larger times become prohibitively expensive). In addition, the discrepancy for Quito is smaller because we use three rounds of experiments within the same calibration to model the noise, which removes short-time fluctuations via bootstrap averaging, whereas for Lima, we used only one repetition. As $T_{1}$ often drifts significantly over hour-long time scales, we would not expect the Lima prediction to line up exactly with the reported $T_{1}$.

VII Summary and Conclusions

This work presents a detailed noise model for transmon qubits consisting of both low ($1/f$-like) and high-frequency noise components based on a hybrid Redfield master equation. We designed an iterative three-step procedure to extract the unknown system-bath and bath parameters from a few simple “free-evolution” and “DD-evolution” experiments, as illustrated in Fig. 2, using the six Pauli matrix eigenstates. In both cases, we used dynamical decoupling (DD) to suppress diagonal ($ZZ$) qubit crosstalk [54] so that the remaining dominant noise effect on the main qubit (the qubit of interest) is decoherence. Our model treats the transmon qubit as a four-level system based on the circuit model description of transmons (Section II) and treats the DD pulses as realistic time-dependent gates subject to quantum control (DRAG).

Once the unknown system-bath and bath parameters are extracted, we compare the model predictions with new experiments and a larger set of initial states and demonstrate that the model predicts the experimental results of free-evolution and DD-evolution with a relative error below $8\%$ and $2\%$ for Quito, and below $9\%$ and $3\%$ for Lima, respectively. This is based on a test with the six Pauli matrix eigenstates and ten random states for a total duration of up to $19.6\;\mu s$. The relative errors are higher for larger times (see Fig. 4), as expected because the Redfield model is based on the weak-coupling approximation, and its accuracy degrades as the simulation time is increased [58].

To test the robustness of our model, we performed a comparison with two simpler, two-level models with instantaneous pulses; while these models capture either the low or the high-frequency noise, the full model captures both types of noise. Furthermore, our method is applicable independently of the particular device-calibration procedure followed, as witnessed by the agreement we find for both Quito and Lima – devices with different drive-frequency calibrations.

The low relative error we find in the case of DD-evolution experiments ($<2\%$ and $<3\%$ for Quito and Lima, respectively) suggests that our full noise model helps model gate dynamics under the influence of decoherence. The model can also be used to study several qubits in parallel as long as there is no direct crosstalk between these qubits. Extending our noise model beyond weak coupling and using it to analyze and improve entangling gates are promising future directions.

We hope this work will benefit experimental groups working with superconducting qubits by helping them understand and learn experimental noise using a first principles approach, which only requires a set of straightforward experiments.

Acknowledgements.

We thank Mostafa Khezri, Ka-Wa Yip, and Humberto Munoz Bauza for several helpful discussions. This material is based upon work supported by the National Science Foundation, the Quantum Leap Big Idea under Grant No. OMA-1936388. We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Quantum team.

Appendix A Derivation of time-dependent drives

We start with Eq. 4 from the main text, and to obtain Eq. 5, we focus on the charge coupling matrix elements for transmon qubits. The selection rules of the transmon qubit due to its cosine potential dictate that $g_{k,k\pm 2}=0$ $\forall k$ [14]. The next order coupling terms, $g_{k,k\pm 3}$, are proportional to the ratio of the anharmonicity and the qubit frequency: $\eta_{q}/\omega_{q}$ [59]. Therefore $g_{k,k\pm 1}$ (the coupling between nearest levels) is dominant, and all other couplings can be ignored, giving Eq. 5.

We can write the approximated charge operator $\hat{p}\equiv\sum_{k\geq 0}\bra{k}\hat{n}\ket{k+1}|{k}\rangle\!\langle k+1|+{\rm h.c.}$ of Eq. 5 in terms of effective creation and annihilation operators for the eigenstates of the transmon, i.e., $\hat{p}=i(a^{\dagger}-a)$, where

[TABLE]

with

[TABLE]

and $g\equiv g_{0,1}=\bra{0}\hat{n}\ket{1}=\tilde{g}_{0}$. We note that $g_{k,k+1}=g_{k,k+1}^{*}$, which follows from the fact that $\hat{n}$ is the number operator for Cooper pairs. We note further that to first order in $\eta_{q}/\omega_{q}$ [59]

[TABLE]

We defined $\tilde{g}_{k}(=\tilde{g}_{k}^{*})$ in Eq. 25b to include all higher order perturbative corrections and do not use the approximation Eq. 26 in our numerical calculations. However, we note that the leading order correction to $\tilde{g}_{k}$ is of the same order as $g_{k,k\pm 3}$, i.e., proportional to $\eta_{q}/\omega_{q}$.

Eq. 4 becomes

[TABLE]

For simplicity, we consider the case when $\phi_{\rm d}=0$ in Eq. 27b, and transform the Hamiltonian $H_{\rm sys}^{\rm eigen}$ into a frame rotating with $U_{\rm d}={e}^{i\omega_{\rm d}\hat{N}t}$. Here $\hat{N}=\sum_{k\geq 0}k|{k}\rangle\!\langle k|$ is the transmon number operator. In this frame, rotating with the drive frequency, the effective Hamiltonian is given by

[TABLE]

where in going from Eq. 28a to Eq. 28b, we used the fact that $\hat{N}$ commutes with $H_{\mathrm{trans}}^{\mathrm{eigen}}$ [Eq. 3].

Note that $a$ and $a^{\dagger}$ should not be confused with the harmonic oscillator raising and lowering operators since

[TABLE]

For the harmonic oscillator case, $\tilde{g}_{k}=1$ for all $k$ and we obtain the usual commutation relation $\left[a,a^{\dagger}\right]=I$. Despite this, we have, as for the harmonic oscillator:

[TABLE]

so that:

[TABLE]

Let us write the last term in Eq. 28c as

[TABLE]

Using Eq. 31, we then have:

[TABLE]

Making the rotating wave approximation, we ignore the fast-oscillating terms with frequencies $\pm 2\omega_{\rm d}$, and Eq. 32 reduces to

[TABLE]

Combining this with Eqs. 24 and 28, the effective Hamiltonian for a drive that results in a rotation about the $x$-axis (in the qubit subspace) is given by

[TABLE]

which is Eq. 6 of the main text. By tuning $\phi_{\rm d}$, we can implement a rotation about any axis in the $(x,y)$ plane in the qubit subspace.

Appendix B Data collection and analysis methodology

We used the IBMQE processor ibmq$\_$quito (Quito) and ibmq$\_$lima (Lima), whose layout is shown schematically in Fig. 7. For both Quito and Lima, we use qubit 1 (Q1) as the main qubit. These are five-qubit processors consisting of superconducting transmon qubits. Various calibration details and hardware specifications relevant to the qubits and gates used in this work are provided in Table 2.

For each initial state, we selected $70$ equidistant time instants ($19.6~{}\mu$s$/70$), with each such instant corresponding to an integer number of cycles of the XY4 DD sequences. We generated a circuit according to the scheme given in Fig. 2 for each such initial state and each such instant. We sent all $70$ such circuits in one job (the maximum allowed number of circuits per job is $75$), and each job was repeated $8192$ times. We ensured that all the jobs were sent consecutively within the same calibration cycle to avoid charge-noise-dependent fluctuations and variations in critical features over different calibration cycles.

We measured only the main qubit in the $Z$ basis, each measurement yielding either [math] or $1$. We computed the empirical fidelity $F^{(e)}$ as the number of favorable outcomes ([math]) to the total number of experiments ($8192$ per initial state and per measurement time instant). This is a proxy for the $F_{\ket{\psi}}={\rm Tr}\left[U^{-1}\mathcal{E}\left(U|{0}\rangle\!\langle 0|U^{-1}\right)U\right]$, where $U$ represents $U3(\theta,\phi,\lambda)$ and $\mathcal{E}$ represents the quantum map of the main qubit corresponding to any of the three types of experiments described in the main text and in Appendix C below.

Error bars were then generated using the standard bootstrapping procedure, where we resample (with replacement) counts out of the experimental counts’ dictionary (i.e., the list of $0/1$ measurement outcomes per state and instant) and create several new dictionaries. The final fidelity and error bars are obtained by calculating the mean and standard deviations over the fidelities of these newly resampled dictionaries. Using ten such resampled dictionaries of the counts sufficed to give small error bars. We report the final fidelity with $2\sigma$ error bars, corresponding to $95\%$ confidence intervals.

Appendix C Experimental fidelity results

For the Quito device, Fig. 8 shows the results of three different types of experiments for the $16$ states consisting of $6$ Pauli states and $10$ Haar-random states. In the first type of experiment, we apply a series of identity gates separated by barriers on all the qubits (main qubits and the spectator qubits). All the qubits always start in the $\ket{0}$ state. The second and third types are the experiments discussed in the main text: free-evolution, where we apply an XY4 DD sequence to the spectator qubits and identity gates on the main qubit, and DD-evolution, where we apply the XY4 DD sequences to the main qubit and identity gates to the spectator qubits (see Fig. 2).

For Lima, we only show two types of experiments in Fig. 9. The first is the free-evolution experiment, where we prepare some initial state of the main qubit, apply a series of identity gates, and measure the computational basis. The second type is the usual DD experiment, where we apply a series of XY4 sequences to the main qubit and the identity operation to all the spectator qubits. In both types of experiments, the spectator qubits are always initialized in the ground state $\ket{0}$. In the Lima case, applying a DD sequence to the spectator qubits is unnecessary.