Detecting and tracking drift in quantum information processors

TL;DR

This paper introduces a spectral analysis method to detect and localize time-dependent errors in quantum information processors, improving stability and characterization accuracy.

Contribution

The paper presents a fast, simple, and statistically sound spectral analysis technique applicable to various quantum characterization protocols to identify and mitigate drift.

Findings

Successfully detected and localized drift in quantum processors.

Implemented drift control techniques to suppress instability.

Validated method on simulations and trapped-ion experiments.

Abstract

If quantum information processors are to fulfill their potential, the diverse errors that affect them must be understood and suppressed. But errors typically fluctuate over time, and the most widely used tools for characterizing them assume static error modes and rates. This mismatch can cause unheralded failures, misidentified error modes, and wasted experimental effort. Here, we demonstrate a spectral analysis technique for resolving time dependence in quantum processors. Our method is fast, simple, and statistically sound. It can be applied to time-series data from any quantum processor experiment. We use data from simulations and trapped-ion qubit experiments to show how our method can resolve time dependence when applied to popular characterization protocols, including randomized benchmarking, gate set tomography, and Ramsey spectroscopy. In the experiments, we detect instability and localize its source, implement drift control techniques to compensate for this instability, and then demonstrate that the instability has been suppressed.

Figures

Tables

Table 1. Table 1: Estimated gate error rates . Error rates for the three gates in the two GST experiments, obtained using time-resolved GST. In the first experiment the error rates of the 𝖦 𝗂 subscript 𝖦 𝗂 \mathsf{G}_{\mathsf{i}} gate vary over time (see Fig.  3 C of the main text). The values reported here are the average over the experiment duration. All other error rates are static over the course of each experiment. See the text for definitions of these three different error rates.

Error rate	${ϵ}_{\diamond }\times {10}^3$		${ϵ}_{\text{f}}\times {10}^3$		${ϵ}_{\mathrm{𝗌𝗉𝖾𝖼},\text{f}}\times {10}^3$
Exp. #	1	2	1	2	1	2
${𝖦}_{𝗂}$	3.45	0.40	0.02	0.07	0.02	0.07
${𝖦}_{𝗑}$	1.25	0.38	0.03	0.06	0.03	0.06
${𝖦}_{𝗒}$	0.82	0.53	0.15	0.04	0.15	0.04

Equations

$$p_l(t)$$

$$F_{\omega i}=\sqrt{\frac{2^{1-{\delta }_{\omega ,0}}}{N}}\cos \left(\frac{\omega \pi }{N}\left(i+\frac{1}{2}\right)\right),$$

$$F_{\omega i}^{(\text{textsc}dct)}=\sqrt{\frac{2^{1-{\delta }_{\omega ,0}}}{N}}\cos \left(\frac{\omega \pi }{N}\left(i+\frac{1}{2}\right)\right),$$

$$\mathbf{s}=\mathbf{p}-\bar{p}1,$$

$$\tilde{\mathbf{z}}=\{\begin{array}{ll}\frac{F(\mathbf{x}-\bar{x}1)}{\sqrt{\bar{x}(1-\bar{x})}}& \text{if}0<\bar{x}<1,\\ (0,1,1,\dots {)}^{\mathsf{T}}& \text{else}.\end{array}$$

$$\mathcal{H}=\{p_{c,i}\in [0,1]\}.$$

$${\tilde{Z}}_{\omega }=\frac{{\tilde{X}}_{\omega }}{\sqrt{\bar{X}(1-\bar{X})}},$$

$$Z_i\approx Y_i=\frac{X_i-\bar{p}}{\sqrt{\bar{p}(1-\bar{p})}},$$

$${\tilde{Y}}_{\omega }=\sum _i{F}_{\omega i}{Y}_i.$$

$${\mu }_{\omega }\equiv \mathbb{E}({\tilde{Y}}_{\omega })=\sum _i\mathbb{E}\left[\frac{F_{\omega i}(X_i-\bar{p})}{\sqrt{\bar{p}(1-\bar{p})}}\right]=\frac{{\tilde{p}}_{\omega }}{\sqrt{\bar{p}(1-\bar{p})}}.$$

$${\nu }_{\omega }\equiv \mathbb{V}({\tilde{Y}}_{\omega })=\frac{1}{\bar{p}(1-\bar{p})}\sum _i{F}_{\omega i}^2{p}_i(1-p_i).$$

$${\nu }_{\omega }=1+\frac{{\Delta }_{\omega }}{\bar{p}(1-\bar{p})},$$

$${\Delta }_{\omega }=\sum _{i=0}^{N-1}{F}_{\omega i}^2\left(s_i[1-2\bar{p}]-s_i^2\right).$$

$${\Delta }_{\omega }=-\frac{\Vert \mathbf{s}{\Vert }_2^2}{N},$$

$$\Vert \mathbf{v}{\Vert }_2={\left(\sum _i{v}_i^2\right)}^{\text{nicefrac}12},$$

$${\Delta }_{\omega }=\frac{1}{\sqrt{2N}}\left({\tilde{s}}_{2\omega }[1-2\bar{p}]-{\tilde{q}}_{2\omega }\right)-\frac{\Vert \mathbf{s}{\Vert }_2^2}{N},$$

$$\underset{n\to \infty }{\lim }\frac{1}{s_n^{2+\delta }}\sum _{i=1}^n\mathbb{E}\left[|A_i-{\mu }_i{|}^{2+\delta }\right]=0,$$

$$\mathbb{E}\left[{\left|Q_{\omega i}-\mathbb{E}[Q_{\omega i}]\right|}^3\right]\le \mathcal{O}(|F_{\omega i}{|}^3),$$

$$\frac{1}{s_n^3}\sum _{i=0}^{N-1}\mathbb{E}\left[|Q_{\omega i}-\mathbb{E}[Q_{\omega i}]{|}^3\right]$$

$$=\mathcal{O}\left(\frac{1}{\sqrt{N}}\right),$$

$${\tilde{z}}_{\text{avg},\omega }^2=\frac{1}{C}\sum _{c=1}^C{\tilde{z}}_{c,\omega }^2,$$

$${\mathcal{H}}_{0,\omega }={\cap }_c{\mathcal{H}}_{0,c,\omega },$$

$$T={\mathsf{CDF}}_1^{-1}(1-\alpha ),$$

$$T_{\mathsf{individual}}={\mathsf{CDF}}_1^{-1}\left[1-\frac{(1-w)\alpha }{(N-1)C}\right],$$

$$T_{\mathsf{average}}=\frac{1}{C}{\mathsf{CDF}}_C^{-1}\left[1-\frac{w\alpha }{N-1}\right].$$

$$p(\gamma ,t)={\gamma }_0+\sum _{\omega \in \mathbb{W}}{\gamma }_{\omega }{f}_{\omega }(t),$$

$$f_{\omega }(t)=\cos \left(\frac{\omega \pi }{N}\left[\frac{t-t_0}{t_{\text{step}}}+\frac{1}{2}\right]\right),$$

$$t_{\text{step}}=\frac{t_{N-1}-t_0}{N-1},$$

$${\gamma }_{\omega }\to {\gamma }_{\omega }-\delta \text{sign}({\gamma }_{\omega }),$$

$${\gamma }_j(t)={\gamma }_{j,0}+\sum _{\omega \in {\mathbb{S}}_j}{\gamma }_{j,\omega }{f}_{\omega }(t),$$

$$k=\sum _j(|{\mathbb{S}}_j|+1),$$

Full text

Introduction

Recent years have seen rapid advances in quantum information processors (QIPs). Testbed processors containing 10s of qubits are becoming commonplace Rol et al. (2017); Otterbach et al. (2017); Friis et al. (2018); Arute et al. (2019) and error rates are being steadily suppressed Rol et al. (2017); Blume-Kohout et al. (2017), fueling optimism that useful quantum computations will soon be performed. Improved theories and models of the types and causes of errors in QIPs have played a crucial role in these advances. These new insights have been made possible by a range of powerful device characterization protocols Blume-Kohout et al. (2017); Merkel et al. (2013); Knill et al. (2008); Magesan et al. (2011); Proctor et al. (2019); Magesan et al. (2012); Cross et al. (2016); Barends et al. (2014); Carignan-Dugas et al. (2015); Gambetta et al. (2012); Kimmel et al. (2015) that allow scientists to probe and study QIP behavior. But almost all of these techniques assume that the QIP is stable — that data taken over a second or an hour reflect some constant property of the processor. These methods can malfunction badly if the actual error mechanisms are time-dependent Dehollain et al. (2016); Epstein et al. (2014); van Enk and Blume-Kohout (2013); Fong and Merkel (2017); Chow et al. (2009); Fogarty et al. (2015); Wan et al. (2019).

Yet temporal instability in QIPs is ubiquitous Wan et al. (2019); Harris et al. (2008); Bylander et al. (2011); Chan et al. (2018); Fogarty et al. (2015); Klimov et al. (2018); Megrant et al. (2012); Müller et al. (2015); Meißner et al. (2018); Klimov et al. (2018); De Graaf et al. (2018); Merkel et al. (2019); Burnett et al. (2019). The control fields used to drive logic gates drift Wan et al. (2019), $T_{1}$ times can change abruptly Burnett et al. (2019), low-frequency $1/f^{\alpha}$ noise is common Bylander et al. (2011), and laboratory equipment produces strongly oscillating noise (e.g., $50\text{\,}\mathrm{Hz}$/$60\text{\,}\mathrm{Hz}$ line noise and \sim$$1\text{\,}\mathrm{Hz} mechanical vibrations from refrigerator pumps). These intrinsically time-dependent error mechanisms are becoming more and more important as technological improvements suppress stable and better-understood errors. As a result, techniques to characterize QIPs with time-dependent behavior are becoming increasingly necessary.

In this article we introduce and demonstrate a general, flexible and powerful methodology for detecting and measuring time-dependent errors in QIPs. The core of our techniques can be applied to time-series data from any set of repeated quantum circuits — so they can be applied to most QIP experiments with only superficial adaptations — and they are sensitive to both periodic instabilities (e.g., $50\text{\,}\mathrm{Hz}$/$60\text{\,}\mathrm{Hz}$ line noise) and aperiodic instabilities (e.g., $1/f^{\alpha}$ noise). This means that they can be used for routine, consistent stability analyses across QIP platforms, and that they can be applied to data gathered primarily for other purposes, e.g., data from running an algorithm or error correction. Moreover, we show how to use our methods to upgrade standard characterization protocols — including randomized benchmarking (RB) Knill et al. (2008); Magesan et al. (2011); Proctor et al. (2019); Magesan et al. (2012); Cross et al. (2016); Barends et al. (2014); Carignan-Dugas et al. (2015); Gambetta et al. (2012) and gate set tomography (GST) Blume-Kohout et al. (2017); Merkel et al. (2013) — into time-resolved techniques. Our methods therefore induce a suite of general-purpose drift characterization techniques, complementing tools that focus on specific types of drift Cortez et al. (2017); Bonato and Berry (2017); Wheatley et al. (2010); Harris et al. (2008); Klimov et al. (2018); Chan et al. (2018); Bylander et al. (2011); Young and Whaley (2012); Gupta and Biercuk (2018); Granade et al. (2016, 2017); Huo and Li (2017); Kelly et al. (2016); Huo and Li (2018); Rudinger et al. (2019). We demonstrate our techniques using both simulations and experiments. In our experiments, we implemented high-precision, time-resolved Ramsey spectroscopy and GST on a 171Yb*+* ion qubit. We detected a small instability in the gates, isolated its source, and modified the experiment to compensate for the discovered instability. By then repeating the GST experiment on the stabilized qubit, we were able to show both improved error rates and that the drift had been suppressed.

Results

Instability in quantum circuits. Experiments on QIPs almost always involve choosing some quantum circuits and running them many times. The resulting data is usually recorded as counts Blume-Kohout et al. (2017); Merkel et al. (2013); Kimmel et al. (2015); Knill et al. (2008); Magesan et al. (2011); Proctor et al. (2019); Magesan et al. (2012); Cross et al. (2016); Barends et al. (2014); Carignan-Dugas et al. (2015); Gambetta et al. (2012) for each circuit — i.e., the total number of times each outcome was observed for each circuit. Dividing these counts by the total number of trials yields frequencies that serve as good estimates of the corresponding probabilities averaged over the duration of the experiment. But if the QIP’s properties vary over that duration, then the counts do not capture all the information available in the data, and time-averaged probabilities do not faithfully describe the QIP’s behavior. The counts may then be irreconcilable with any model for the QIP that assumes that all operations (state preparations, gates, and measurements) are time independent. This discrepancy results in failed or unreliable tomography and benchmarking experiments Dehollain et al. (2016); Epstein et al. (2014); van Enk and Blume-Kohout (2013); Fong and Merkel (2017); Fogarty et al. (2015); Chow et al. (2009); Wan et al. (2019).

Time-resolved analysis of the data from any set of circuits can be enabled by simply recording the observed outcomes (“clicks”) for each circuit in sequence, rather than aggregating this sequence into counts. We call the sequence of outcomes $\mathbf{x}=(x_{1},x_{2},\ldots,x_{N})$ obtained at $N$ data collection times $t_{1},t_{2},\ldots,t_{N}$ a “clickstream.” There is one clickstream for each circuit. We focus on circuits with binary 0/1 outcomes (see Appendix A for discussion of the general case), and on data obtained by “rastering” through the circuits. Rastering means running each circuit once in sequence, then repeating that process until we have accumulated $N$ clicks per circuit (see Fig. 1B). Under these conditions, the clickstream associated with each circuit is a string of bits, at successive times, each of which is sampled from a probability distribution over $\{0,1\}$ that may vary with time. If this probability distribution does vary over time, then we say that the circuit is temporally unstable. In this article we present methods for detecting and quantifying temporal instability, using clickstream data from any circuits, which are summarized in the flowchart of Fig. 1A.

Our methodology is based on transforming the data to the frequency domain and then thresholding the resultant power spectra. From this foundation we generate a hierarchy of outputs: 1) yes/no instability detection; 2) a set of drift frequencies; 3) estimates of the circuit probability trajectories; and 4) estimates of time-resolved parameters in a device model. To motivate this strategy, we first highlight some unusual aspects of this data analysis problem.

Formally, a clickstream $\mathbf{x}$ is a single draw from a vector of independent Bernoulli (“coin”) random variables $\mathbf{X}=(X_{1},X_{2},\dots,X_{N})$ with biases $\mathbf{p}=(p_{1},p_{2},\dots,p_{N})$. Here $p_{i}=p(t_{i})$ is the instantaneous probability to obtain 1 at the $i^{\rm th}$ repetition time of the circuit, and $p(\cdot)$ is the continuous-time probability trajectory. The naive strategy for quantifying instability is to estimate $\mathbf{p}$ from $\mathbf{x}$ assuming nothing about its form. However, $\mathbf{p}$ consists of $N$ independent probabilities and there are only $N$ bits from which to estimate them, so this strategy is flawed. The best fit is always $\mathbf{p}=\mathbf{x}$, which is a probability jumping between 0 and 1, even if the data seems typical of draws from a fixed coin. This is overfitting.

To avoid overfitting, we must assume that $\mathbf{p}$ is within some relatively small subset of all possible probability traces. Common causes of time variation in QIPs are not restricted to any particular portion of the frequency spectrum, but they are typically sparse in the frequency domain, i.e., their power is concentrated into a small range of frequencies. For example, step changes and $1/f^{\alpha}$ noise have power concentrated at low frequencies, while $50\text{\,}\mathrm{Hz}$/$60\text{\,}\mathrm{Hz}$ line noise has an isolated peak, perhaps accompanied by harmonics. Broad-spectrum noise does appear in QIP systems, but because it has an approximately flat spectrum, it acts like white noise — which produces uncorrelated stochastic errors that are accurately described by time-independent models. So, we model variations as sparse in the frequency domain, but otherwise arbitrary. Note that we do not make any other assumptions about $p(t)$. We do not assume that it is sampled from a stationary stochastic process, or that the underlying physical process is, e.g., strongly periodic, deterministic, or stochastic.

Detecting instability. The expected value of a clickstream is the probability trajectory, and this also holds in the frequency domain. That is, $\mathbb{E}[\tilde{\mathbf{X}}]=\tilde{\mathbf{p}}$ where $\mathbb{E}[\cdot]$ is the expectation value and $\tilde{\mathbf{v}}$ denotes the Fourier transform of the vector $\mathbf{v}$ (see the Methods for the particular transform that we use). In the time domain, each $x_{i}$ is a very low-precision estimate of $p_{i}$. In the frequency domain, each $\tilde{x}_{\omega}$ is the weighted sum of $N$ bits, so the strong, independent shot noise inherent in each bit is largely averaged out and any non-zero $\tilde{p}_{\omega}$ is highlighted. Of course, simply converting to the frequency domain cannot reduce the total amount of shot noise in the data. To actually suppress noise we need a principled method for deciding when a data mode $\tilde{x}_{\omega}$ is small enough to be consistent with $\tilde{p}_{\omega}=0$. One option is to use a regularized estimator inspired by compressed sensing Donoho (2006). But we take a different route, as this problem naturally fits within the flexible and transparent framework of statistical hypothesis testing Lehmann and Romano (2006); Shaffer (1995).

We start from the null hypothesis that all the probabilities are constant, i.e., $\tilde{p}_{\omega}=0$ for every $\omega>0$ and every circuit. Then, for each $\omega$ and each circuit we conclude that $|\tilde{p}_{\omega}|>0$ only if $|\tilde{x}_{\omega}|$ is so large that it is inconsistent with the null hypothesis at a pre-specified significance level $\alpha$. If we standardize $\mathbf{x}$, by subtracting its mean and dividing by its variance, then this procedure becomes particularly transparent: if the probability trace is constant, then the marginal distribution of each Fourier component $\tilde{X}_{\omega}$ for $\omega>0$ is approximately normal, and so its power $|\tilde{X}_{\omega}|^{2}$ is $\chi_{1}^{2}$ distributed. So if $|\tilde{x}_{\omega}|^{2}$ is larger than the $(1-\alpha)$-percentile of a $\chi^{2}_{1}$ distribution, then it is inconsistent with $\tilde{p}_{\omega}=0$. To test at every frequency in every circuit requires many hypothesis tests. Using standard techniques Lehmann and Romano (2006); Shaffer (1995), we set an $\alpha$-significance power threshold such that the probability of falsely concluding that $|\tilde{p}_{\omega}|>0$ at any frequency and for any circuit is at most $\alpha$ (i.e., we seek strong control of the family-wise error rate; see Appendix A).

We now demonstrate this drift detection method with data from a Ramsey experiment on a 171Yb*+* ion qubit suspended above a linear surface-electrode trap Blume-Kohout et al. (2017) and controlled using resonant microwaves. Shown in Fig. 1B, these circuits consist of preparing the qubit on the $\hat{x}$ axis of the Bloch sphere, waiting for a time $lt_{w}$ ($l=1,2,4,\dots,8192$, $t_{w}\approx$400\text{,}\mathrm{\SIUnitSymbolMicro s}$$), and measuring along the $\hat{y}$ axis. We performed $6000$ rasters through these circuits, over approximately 8 hours. A representative subset of the power spectra for this data are shown in Fig. 1C, as well as the $\alpha$-significance threshold for $\alpha=5\%$. The spectra for circuits containing long wait times exhibit power above the detection threshold, so instability was detected. These data are inconsistent with constant probabilities. Ramsey circuits are predominantly sensitive to phase accumulation, caused by detuning between the qubit and the control field frequencies, so it is reasonable to assume that it is this detuning that is drifting. The detected frequencies range from the lowest Fourier basis frequency for this experiment duration — approximately $15\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$ — up to approximately $250\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$. The largest power is more than $1700$ standard deviations above the expected value under the null hypothesis, which is overwhelming evidence of temporal instability.

Quantifying instability. Statistically significant evidence in data for time-varying probabilities does not directly imply anything about the scale of the detected instability. For instance, even the weakest periodic drift will be detected with enough data. We can quantify instability in any circuit by the size of the variations in its outcome probabilities. We can measure this size by estimating the probability trajectory $\mathbf{p}$ for each circuit (step 3, Fig. 1A). As noted above, the unregularized best-fit estimate of $\mathbf{p}$ is the observed bit-string $\mathbf{x}$, which is overfitting. To regularize this estimate, we use model selection. Specifically, we select the time-resolved parameterized model $p(t)=\gamma_{0}+\sum_{k}\gamma_{k}f_{k}(t)$, where $f_{k}(t)$ is the $k^{\rm th}$ basis function of the Fourier transform, the summation is over those frequencies with power above the threshold in the power spectrum, and the $\gamma_{k}$ are parameters constrained only so that each $p(t)$ is a valid probability. We can then fit this model to the clickstream for the corresponding circuit, using any standard data fitting routine, e.g., maximum likelihood estimation.

Estimates of the time-resolved probabilities for the Ramsey experiment are shown in Fig. 1D (unbroken lines). Probability traces are sufficient for heuristic reasoning about the type and size of the errors, and this is often adequate for practical debugging purposes. For example, these probability trajectories strongly suggest that the qubit detuning is slowly drifting. To draw more rigorous conclusions, we can implement time-resolved parameter estimation.

Time-resolved benchmarking & tomography. The techniques presented so far provide a foundation for time-resolved parameter estimation, e.g., time-resolved estimation of gate error rates, rotation angles, or process matrices. We introduce two complementary approaches, which we refer to as “non-intrusive” and “intrusive”, that can add time resolution to any benchmarking or tomography protocol. The non-intrusive approach is to replace counts data with instantaneous probability estimates in existing benchmarking/tomography analyses (step 4, Fig 1A). It is non-intrusive because it doesn’t require modifications to existing analysis codes. In constrast, the intrusive approach builds an explicitly time-resolved model and fits its parameters to the time-series data. We now detail and demonstrate these two techniques.

All standard characterization protocols, including all forms of tomography Blume-Kohout et al. (2017); Merkel et al. (2013) and RB Knill et al. (2008); Magesan et al. (2011); Proctor et al. (2019); Magesan et al. (2012); Cross et al. (2016); Barends et al. (2014); Carignan-Dugas et al. (2015); Gambetta et al. (2012), are founded on some time-independent parameterized model that describes the outcome probabilities for the circuits in the experiment, or a coarse-graining of them (e.g., mean survival probabilities in RB). When analyzing data from these experiments, the counts data from these circuits are fed into an analysis tool that estimates the model parameters, which we denote $\{\gamma_{i}\}$. To upgrade such a protocol using the non-intrusive method, we: (i) use the spectral analysis tools above to construct time-resolved estimates of the probabilities; (ii) for a given time, $t_{j}$, input the estimated probabilities directly into the analysis tool in place of frequencies; (iii) recover an estimate of the model parameters, $\{\gamma_{i}(t_{i})\}$ at that time; and (iv) repeat for all times of interest $\{t_{j}\}$. This non-intrusive approach is simple, but statistically ad hoc.

The intrusive approach permits statistical rigor at the cost of a more complex analysis. It consists of (i) selecting an appropriate time-resolved model for the protocol and (ii) fitting that model to the time-series data (steps 5a-5b, Fig. 1A). In the model selection step, we expand each model parameter $\gamma$ into a sum of Fourier components: $\gamma\to\gamma_{0}+\sum_{\omega}\gamma_{\omega}f_{\omega}(t)$, where the $\gamma_{\omega}$ are real-valued amplitudes, and the summation is over some set of non-zero frequencies. This set of frequencies can vary from one parameter to another, and may be empty if the parameter in question appears to be constant. To choose these expansions we need to understand how any drift frequencies in the model parameters would manifest in the circuit probability trajectories, and thus in the data.

To demonstrate the intrusive approach, we return to the Ramsey experiment. In the absence of drift the probability of “1” in a Ramsey circuit with a wait time of $lt_{w}$ is $p_{l}=A+B\exp(-l/l_{0})\sin(2\pi lt_{w}\Omega)$, where $\Omega$ is the detuning between the qubit and the control field, $\nicefrac{{1}}{{l_{0}}}$ is the rate of decoherence per idle, and $A,B\approx\nicefrac{{1}}{{2}}$ account for any state preparation and measurement errors. In our Ramsey experiment, the probability trace estimates shown in Fig. 1C suggest that the state preparation, measurement and decoherence error rates are approximately time-independent, as the contrast is constant over time. So we define a time-resolved model that expands only $\Omega$ into a time-dependent summation:

[TABLE]

where $\Omega(t)=\gamma_{0}+\sum_{\omega}\gamma_{\omega}f_{\omega}(t)$. To select the set of frequencies in the summation, we observe that the dependence of the circuit probabilities on $\Omega$ is approximately linear for small $l$ (e.g., expand Eq. (1) around $lt_{w}\Omega(t)\approx 0$). Therefore, the oscillation frequencies in the model parameters necessarily appear in the circuit probabilities. So in our expansion of $\Omega$, we include all 13 frequencies detected in the circuit probabilities (i.e., the ones with power above the threshold in Fig. 1C). The circuit probabilities will also contain sums, differences and harmonics of the frequencies in the true $\Omega$ — Fig. 1D shows clearly that the phase is wrapping around the Bloch sphere in the circuits with the longest wait times ($l\geq 2048$), so these harmonic contributions will be significant in our data. Therefore, this frequency selection strategy could result in erroneously including some of these harmonics in our model. We check for this using standard information-theoretic criteria Akaike (1974), and then discard any frequencies that should not be in the model (see Appendix B). This avoids over-fitting the data. Once the model is selected, we have a time-resolved parameterized model that we can directly fit to the time-series data. We do this with maximum likelihood estimation.

Fig. 1E shows the estimated qubit detuning $\Omega(t)$ over time. It varies slowly between approximately $-$0.5\text{,}\mathrm{Hz}$$ and $0.5\text{\,}\mathrm{Hz}$. The detuning is correlated with an ancillary measurement of the ambient laboratory temperature (the Spearman correlation coefficient magnitude is 0.92), which fluctuates by $\sim$1.5 C∘ over the course of the experiment. This suggests that temperature fluctuations are causing the drift in the qubit detuning (this conclusion is supported by further experiments: see later and the Methods). The detuning has been estimated to high precision, as highlighted by the $2\sigma$ confidence regions in Fig. 1E. As with all standard confidence regions, these are in-model uncertainties, i.e., they do not account for any inadequacies in the model selection. However, we can confirm that the estimated detuning is reasonably consistent with the data by comparing the $p_{l}(t)$ predicted by the estimated model (dotted lines, Fig. 1D) with the model-independent probability estimates obtained earlier (unbroken lines, Fig 1D). These probabilities are in close agreement.

Demonstration on simulated data. RB Knill et al. (2008); Magesan et al. (2011); Proctor et al. (2019); Magesan et al. (2012); Cross et al. (2016); Barends et al. (2014); Carignan-Dugas et al. (2015); Gambetta et al. (2012) and GST Blume-Kohout et al. (2017); Merkel et al. (2013) are two of the most popular methods for characterizing a QIP. Both methods are robust to state preparation and measurement errors; RB is fast and simple, whereas GST provides detailed diagnostic information about the types of errors afflicting the QIP. We now demonstrate time-resolved RB and GST on simulated data, using the general methodology introduced above. The number of circuits and circuit repetitions in these simulated experiments are in line with standard practice for these techniques, so they demonstrate that our techniques can be applied to RB and GST without additional experimental effort.

We simulated data from 2000 rasters through 100 randomly sampled RB circuits Knill et al. (2008); Magesan et al. (2011); Proctor et al. (2019) on two qubits. The error model consisted of 1% depolarization on each qubit and a time-dependent coherent $\hat{z}$-rotation that is shown in inset of Fig. 2A (see Appendix B for details). The general instability analysis was implemented on this simulated data, after converting the 4-outcome data to the standard “success”/“fail” format of RB. This analysis yielded a time-dependent success probability for each circuit. Following our non-intrusive framework, instantaneous success probabilities at each time of interest were then fed into the standard RB data analysis (fitting an exponential) as shown for three times in Fig. 2B. The instantaneous RB error rate estimate is then (up to a constant Proctor et al. (2019)) the decay rate of the fitted exponential at that time. The resultant time-resolved RB error rate is shown in Fig. 2A. It closely tracks the true error rate.

GST is a method for high-precision tomographic reconstruction of a set of time-independent gates, state preparations and measurements Blume-Kohout et al. (2017); Merkel et al. (2013). We consider GST on a gate set comprised of standard $\hat{z}$-axis preparation and measurement, and three gates $\mathsf{G}_{\mathsf{x}}$, $\mathsf{G}_{\mathsf{y}}$, and $\mathsf{G}_{\mathsf{i}}$. Here $\mathsf{G}_{\mathsf{x/y}}$ are $\nicefrac{{\pi}}{{2}}$ rotations around the $\hat{x}/\hat{y}$ axes and $\mathsf{G}_{\mathsf{i}}$ is the idle gate. The GST circuits have the form $\mathsf{S}_{\text{prep}}\mathsf{S}_{\text{germ}}^{k}\mathsf{S}_{\text{meas}}$ (circuits are written in operation order where the leftmost operation occurs first). In this circuit: $\mathsf{S}_{\text{prep}}$ and $\mathsf{S}_{\text{meas}}$ are each one of six short sequences chosen to generate tomographically complete state preparations and measurements; $\mathsf{S}_{\text{germ}}$ is one of twelve short “germ” sequences, chosen so that powers (repetitions) of these germs amplify all coherent, stochastic and amplitude-damping errors; $k$ runs over an approximately logarithmically spaced set of integers, given by $k=\lfloor L/|\mathsf{S}_{\text{germ}}|\rfloor$ where $|\mathsf{S}_{\text{germ}}|$ is the length of the germ and $L=2^{0},2^{1},2^{2},\dots,L_{\max}$ for some maximum germ power $L_{\max}$.

We simulated data from 1000 rasters through these GST circuits (with $L_{\max}=128$). The error model consisted of 0.1% depolarization on each gate. Additionally, $\mathsf{G}_{\mathsf{x}}$ and $\mathsf{G}_{\mathsf{y}}$ are subject to over/under-rotation errors that oscillate both quickly and slowly, while $\mathsf{G}_{\mathsf{i}}$ is subject to slowly varying $\hat{z}$-axis coherent errors. We used our intrusive approach to time-resolved tomography: the general instability analysis was implemented on this simulated data, the results were used to select a time-resolved model for the gates, and this model was then fit to the time-series data using maximum likelihood estimation (see Appendix B for details). The resulting time-resolved estimates of the gate rotation angles are shown in Fig. 2C-D. The estimates closely track the true values.

Demonstration on experimental data. Having verified that our methods are compatible with data from GST circuits, we now demonstrate time-resolved GST on two sets of experimental data, using the three gates $\mathsf{G}_{\mathsf{x}}$, $\mathsf{G}_{\mathsf{y}}$ and $\mathsf{G}_{\mathsf{i}}$. These experiments comprehensively quantify the stability of our 171Yb*+* qubit, because the GST circuits are tomographically complete and they amplify all standard types of error in the gates. The $\mathsf{G}_{\mathsf{x}}$ and $\mathsf{G}_{\mathsf{y}}$ gates were implemented with BB1 compensated pulses Wimperis (1994); Merrill and Brown (2012), and $\mathsf{G}_{\mathsf{i}}$ was implemented with a dynamical decoupling $X_{\pi}Y_{\pi}X_{\pi}Y_{\pi}$ sequence Khodjasteh and Viola (2009), where $X_{\pi}$ and $Y_{\pi}$ represent $\pi$ pulses around the $\hat{x}$ and $\hat{y}$ axes. The first round of data collection included the GST circuits to a maximum germ power of $L_{\text{max}}=2048$ (resulting in 3889 circuits). These circuits were rastered 300 times over approximately $5.5$ hours.

Fig. 3A-B summarizes the results of our general instability assessment on this data, using a representation that is tailored to GST circuits. Each pixel in this plot corresponds to a single circuit, and summarizes the evidence for instability by $\lambda_{\mathsf{p}}=-\log_{10}(\mathsf{p})$ where $\mathsf{p}$ is the p-value of the largest power in the spectrum for that circuit ($\lambda_{\mathsf{p}}$ is 5% significant when it is above the multi-test adjusted threshold $\lambda_{\mathsf{p},\text{threshold}}\approx 7$). The only circuits that displayed detectable instability are those that contain many sequential applications of $\mathsf{G}_{\mathsf{i}}$. Fig. 3B further narrows this down to generalized Ramsey circuits, whereby the qubit is prepared on the equator of the Bloch sphere, active idle gates are applied, and then the qubit is measured on the equator of the Bloch sphere. These circuits amplify erroneous $\hat{z}$-axis rotations in $\mathsf{G}_{\mathsf{i}}$. Other GST circuits amplify all other errors, but none of those circuits exhibit detectable drift. This is conclusive evidence that the angle of these $\hat{z}$-axis rotations is varying over the course of the experiment.

The instability in $\mathsf{G}_{\mathsf{i}}$ can be quantified by implementing time-resolved GST, with the $\hat{z}$-axis error in $\mathsf{G}_{\mathsf{i}}$ expanded into a summation of Fourier coefficients (see Appendix B for details). The results are summarized in Fig. 3C-D (dotted lines). Fig. 3D shows the diamond distance error rate ($\epsilon_{\diamond}$) Aharonov et al. (1998) in the three gates over time. It shows that $\mathsf{G}_{\mathsf{i}}$ is the worst performing gate, and that the error rate of $\mathsf{G}_{\mathsf{i}}$ drifts substantially over the course of the experiment ($\epsilon_{\diamond}$ varies by \sim$$25\%). The gate infidelities are an order of magnitude smaller (see Table 1 in Appendix B). Fig. 3C shows the coherent component of the $\mathsf{G}_{\mathsf{i}}$ gate over time, resolved into rotation angles $\theta_{x}$, $\theta_{y}$ and $\theta_{z}$ around the three Bloch sphere axes $\hat{x}$, $\hat{y}$ and $\hat{z}$. The varying $\hat{z}$-axis component is the dominant source of error.

This first round of experiments revealed instability, so we changed the experimental setup. Changes included the addition of periodic recalibration of the microwave drive frequency, the $\pi$-pulse duration, and the pointing of the detection laser (details in the Methods). We then repeated this GST experiment. To increase sensitivity to any instability, we collected more data, over a longer time period, and we included longer circuits. We ran the GST circuits out to a maximum germ power of $L_{\text{max}}=16384$, rastering $328$ times through this set of 5041 circuits over approximately 40 hours. The purpose of running such a comprehensive experiment was to maximize sensitivity — our methods need much fewer experimental resources for useful results (see below). Repeating the above analysis on this data, we found that none of the $\lambda_{\mathsf{p}}$ were statistically significant, i.e., no instability was detected in any circuit, including circuits containing over $10^{5}$ sequential $\mathsf{G}_{\mathsf{i}}$ gates. Again, we performed time-resolved GST. Since no time dependence was detected, this reduces to standard time-independent GST. The results are summarized in Fig. 3C-D (unbroken lines). The gate error rates have been substantially suppressed ($\epsilon_{\diamond}$ decreased by \sim$$10\times for $\mathsf{G}_{\mathsf{i}}$), and the $\hat{z}$-axis coherent error in $\mathsf{G}_{\mathsf{i}}$ reduced and stabilized. This is a comprehensive demonstration that the recalibrations are stabilizing the qubit. Furthermore, the recalibrated parameters versus time are strongly correlated with ambient laboratory temperature (see the Methods), suggesting temperature stabilization as an alternative route to qubit stabilization, and supporting the conclusions of our Ramsey experiments.

No individual circuit exhibited signs of drift in this second GST experiment, but we can also perform a collective test for instability on the clickstreams from all the circuits. In particular, we can average the per-circuit power spectra, and look for statistically significant peaks in this single spectrum. This suppresses the shot noise inherent in each individual clickstream, so it can reveal low-power drift that would otherwise be hidden in the noise (see Appendix A). This average spectrum is shown for both experiments in Fig. 3E. The power at low frequencies decreases substantially from the first to the second experiment, further demonstrating that our drift compensation is stabilizing the qubit. However, there is power above the $5\%$ significance threshold for both experiments. So there is still some residual instability after the experimental improvements. But this residual drift is no longer a significant source of errors, as demonstrated by the low and stable error rates shown in Fig. 3D.

Experiment design. Our method is an efficient way to identify time dependence in the outcome probability distribution of any quantum circuit. In its most basic application, it can verify the stability of application or benchmarking data. No special-purpose circuits are required, as the drift detection can be applied to data that is already being taken. The analysis will then be sensitive to any drifting errors that impact this application, in proportion to their effect on the application.

As we have demonstrated, our method can also be used to create dedicated drift characterization protocols. This mode requires a carefully chosen set of quantum circuits that are sensitive to the specific parameters under study. Without a priori knowledge about what may be drifting, this circuit set should be sensitive to all of the parameters of a gate set. The GST circuits are a good choice. However, if only a few parameters are expected to drift, a smaller set of circuits sensitive only to these parameters can be used, resulting in a more efficient experiment. For example, Ramsey circuits serve as excellent probes of time variation in qubit phase rotation rates. Many of the most sensitive circuits, such as those used in GST, Ramsey spectroscopy, and robust phase estimation Kimmel et al. (2015), are periodic and extensible. These circuits achieve $\mathcal{O}(1/L)$ precision scaling, with $L$ the maximum circuit length, up until decoherence dominates. So, by choosing a suitably large $L$, very high-precision drift tracking can be achieved, as in our experiments.

Interleaving dedicated drift characterization circuits with application circuits combines the two use cases for our methods — dedicated drift characterization and auxiliary analysis. This reduces the data acquisition rate for both the application and characterization circuits, but it directly probes whether time variation in a parametric model is correlated with drift in the outcomes of an application circuit. While this reduces sensitivity to high-frequency instabilities, much of the drift seen in the laboratory is on timescales that are long compared to the data acquisition rate. As a simple demonstration of this, we note that discarding 80% of our Ramsey data — keeping only every fifth bit for each circuit — still yields a high-precision time-resolved phase estimate, as shown in Fig. 1E (gray dashed line).

The sensitivity of our analysis depends on both the number of times a circuit is repeated ($N$) and the sampling rate ($t_{\rm gap}$). As in all signal analysis techniques, the sampling rate sets the Nyquist limit — the highest frequency the analysis is sensitive to without aliasing — while $(N-1)t_{\rm gap}$ sets the lowest frequency drift that will be visible. While the sensitivity of our methods increases with more data, statistically significant results can be achieved without dedicating hours or days to data collection. For example, both the simulated GST and RB experiments (see Fig. 2) used a number of circuits and repetitions consistent with standard practices. Further details relating the sampling parameters and the analysis sensitivity are provided in Appendix C.

Discussion

Reliable quantum computation demands stable hardware. But current standards for characterizing QIPs assume stability — they cannot verify that a QIP is stable, nor can they quantify any instabilities. This is becoming a critical concern as stable sources of errors are steadily reduced. For example, drift significantly impacted the recent tomographic experiments of Wan et al. Wan et al. (2019) but this was only verified using a complex analysis. In this article we have introduced a general, flexible, and powerful methodology for diagnosing instabilities in a QIP. We have applied these methods to a trapped-ion qubit, demonstrating both time-resolved phase estimation and time-resolved tomographic reconstructions of logic gates. Using these tools, we were able to identify the most unstable gate, confirm that periodic recalibration stabilized the qubit to an extent that drift is no longer a significant source of error, and isolate the probable source of the instabilities (temperature changes).

Our methods are widely-applicable, platform-independent, and do not require special-purpose experiments. This is because the core techniques are applicable to data from any set of quantum circuits — as long as it is recorded as a time series — and the data analysis is fast and simple (speed is limited only by the fast Fourier transform). These techniques enable routine stability analysis on data gathered primarily for other purposes, such as data from algorithm, benchmarking or error correction circuits. These techniques are even applicable outside of the context of quantum computing — they could be used for time-resolved quantum sensing. We have incorporated these tools into an open-source software package Nielsen et al. (2020a, b), making it easy to check any time-series QIP data for signs of instability. Because of the disastrous impact of drift on characterization protocols Dehollain et al. (2016); Epstein et al. (2014); van Enk and Blume-Kohout (2013); Fong and Merkel (2017); Fogarty et al. (2015); Chow et al. (2009); Wan et al. (2019), its largely unknown impact on QIP applications, and the minimal overhead required to implement our methods, we hope to see this analysis broadly and quickly adopted.

Methods

Experiment details. We trap a single 171Yb*+* ion $\sim$34\text{,}\mathrm{\SIUnitSymbolMicro m} above a Sandia multi-layer surface ion trap with integrated microwave antennae, shown in Fig. [**5**](#A1.F5). The radial trapping potential is formed with $170\text{\,}\mathrm{V}$ of rf-drive at $88\text{\,}\mathrm{MHz}$; the axial field is generated by up to $2\text{\,}\mathrm{V}$ on the segmented dc control electrodes. This yields secular trap frequencies of $0.7\text{\,}\mathrm{MHz}$, $5\text{\,}\mathrm{MHz}$, and $5.5\text{\,}\mathrm{MHz}$ for the axial and radial modes, respectively. An electromagnetic coil aligned with its axis perpendicular to the trap surface creates the quantization field of approximately $5\,\textrm{G}$ at the ion. The field magnitude is calibrated using the qubit transition frequency, which has a second-order dependence on the magnetic field of $f=$12.642\,812\,118\text{\,}\mathrm{GHz}$+310.8B^{2}$\text{\,}\mathrm{Hz}, where $B$ is the externally applied magnetic field in Gauss Fisk et al. (1997). The qubit is encoded in the hyperfine clock states of the ${}^{2}\textrm{S}_{1/2}$ ground state of 171Yb*+*, with logical 0 and 1 defined as ${\left|F=0,m_{F}=0\right\rangle}$ and ${\left|F=1,m_{F}=0\right\rangle}$, respectively.

Each run of a quantum circuit consists of four steps: cooling the ion, preparing the input state, performing the gates, and then measuring the ion. First, using an adaptive length Doppler cooling scheme, we verify the presence of the ion. The ion is Doppler cooled for $1\text{\,}\mathrm{ms}$, during which fluorescence events are counted. If the number of detected photons is above a threshold ($\sim 85\%$ of the average fluorescence observed for a cooled ion) Doppler cooling is complete, otherwise, the cooling is repeated. If the threshold is not reached after 300 repetitions, the experiment is halted to load a new ion. This ensures that an ion is present in the trap and that it is approximately the same temperature for each run. After cooling and verifying the presence of the ion, it is prepared in the ${\left|F=0,m_{F}=0\right\rangle}$ ground state using an optical pumping pulse Olmschenk et al. (2007). All active gates are implemented by directly driving the $12.6428\text{\,}\mathrm{GHz}$ hyperfine qubit transition, using a near-field antenna integrated into the trap (see Fig. 5). The methods used for generating the microwave radiation are discussed in Ref. Blume-Kohout et al. (2017). A standard state fluorescence technique Olmschenk et al. (2007) is used to measure the final state of the qubit.

The gates we use are $\mathsf{G}_{\mathsf{x}}$, $\mathsf{G}_{\mathsf{y}}$ and $\mathsf{G}_{\mathsf{i}}$, which are $\nicefrac{{\pi}}{{2}}$ rotations around the $\hat{x}$- and $\hat{y}$-axes, and an idle gate. The $\mathsf{G}_{\mathsf{x}}$ and $\mathsf{G}_{\mathsf{y}}$ gates, used in both the Ramsey and GST experiments, are implemented using BB1 pulse sequences Wimperis (1994); Merrill and Brown (2012). The $\mathsf{G}_{\mathsf{i}}$ gate, used only in the GST experiments, is a second-order compensation sequence: $\mathsf{G}_{\mathsf{i}}=X_{\pi}Y_{\pi}X_{\pi}Y_{\pi}$, where $X_{\pi}$ and $Y_{\pi}$ denote $\pi$ pulses about the $\hat{x}$- and $\hat{y}$-axis, respectively Khodjasteh and Viola (2009). To maintain a constant power on the microwave amplifier and reduce the errors from finite on/off times, active gates are performed gapless, i.e., we transition from one pulse to the next by adjusting the phase of the microwave signal without changing the amplitude of the microwave radiation. In the first GST experiment the Rabi frequency was $119\text{\,}\mathrm{kHz}$.

Changing the phase in the analog output signal takes $\sim$5\text{,}\mathrm{ns}$$ and this causes errors, because the pulse sequences are performed gapless. These errors are larger for shorter $\pi$-pulse times. To reduce this error, in this second GST experiment the Rabi frequency was decreased to $74\text{\,}\mathrm{kHz}$. To compensate for the drift that we observed in both the Ramsey experiment and the first GST experiment, in the second GST experiment we incorporated three forms of active drift control. The detection laser position was recalibrated every 45 minutes, and both the $\pi$-time ($\tau_{\pi}$) and the microwave drive frequency ($f$) were updated based on the results of interleaved calibration circuits. After every $4^{\rm th}$ circuit, a circuit consisting of a $10.5\pi$ pulse was performed. If the outcome was 0 (resp., 1) then $1.25\text{\,}\mathrm{ns}$ was added to $\tau_{\pi}$ (resp., subtracted from $\tau_{\pi}$). The applied $\pi$-time is $\tau_{\pi}$ rounded to an integer multiple of $20\text{\,}\mathrm{ns}$, so only consistently bright or dark measurements result in changes of the pulse time. After every $16^{\rm th}$ circuit a $10\text{\,}\mathrm{ms}$ wait Ramsey circuit was performed. If the outcome was 0 (resp., 1) $10\text{\,}\mathrm{mHz}$ is added to $f$ (resp., subtracted from $f$).

Fig. 4 shows the detection beam position, $\tau_{\pi}$, $f$, and ambient laboratory temperature over the course of the second GST experiment. The calibrated $f$ is correlated with the ambient temperature. This is consistent with the observed correlation between the ambient temperature and the estimated detuning in the Ramsey experiment (see Fig. 1E). The temperature is also strongly correlated with the calibrated detection beam location points, suggesting that thermal expansion is a plausible underlying cause of the frequency shift.

Data analysis details. To generate a power spectrum from a clickstream we use the Type-II discrete cosine transform with an orthogonal normalization. This is the matrix $F$ with elements

[TABLE]

where $\omega,i=0,\dots,N-1$ Ahmed et al. (1974). However, note that the exact transform used is not important: we only require that $F$ is an orthogonal and Fourier-like matrix (see Appendix A). Our hypothesis testing is all at a statistical significance of 5%, and uses a Bonferroni correction to maintain this significance when implementing many hypothesis tests (see Appendix A). All data fitting uses maximum likelihood estimation, except for the $p(t)$ estimation in the time-resolved RB simulations. In that case, we use a simple form of signal filtering (see Appendix A), so that the entire analysis chain maintains the speed and simplicity inherent to RB. When choosing between multiple time-resolved models, as in the time-resolved Ramsey tomography and GST analyses, we use the Akaike information criteria Akaike (1974) to avoid overfitting (see Appendix B). Further details on these methods, and a supporting theory, is provided in Appendices A-C.

Data availability

All experimental and simulated data presented in this paper are available at http://doi.org/10.5281/zenodo.4033077.

Code availability

Code for implementing the general drift charactierization methods introduced in this paper has been incorporated into the open-source Python package pyGSTi Nielsen et al. (2020a, b). The pyGSTi-based Python scripts and notebooks used for the data analysis reported in this paper are available at http://doi.org/10.5281/zenodo.4033077.

Acknowledgments

This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research Quantum Testbed Program; the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA); and the Laboratory Directed Research and Development program at Sandia National Laboratories. Sandia National Laboratories is a multi-program laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525. All statements of fact, opinion or conclusions contained herein are those of the authors and should not be construed as representing the official views or policies of IARPA, the ODNI, the U.S. Department of Energy, or the U.S. Government.

Appendix A Drift characterization using time-series data

Here we provide supporting theory and information for the methods that we propose for detecting and characterizing instability in the outcome distribution of arbitrary quantum circuits. This is the methods summarized in steps 1-3 of the flowchart in Fig. 1.

A.1 Power spectra

To begin, we detail how a clickstream is transformed into the frequency domain and then converted into a power spectrum (step 1 in the flowchart of Fig. 1A). To transform into the frequency domain we apply a Fourier-like transform $F$ to rescaled data. All of the data analysis in this article uses the Type-II discrete cosine transform (DCT), with an orthogonal normalization. It is define by the matrix $F^{(\textsc{dct})}$ with the elements

[TABLE]

where $\omega,i=0,\dots,N-1$ Rao and Yip (2014); Ahmed et al. (1974). The Type-II DCT typically results in a sparser representation of signals with differing values at the boundaries than many of the alternative Fourier-like transforms Rao and Yip (2014); Ahmed et al. (1974). Functions of this sort are typical in our problem, and so this is our motivation for using this transform.

Although all our data analysis uses the Type-II DCT, our methods and our theory (except where stated) are all directly applicable for any transform matrix $F$ that satisfies the following criteria: (i) $F$ is orthogonal, (ii) every matrix element of $F$ satisfies $F_{\omega i}^{2}\leq b_{F}/N$ for a small constant $b_{F}$, and (iii) the top row of $F$ is proportional to the vector of all ones. This includes, e.g., sine transforms, other cosine transforms, and the Walsh-Hadamard transform Feldman (1987); Rao and Yip (2014); Ahmed et al. (1974) (all with the appropriate normalizations). Note that the standard discrete Fourier transform does not satisfy these requirements (with an appropriate normalization it is unitary, but not orthogonal). However, the methods described here are compatible with the discrete Fourier transform as long as some minor adjustments are made to guarantee that, e.g., the probability trace estimates are real. Throughout the following theory, we refer to the indices of the $F$-domain vector as frequencies, but note that (i) when we report frequencies for experimental data these indices have been converted to Hertz, and (ii) for a generic transformation an index should only be interpreted as indicating the corresponding basis function of $F$, and it may or may not constitute a frequency.

As discussed in the main text, the motivation for using spectral analysis is that many physically plausible $p(t)$ are sparse in the frequency domain. But note that the relevant quantity to the data analysis is the discrete-time probability trajectory $\mathbf{p}$, where $p_{i}=p(t_{i})$ and $t_{i}$ is the $i^{\rm th}$ data collection time for the circuit in question. Moreover, as we are concerned with deviations from time-independence, it is the “signal” vector

[TABLE]

that is of most relevance, where $\bar{v}$ denotes the mean of the vector $\mathbf{v}$ and $\mathbf{1}=(1,1,\dots,1)^{\mathsf{T}}$. It is sparsity of $\tilde{\mathbf{s}}$ that is important for our data analysis. Sparsity in the frequency domain of the continuous-time probability $p(t)$ implies that $\tilde{\mathbf{s}}=F\mathbf{s}$ is a sparse vector if: (i) a sensible choice is made for the transformation $F$, such as the Type-II DCT, (ii) the number of sample times $N$ is sufficiently high, and (iii) the sample times $t_{0},t_{1},t_{2},\dots,t_{N-1}$ for the circuit in question are spaced reasonably uniformly.

All the experiments and simulations herein consist of rastering through a set of circuits — i.e., running each circuit once, and then looping through this until each circuit has been run $N$ times. So the time between sequential runs of a circuit is exactly constant (and the same for all the circuits) in an idealized experiment. Violations of this ideal are unavoidable, but if they are small (as they are in our experiments) they only mildly degrade the effectiveness of our techniques. This can be understood by noting that small perturbations on equally-spaced times (or a small number of erroneously long time-steps of any size) will only slightly reduce the sparsity of $\tilde{\mathbf{s}}$, for any Fourier-sparse $p(t)$. Large deviations from the equally-spaced-times ideal can substantial degrade the sensitivity of the methods we presented in the main text (although this will not cause instability to be spuriously detected). In this case, we can instead generate power spectra using a technique that explicitly takes into account the sample times. One such option is the floating-mean Lomb-Scargle periodogram Zechmeister and Kürster (2009); VanderPlas (2018), and all our methods can be adapted to this spectral analysis technique. As there are technical differences in this case, and this method is not needed to analyze our experimental data, we do not include the details here.

Before transforming a clickstream ($\mathbf{x}$) to the frequency domain, we first subtract its mean and divide by its variance. This simplifies the statistics of the power spectrum, but the spectrum is not well-defined if the clickstream consists of only zeros or only ones (because then the variance of $\mathbf{x}$ is zero). So we define the frequency-domain amplitudes by:

[TABLE]

The power at frequency $\omega$ is then given by $|\tilde{z}_{\omega}|^{2}$. This convention for the constant clickstream case is convenient, particularly when averaging the power spectra from multiple clickstreams, as there are then fewer situations in which we need to explicitly take this special case into account. Moreover, $(0,1,1,\dots)$ roughly corresponds to the $N\to\infty$ limit of $\tilde{\mathbf{z}}$ for a clickstream containing a single 1 or 0 (the correspondence is not precise, because the spectrum of a delta function is not completely flat for a general transform $F$).

A.2 Statistical hypothesis testing

All of our methods are built upon statistical hypothesis testing Lehmann and Romano (2006); Shaffer (1995). We now review the relevant aspects of this field and formulate our hypothesis testing problem. Consider some set of random variables $\{A_{1},A_{2},\dots\}$ with a joint distribution that is parameterized by $\theta\in\mathcal{H}$ for some space $\mathcal{H}$. A statistical hypothesis is the conjecture that $\theta\in\mathcal{H}_{0}\subset\mathcal{H}$. In this article, data consists of a clickstream of $N$ bits for each of $C\geq 1$ circuits. Generalizing the notation used in the main text so that we explicitly denote the circuit index, let $x_{c,i}$ be the bit output by the $i^{\rm th}$ repetition of the $c^{\rm th}$ circuit, and denote the corresponding Bernoulli random variable from which that datum is drawn by $X_{c,i}$. The set of random variables about which we state statistical hypotheses is $\{X_{c,i}\}$. The random variable $X_{c,i}$ has some unknown probability $p_{c,i}$ of taking the value 1, so the set of random variables $\{X_{c,i}\}$ is parameterized by the space

[TABLE]

Statistical hypothesis testing starts from some set of null hypotheses $\{\mathcal{H}_{0,i}\}$ and uses statistical hypothesis tests on data drawn from the random variables to attempt to reject one or more of these null hypotheses (and/or intersections of these null hypotheses). Our methods look for evidence that $|\tilde{p}_{\omega,c}|>0$ at each non-zero frequency $\omega$ and for each circuit $c$. This is formalized by the set of null hypotheses $\{\mathcal{H}_{0,c,\omega}\}$ where $\mathcal{H}_{0,c,\omega}$ is the conjecture that $\tilde{p}_{c,\omega}=0$, with $\omega=1,2,\dots,N-1$ and $c=1,2,\dots,C$. For our purposes, testing a null hypothesis consists of the following four steps: (i) pick a significance level $\alpha\in(0,1)$, with $\alpha=5\%$ a common choice; (ii) select a test statistic $\Lambda$, which is a function from data to $\mathbb{R}$; (iii) find a threshold such that the probability of observing $\Lambda$ larger than this threshold if the null hypothesis is true is at most $\alpha$; (iv) collect data, evaluate $\Lambda$, and reject the null hypothesis if and only if $\Lambda$ is larger than the threshold. This procedure guarantees that the null hypothesis is falsely rejected with probability at most $\alpha$. Our test statistic for testing the null hypothesis $\mathcal{H}_{0,c,\omega}$ is simply the power of the normalized data at the frequency $\omega$, i.e., $|\tilde{z}_{c,\omega}|^{2}$.

Implementing multiple hypothesis tests that use the above procedure will typically cause the probability of falsely rejecting at least one null hypothesis — the family-wise error rate (FWER) — to increase quickly with the number of tests (for $T$ tests $\text{FWER}\leq T\alpha$). The standard approach when implementing multiple tests is to adapt the testing procedure to, at a minimum, maintain weak control of the FWER. Weak control of the FWER means that if all the null hypotheses are true then, for some pre-specified global significance $\alpha$, the probability of falsely rejecting one or more null hypotheses is at most $\alpha$. It is common to choose a testing procedure that also maintains strong control of the FWER. Strong control of the FWER means that the probability of falsely rejecting one or more true null hypotheses is at most $\alpha$ and this holds even if one or more of the null hypotheses are false. We seek to maintain strong control of the FWER, so that whenever we reject $\mathcal{H}_{0,c,\omega}$ we can be $(1-\alpha)$ confident that $\omega$ is a frequency with a non-zero amplitude in the Fourier decomposition of $\mathbf{p}_{c}$.

Strong control of the FWER at a global significance of $\alpha$ can be maintained by implementing the $i^{\rm th}$ of $T$ tests at a local significance of $\alpha w_{i}$ for any $w_{i}\geq 0$ satisfying $\sum_{i}w_{i}=1$. This is known as the (generalized) Bonferroni correction Lehmann and Romano (2006); Shaffer (1995). The Bonferroni correction can be conservative; there are other methods for maintaining strong control of the FWER that are more powerful than the Bonferroni correction, such as the Holms procedure Lehmann and Romano (2006); Shaffer (1995). Moreover, strong control of the FWER is not the only way to ensure “confidence” in the results — e.g., we could instead choose to control the false discovery rate, which is the expected ratio of the number of falsely rejected null hypotheses to the total number of rejected null hypotheses Benjamini and Hochberg (1995). Variants of our techniques that control the false discovery rate and/or use a different multi-test correction procedure are possible, and essentially involve setting a slightly lower power significance threshold in the analysis.

A.3 Statistics of the power spectra

We now derive the distribution of the power spectrum of the data, both when $p(t)$ is and is not constant. We consider the frequency-domain data $\tilde{z}_{\omega}$, defined in Eq. (5), for an arbitrary non-zero frequency $\omega$, and we derive an approximation to the distribution of the corresponding random variable,

[TABLE]

where $\omega>0$ (we are again dropping the circuit indexing, and we have ignored the edge-case of the constant clickstream instance in this equation). Note that $\bar{X}=\frac{1}{N}\sum_{i}X_{i}$ is the average of random variables, and so it is itself a random variable. In the main text we asserted that $\tilde{Z}^{2}_{\omega}$ is approximately $\chi_{1}^{2}$ distributed whenever the probability trajectory is a constant. This is a sub-case of what we derive below, but we also study the marginal distribution of $\tilde{Z}_{\omega}^{2}$ when $\tilde{p}_{\omega}=0$ but the probability trajectory is not necessarily constant (i.e., there is perhaps power at one or more other non-zero frequencies). We do so in order to be able to claim that we are maintaining strong control of the FWER. By itself, the fact that $\tilde{Z}_{\omega}^{2}$ is approximately $\chi_{1}^{2}$ distributed when the probability trajectory is constant can only be used to ensure weak control of the FWER.

The random variable $\tilde{Z}_{\omega}$ is rather complex to directly analyze, because it depends on both $\tilde{X}_{\omega}$ and $\bar{X}$. So we first make an approximation. As $\bar{X}$ is the mean of the $X_{i}$ we have that $\mathbb{E}[\bar{X}]=\bar{p}$ and $\mathbb{V}[\bar{X}]=O(\nicefrac{{1}}{{N}})$, where $\mathbb{E}(A)$ and $\mathbb{V}(A)$ denote the expectation value and variance of the random variable $A$, respectively. Therefore

[TABLE]

where we assume that $0<\bar{p}<1$. So we now derive the distribution of $\tilde{Y}_{\omega}$, as an approximation to the distribution of $\tilde{Z}_{\omega}$.

The random variable $\tilde{Y}_{\omega}$ is the sum of $N$ independent random variables:

[TABLE]

A heuristic appeal to the central limit theorem therefore suggests that $\tilde{Y}_{\omega}$ is normally distributed. That is, $\tilde{Y}_{\omega}$ is approximately $\mathcal{N}(\mu_{\omega},\nu_{\omega})$ distributed for $N\gg 1$, with some mean $\mu_{\omega}$ and variance $\nu_{\omega}$. Later we formally apply the central limit theorem to $\tilde{Y}_{\omega}$, but first we derive formulae for $\mu_{\omega}$ and $\nu_{\omega}$. For any non-zero $\omega$, $\mu_{\omega}$ is simply a rescaling of the component of the probability trajectory at this frequency:

[TABLE]

So, $\mu_{\omega}=0$ under the null hypothesis that $\tilde{p}_{\omega}=0$. The variance $\nu_{\omega}$ has a more subtle dependence on the probability trajectory. By noting that $\tilde{Y}_{\omega}$ is a summation of independent Bernoulli random variables multiplied by constants, and by noting that the variance of a Bernoulli random variable with bias $b$ is $b(1-b)$, we obtain

[TABLE]

It then follows from the orthogonality of $F$ that $\nu_{\omega}=1$ when the probability trajectory is constant. So $\tilde{Z}_{\omega}$ is approximately $\mathcal{N}(0,1)$ distributed when the probability trajectory is constant. Therefore, calculating significance thresholds using the model that $\tilde{Z}_{\omega}^{2}$ is $\chi^{2}_{1}$ distributed under the null hypothesis will maintain weak control of the FWER.

Using the relation $p_{i}=s_{i}+\bar{p}$, where $\mathbf{s}$ is the signal vector, we can rewrite the variance as

[TABLE]

where

[TABLE]

So far, $F$ is a general Fourier-like transformation that is an orthogonal matrix (see earlier). For some transformation, such as the Walsh-Hadamard transformation (but not the DCT), $F_{\omega i}^{2}=1/N$ and so Eq. (13) simplifies to

[TABLE]

because $\sum_{i}s_{i}=0$. Here

[TABLE]

is the 2-norm. That is, the variance is upper-bounded by 1, and it decreases in proportion to the signal power per time-step (which is equal to the variance of $\mathbf{p}$). So the maximum variance over the null hypothesis space of $\tilde{p}_{\omega}=0$ is $\nu_{\omega}=1$. This means that modeling $\tilde{Z}_{\omega}^{2}$ with a $\chi_{1}^{2}$ distribution is sufficient for strong control of the FWER if the power spectra are generated using, e.g., the Walsh-Hadamard transformation.

For a general transformation the situation is more complicated. In particular, even with $p_{\omega}=0$ it is possible to obtain $\nu_{\omega}>1$ for a matrix $F$ containing elements of magnitude greater than $\nicefrac{{1}}{{\sqrt{N}}}$. Because all of our data analysis uses the Type-II DCT, defined in Eq. (3), at this point we specialize to this transformation. In this case, simple algebra can be used to show that

[TABLE]

where $\mathbf{q}=(s_{0}^{2},s_{1}^{2},s_{2}^{2},\dots)$ is the vector of time-domain signal powers, and we are allowing the frequency index to extend outside its range as specified in the definition of the DCT. This equation implies that the variance at $\omega$ is increased by signal power at $2\omega$ if the oscillations are not centered on $\nicefrac{{1}}{{2}}$. This can only result in a variance above 1 if the amplitude at $2\omega$ outweighs the decrease in the variance caused by the total power per time-step in $\mathbf{s}$. So the variance at $\omega$ is maximized by a pure-tone probability trajectory oscillating between $\nicefrac{{1}}{{2}}$ and 1 or 0 at a frequency of $2\omega$, resulting in the saturable bound $\nu_{\omega}\leq 1+\nicefrac{{1}}{{6}}$. Therefore the maximal variance over the null hypothesis space of $\tilde{p}_{\omega}=0$ is $\nu_{\omega}=\nicefrac{{7}}{{6}}$. So, ignoring the two approximations made so far, there exists a probability trajectory within the null hypothesis subspace defined by $\tilde{p}_{\omega}=0$ such that $\tilde{Z}_{\omega}$ is $\mathcal{N}(0,\nicefrac{{7}}{{6}})$ distributed. This implies that calculating significance thresholds using a $\chi_{1}^{2}$ does not strictly provide strong control of the FWER.

We could guarantee strong control of the FWER by calculating thresholds using a rescaled $\chi_{1}^{2}$ distribution (or, equivalently, by normalizing the data differently). But we choose not to do this, because for a broad class of probability trajectories $\nu_{\omega}\leq 1$ for all $\omega$. So accounting for the possibility of $\nu_{\omega}>1$ would result in a reduction in test power — i.e., higher significance thresholds — to gain statistical rigor in some edge cases. In our opinion this is not a good tradeoff, and we instead settle for maintaining weak control of the FWER and “almost” maintaining strong control of the FWER. This is true in the sense that although there are some probability trajectories for which the FWER is slightly above the desired value, (i) the increase in the FWER is small, (ii) it can only occur for extremal points in the parameter space, and (iii) it is counteracted by the conservative nature of the Bonferroni correction that we use when calculating thresholds from a $\chi^{2}_{1}$ distribution (see later). All further statements about strong control of the FWER should be understood to implicitly contain this minor caveat.

The above reasoning was all based on the assumption that $\tilde{Y}_{\omega}$ is approximately $\mathcal{N}(\mu_{\omega},\nu_{\omega})$ distributed for large $N$. So, we now prove that the $N\to\infty$ limiting distribution of $(\tilde{Y}_{\omega}-\mu_{\omega})/\sqrt{\nu_{\omega}}$ is $\mathcal{N}(0,1)$, under a reasonable notion of the limit of a probability trajectory. In particular, we assume that $\epsilon_{1}<\bar{p}<1-\epsilon_{1}$ and $\nu_{\omega}>\epsilon_{2}$ for some constants $\epsilon_{1},\epsilon_{2}>0$ and all $N$. The proof will use Lyapunov’s central limit theorem Billingsley (1995), which we now state. Let $\{A_{1},A_{2},\dots,A_{n}\}$ be a sequence of independent random variables where $A_{i}$ has a finite expected value $\mu_{i}$ and variance $\nu_{i}$, and define $s_{n}^{2}=\sum_{i}\nu_{i}$. Then, if Lyapunov’s condition

[TABLE]

holds for some $\delta>0$, the distribution of $A=\sum_{i=1}^{n}(A_{i}-\mu_{i})/s_{n}$ converges to $\mathcal{N}(0,1)$ as $N\to\infty$. The random variable $\tilde{Y}_{\omega}$ may be written in the format of this central limit theorem: $\tilde{Y}_{\omega}=\sum_{i}Q_{\omega i}$ where $Q_{\omega i}=F_{\omega i}Y_{i}$, and the expected value and variance of every $Q_{\omega i}$ is bounded (due to our assumptions). It remains to prove that Lyapunov’s condition holds. We have that

[TABLE]

because $Q_{\omega i}=F_{\omega i}Y_{i}$ and the moments of $Y_{i}$ are bounded from above by a constant (as $\epsilon_{1}<\bar{p}<1-\epsilon_{1}$). Therefore

[TABLE]

with the equality holding because $s_{n}^{3}=\nu_{\omega}^{\nicefrac{{3}}{{2}}}$ and $\nu_{\omega}>\epsilon_{2}$. So Lyapunov’s condition holds for $\delta=1$.

In addition to performing hypothesis testing on the power spectrum for each circuit, we also test the single power spectrum obtained by averaging over the power spectra for the different circuits, e.g., see Fig. 3. The power at frequency $\omega$ in this spectrum is given by

[TABLE]

where $C$ is the number of circuits. The clickstreams for different circuits are independent, and so the corresponding random variable $\tilde{Z}_{\text{avg},\omega}^{2}$ is the average of independent random variables that are approximately $\chi_{1}^{2}$ distributed under the (intersection) null hypothesis

[TABLE]

that $\tilde{p}_{c,\omega}=0$, for all $c$. Therefore $\tilde{Z}_{\text{avg},\omega}^{2}$ is approximately $\nicefrac{{\chi^{2}_{C}}}{{C}}$ under this null hypothesis.

A.4 Hypothesis testing of the power spectra

Here we use the statistical model developed above to calculate the power significance threshold. We argued above that $\tilde{Z}_{c,\omega}^{2}$ can be modeled as a $\chi^{2}_{1}$ random variable under the null hypothesis that $\tilde{p}_{c,\omega}=0$ (up to minor caveats). So if we are only testing a single circuit $c$ and a single frequency $\omega$, then the $\alpha$-significance power threshold is simply

[TABLE]

where $\mathsf{CDF}_{k}$ denotes the cumulative distribution function for the $\chi^{2}_{k}$ distribution. However, we test every circuit at every non-zero frequency, and we also test the average power spectrum. So, to maintain strong control of the FWER with a global significance of $\alpha$ in this set of $(N-1)(C+1)$ tests, we use a Bonferroni correction. In general, this means setting the local significance of the test on circuit $c$ at frequency $\omega$ to $w_{c,\omega}\alpha$, and setting the local significance of the test on the average power spectrum at frequency $\omega$ to $w_{\omega}\alpha$ for some non-negative weightings $w_{c,\omega}$ and $w_{\omega}$ satisfying $\sum_{\omega}(w_{\omega}+\sum_{c}w_{c,\omega})=1$. The test weightings can be optimized for different usages.

Often it is equally important to identify any drift at every frequency and in every circuit. A one-parameter subclass of weightings that is useful under this circumstance is given by $w_{c,\omega}=(1-w)/(C(N-1))$ and $w_{\omega}=w/(N-1)$ for some weighting factor $w$. This choice of weightings results in a single $\omega$-independent significance threshold for the power in the individual power spectra for the $C$ circuits of

[TABLE]

and an $\omega$-independent significance threshold for the power in the averaged spectrum of

[TABLE]

The weighting factor $w$ controls the proportion of the test statistical significance that is allocated to the test on the averaged power spectrum versus the significance allocated to the individual power spectra. The averaged power spectrum does not contain information about which circuits (if any) exhibit drift, but has (under most circumstances) increased sensitivity to any drift. So the best choice for this weighting depends on the relative importance of detecting any drift versus being able to assigning it to a particular circuit. All our data analysis uses the thresholds given in Eqs. (24-25), with $\alpha=5\%$ and $w=\nicefrac{{1}}{{2}}$, with the exception of time-resolved RB, wherein $w=1$.

A.5 Estimating probability trajectories

We now explain our methods for estimating the circuit outcome probability trajectories (step 3, Fig. 1A). The methods are based on model selection, to avoid overfitting. The general form of the model that we use for a probability trajectory is:

[TABLE]

where $f_{\omega}(t)$ is the $\omega^{\rm th}$ basis function of the chosen Fourier transform, the summation is over some set of frequencies $\mathbb{W}$ that are to be chosen using model selection, and the $\gamma_{\omega}$ are parameters constrained only so that $p(\mathbf{\gamma},t)$ is a valid probability at all data acquisition times (the condition that $p(\mathbf{\gamma},t)$ is within $[0,1]$ for all real-valued $t$ can cause overly severe compression of the amplitudes). In the case of the Type-II DCT the basis functions are

[TABLE]

where, as throughout, $t_{i}$ is the $i^{\rm th}$ data collection time for the circuit in question, and

[TABLE]

is the size of the time-step.

The first step is model selection to identify a set of frequencies ($\mathbb{W}$). These frequencies can be chosen using the results of the hypothesis testing on the data spectra. It is useful to allow flexibility in exactly how they are chosen from the hypothesis testing results, because this allows the techniques to be adapted to different applications. A good general-purpose option, which results in estimates with a simple statistical meaning, is to set $\mathbb{W}$ to contain those and only those frequencies that are found to be statistically significant in the power spectrum for the circuit in question. With this method we are only including frequencies that contribute a nonzero component to the true $p(t)$ with confidence $(1-\alpha)$, where $\alpha$ is the chosen global significance level. All of the data analysis in this article uses this method, with the exception of the model selection in time-resolved RB (see later).

It remains to estimate the model parameters, i.e., the amplitudes $\mathbf{\gamma}$. Maximum likelihood estimation (MLE) is a statistically sound methodology. But MLE requires a numerical optimization, which can sometimes be slow. So we also provide an optimization-free method. The data amplitude $\tilde{x}_{\omega}$ is a well-motivated estimate of the probability trajectory amplitude at $\omega$, since $\mathbb{E}(\tilde{X}_{\omega})=\tilde{p}_{\omega}$. In terms of the rescaled amplitudes $\tilde{\mathbf{z}}$, this suggests setting $\gamma_{\omega}=S(\mathbf{x})\tilde{z}_{\omega}$ where $S(\mathbf{x})$ is a rescaling factor that inverts the normalization of the $\tilde{\mathbf{z}}$ amplitudes by the data standard deviation (and also rescales for the chosen normalization of the basis functions). This is intuitive but flawed: the resulting probability trajectory is not necessarily a valid probability at all times. To address this, we add a regularization to this estimator. We map each amplitude via

[TABLE]

with $\delta$ the smallest constant such that the estimated probability trajectory $p(\mathbf{\gamma},t)$ is within $[\epsilon,1-\epsilon]$ at all data acquisition times $t_{i}$, for some chosen $\epsilon\geq 0$ (e.g., $\epsilon=0$ if estimates on the parameter space boundary are considered acceptable). We refer to this estimation method as the “Fourier filter,” because it is a form of signal filtering. The Fourier filter returns very similar estimates to MLE under most conditions, although note that if $\epsilon=0$ it is possible for the Fourier filter to return estimates that have a likelihood of zero.

All data analysis in this paper uses MLE, except the time-resolved RB analysis. The motivation for using the Fourier filter in time-resolved RB is that then the speed of the analysis is limited only by the fast Fourier transform and curve fitting at the times of interest. This is particularly relevant to RB, as one motivation for using RB is its speed and simplicity. Moreover, note that for all the data herein the difference between the MLE and the Fourier filter estimates is negligible.

A.6 Analyzing data from circuits with more than two outcomes

Throughout this article we have only considered analysis of circuits with two possible outcomes. This is sufficient for all the experiments and simulations herein. However, the techniques we present can be easily generalized to $M$-outcome circuits, as we now briefly outline. In this more general setting, the aim is still to detect and quantify time variation in the probability distribution over outcomes for each circuit. However $M$ can be very large in commonly encountered circumstances, as $M=2^{n}$ for an $n$-qubit circuit terminated with the standard measurement. Even accurate estimation of a time-independent distribution over a very large number of outcomes is infeasible with reasonable quantities of data, so it will not be possible to accurately estimate a time-dependent distribution with many outcomes. For this reasons, it is useful to allow for a coarse-graining of the outcomes into $\check{M}\leq M$ bins (e.g., we can marginalize over some subset of the qubits in an $n$-qubit circuit or bin the outcomes into two sets). The analysis will then interrogate time-dependence in the coarse-grained probability distribution. So coarse-graining can make the analysis feasible with a reasonable amount of data, but the analysis will be insensitive to any time dependence that disappears under the chosen coarse-graining.

Each coarse-grained outcome has a bit-string time-series associated with it, corresponding to whether or not that measurement outcome was observed at each time. A power spectrum can be calculated for each such bit-string, and exactly the same hypothesis testing performed on each such power spectrum as in the 0/1 outcomes cases (with an appropriate correction for the increased number of hypothesis tests). Moreover, the spectrum obtained by averaging the spectra for a circuit over course-grained measurement outcomes can also be tested (which will be approximately $\chi^{2}_{\check{M}-1}$ distributed under the null hypothesis), or it can be tested instead of the per-outcome spectra — to avoid significance dilution by implementing too many hypothesis tests. The probability trajectory estimation (step 3, Fig. 1A) then proceeds exactly as in the 0/1 outcome case, except that now the estimation is of $\check{M}$ trajectories that must sum to 1 at all times. Any time-resolved parameter estimation proceeds exactly as in the 0/1 outcome case.

Appendix B Time-resolved tomography and benchmarking

Here we provide supporting theory and information for our time-resolved tomography and benchmarking methods. This is the methods summarized in steps 4 and 5 of the flowchart in Fig. 1.

B.1 Model selection in time-resolved tomography

We now explain the model selection step in our “intrusive” approach to time-resolved tomography and benchmarking (step 5a, Fig. 1A). Our methods for time-resolved tomography and benchmarking start from a time-independent $\{\gamma_{j}\}$-parameterized model $\mathcal{M}$ that predicts circuit outcomes. In the intrusive approach we expand each parameter $\gamma_{j}$ into the time-dependent form

[TABLE]

where the $\gamma_{j,\omega}$ are restricted only so that they satisfy any constraints on $\gamma_{j}$, the summation is over some set of non-zero frequencies $\mathbb{S}_{j}$, and $j$ indexes all the parameters in the time-independent model. So, we obtain a new model $\mathcal{M}_{\rm{tr}}$ containing

[TABLE]

free parameters, where $|\mathbb{S}|$ denotes the size of the set $\mathbb{S}$. The parameterized model is entirely defined by the set $\mathbb{M}=\{\mathbb{S}_{j}\}$, and the model selection consists of choosing this set.

For the related problem of choosing the frequencies in the model for a probability trajectory (see above) there is a clear solution: our instability detection methods are directly testing hypotheses about the frequencies in $p(t)$, and so there is a clear statistical justification for choosing the frequencies in our model for $p(t)$ using the results of these tests. In contrast, the parameters in the model $\mathcal{M}$ do not necessarily have a simple relationship to the circuit probabilities — for an arbitrary model the relationship is entirely arbitrary. So, the relationship between the frequencies detected in the probability trajectories and the frequencies in each $\gamma_{j}$ is also, for a general model, entirely arbitrary. However, for typical models associated with tomographic methods and circuits, under some approximation there is a simple relationship between the model parameters and the predicted probabilities.

To illustrate this, we turn to the Ramsey model with a time-dependent phase that was introduced in the main text:

[TABLE]

If we expand around $lt_{w}\Omega(t)=0$ we obtain

[TABLE]

So, if $lt_{w}\Omega(t)$ is small then the dominant frequencies in $p_{l}(t)$ are also necessarily frequencies in $\Omega(t)$. Equivalent expansions around the zero-error parameter values hold for the general process matrix model of GST — and most other tomographic techniques use a model that is a restriction of that in GST. Therefore, we can use the results of the hypothesis tests on the data from sufficiently short circuits to select the frequencies for the model parameters. But longer circuits are more sensitive to small variations in the model parameters, so we do not want to entirely discount that data in the model selection. To address this problem we can pick a small set of candidate models and then choose between them using the Akaike information criterion (AIC) Akaike (1974). For a parameterized model $\mathcal{M}$ with $k$ free parameters, the AIC score for that model is

[TABLE]

where $\mathcal{L}_{\max}$ is the maximum of the likelihood function for that model. For a set of candidate parameterized models, the preferred model is the one with the minimum AIC score. Moreover, $\mathsf{AIC}(\mathcal{M}_{a})-\mathsf{AIC}(\mathcal{M}_{b})$ is twice the log relative likelihood of model $\mathcal{M}_{b}$ with respect to model $\mathcal{M}_{a}$, so this difference quantifies how much one model is preferred over another. This allows us to decide between multiple time-resolved models.

For the Ramsey experiment presented in the main text we select the expansion

[TABLE]

where $\mathbb{W}=\{1,2,3,4,5,6,7,9,11,12,13,14,15\}$ is the set of all the frequency indices found to be significant in any of the power spectra. We then check that the data justifies including the higher frequencies in the model — because the lowest frequencies appear for circuits where $l$ is small enough so that $lt_{w}\Omega(t)\approx 0$, and so they are almost certainly components in the true $\Omega(t)$. In particular, we select the 10 candidate parameterized models defined by the frequency sets $\mathbb{W}_{0}=\{1,2,3,4\}$, $\mathbb{W}_{1}=\mathbb{W}_{0}\cup\{5\}$, $\mathbb{W}_{2}=\mathbb{W}_{1}\cup\{6\}$ through to $\mathbb{W}_{9}=\mathbb{W}$, and compare their AIC scores, $\mathsf{AIC}_{k}$ for $k=1,2,\dots,9$. The values of $\mathsf{AIC}_{k}-\mathsf{AIC}_{9}$ are all well above zero (ranging from around 10 to around 1000). So the model containing all the frequencies is clearly favored according to the AIC. The detuning estimate reported in the main text is the MLE over this parameterized model.

There is no guarantee that $\mathbb{W}$ contains all of the frequencies that are components in the true $\Omega(t)$. We could also select some larger models to consider, but we do not do this because it would be ad hoc: there are $N-14=5986$ remaining frequencies that we could also include, and there is no clear method for choosing the next frequency to add. We could calculate the AIC for every possible model, but this is clearly impractical. In general, within this framework of time-resolved models based on Fourier decompositions there are $2^{Nk-k}$ possible time-resolved models, where $k$ is the number of parameters in the time-independent model. Moreover, this would not be a reasonable use of the AIC. When adding an erroneous parameter to a model sometimes the AIC score will decrease slightly due to random fluctuations in the data slightly favoring this model, and if we consider adding many different parameters this will almost certainly occur in at least one case.

In the main text we also presented the results of time-resolved Ramsey tomography after we had discarded 80% of our data (i.e., we keep every $5^{\rm th}$ click for each circuit). This analysis followed the same procedure as the analysis on all the data. In particular, note that we implemented independent model selection, i.e., we did not “cheat” by using the time-resolved model selection that was performed using all the data. In particular, using only 80% of the data resulted in a smaller model whereby $\mathbb{W}=\{1,2,3,5,6,7,9\}$. We therefore would not expect the two estimates of $\Omega(t)$ to be equal up to statistical uncertainty (which they are not — see Fig. 1E), as both these statistical uncertainties are standard in-model uncertainties and the two estimates are fits to different models.

B.2 Time-resolved RB

RB is a set of protocols for benchmarking a set of gates on $n$ qubits Knill et al. (2008); Magesan et al. (2011); Proctor et al. (2019). The core protocols all follow the same general procedure:

• (i)

For a range of lengths $m$, run $K$ circuits sampled from some $\Phi_{m}$ distribution whereby (a) the average length of a sampled circuit scales linearly with $m$, and (b) each circuit will deterministically output some $n$-bit string if implemented perfectly.

• (ii)

Process the data to obtain the average success probability ($P_{m}$) at each length $m$, meaning the probability that the expected bit-string was successfully output.

• (iii)

Fit $P_{m}$ to $P_{m}=A+B\lambda^{m}$, and convert $\lambda$ to an error rate via $r=\mathcal{N}(1-\lambda)$ where $\mathcal{N}$ is a constant that we set to $\mathcal{N}=(4^{n}-1)/4^{n}$ where $n$ is the number of qubits.

Within this framework there are a range of different RB protocols, e.g., those of Refs. Knill et al. (2008); Magesan et al. (2011); Proctor et al. (2019), that vary in how the random circuits are chosen. Our methods can be used to upgrade any of these RB protocols to time-resolved techniques.

We focus on the non-intrusive approach to time-resolved RB, which is the following method:

Obtain data by rastering through the sampled circuits. 2. 2.

Use the general stability analysis to estimate the time-resolved success probability $p(t)$ for each circuit. 3. 3.

Average the time-resolved success probabilities for all circuits of the same RB length $m$ to obtain time-resolved average success probabilities $P_{m}(t)$. 4. 4.

Apply the standard RB curve fitting analysis [step (iii) above] to $P_{m}(t_{i})$ at every time of interest $t_{i}$.

The general instability analysis step in this process has some useful flexibility. It has two aspects that can be optimized to a particular task: the exact statistical hypothesis testing on the power spectra, and the selection of the frequencies to include in the time-resolved $p(t)$ models.

Because the RB circuits are random, they are broadly all sensitive to the same physical parameters, with the exact sensitivities depending on the precise details of the random circuit. So the frequencies in $p(t)$ will typically be the same in all the sampled circuits, even if the power in each frequency will generally vary between circuits. In this circumstance the averaged power spectrum is more powerful for detecting these frequencies, because the averaging suppresses measurement shot noise and amplifies signal if the same frequencies appear in most of the spectra. Moreover, here there is little penalty for erroneously including a frequency in the success probability trajectory in a few circuits — we are implementing a global analysis of the data, rather than trying to extract detailed diagnostic information. So we only perform hypothesis testing on the average power spectrum (i.e., we set $w=1$ in Eqs. (24–25)), and then we include every statistically significant frequency in the model for each success probability.

B.3 Time-resolved RB simulation

In the main text we demonstrate time-resolved RB on simulated data. This demonstration is of time-resolved direct randomized benchmarking (DRB) Proctor et al. (2019) on two qubits, with a gate set consisting of a cphase gate between the qubits (i.e., a controlled $\sigma_{z}$), rotations around the $\hat{x}$ and $\hat{y}$ axes by $\pm\frac{\pi}{2}$ and $\pi$ on each qubit, and the Hadamard and idle gates on each qubit. The details of how DRB circuits are sampled (i.e., the $\Phi_{m}$ of DRB) can be found in Ref. Proctor et al. (2019) — they are omitted here as they are not relevant to assessing the performance of time-resolved RB.

We simulated gates with time-dependent coherent phase errors and time-independent depolarizing errors. Each gate is modeled as the perfect gate followed by the error map

[TABLE]

on each qubit, where

[TABLE]

is a one-qubit depolarizing map, and

[TABLE]

We set $\gamma=1-0.04/3$, making the probability of Pauli error 1%. The rotation angle is set to

[TABLE]

where $a=2\pi\times 10^{-5}$, $b=1.5\times 10^{-2}$, $\phi=2\pi\times 10^{-2}$, and $c_{t}$ is a Brownian stochastic process given by $c_{0}=0$ and $c_{t+1}=c_{t}+\mathcal{N}(0,\nu)$ where $\nu=\nicefrac{{3}}{{200}}$. The parameters in this error model have been chosen to give physically plausible phase trajectories, and to not be artificially easy for our methods (e.g., they do not induce $p(t)$ that are exceptionally sparse in the Type-II DCT basis). They have no further significance. The time $t$ starts at 0 and increases by 1 after a single run of any circuit.

In Fig. 2A the estimated time-resolved error rate is compared to the true time-dependent $r(t)$. The value of $r(t)$ is obtained by fitting the exact instantaneous average success probabilities for the set of 100 sampled circuits. This $r(t)$ is the target that time-resolved DRB with this particular set of random circuits is aiming for — it does not correspond to the large-sample-limit error rate that finite-sampling DRB is estimating. We choose this comparison because analysis of time-series data from one particular set of $K$ randomly sampled circuits cannot correct for any finite sampling error introduced when sampling this set of circuits. By defining $r(t)$ in this way, we avoid conflating this finite-circuit-sampling error with any inaccuracies that are intrinsic to time-resolved benchmarking.

B.4 Gate set tomography

GST is a method for jointly reconstructing a set of gates, a state preparation and a measurement Blume-Kohout et al. (2017); Merkel et al. (2013). GST is predicated on the standard Markovian model for errors in QIPs, whereby the prepared state is a fixed density operator, the gates are fixed completely positive and trace preserving (CPTP) maps — which are the linear transformations that map density operators to density operators — and the measurement is a fixed positive-operator valued measure. For our implementation of GST, we parameterize a CPTP gate by $\mathsf{G}=\Lambda\circ\mathsf{G}_{\mathsf{target}}$ where $\mathsf{G}_{\mathsf{target}}$ is the fixed target map and $\Lambda$ is a CPTP map encompassing the errors in the gate. This error map is then parameterized as $\Lambda=\exp(G)$ with $G$ defined by

[TABLE]

where the summations are over $x$, $y$ and $z$, and

[TABLE]

Here, $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{z}$ are the Pauli operators. The $\{H_{k}\}$ maps generate all unitary errors. If $(\theta_{x},\theta_{y},\theta_{z})=\phi\hat{v}$ where $\hat{v}$ is a unit vector then the generated unitary is a rotation around the $\hat{v}$ axis of the Bloch sphere by an angle $\phi$. So $\theta_{x}$, $\theta_{y}$ and $\theta_{z}$ are the components of the unitary error along the $\hat{x}$, $\hat{y}$ and $\hat{z}$ axes in the standard sense. The $\{S_{ij}\}$ maps generate all non-unitary Markovian errors, so the matrix $s$ dictates the stochastic error rates and the size of any amplitude damping errors. In this parameterization, the CPTP constraint is that $s$ is positive semi-definite. Finally, note that this parameterization of $\Lambda$ includes only infinitesimally-generated CPTP maps; it does not strictly permit every possible CPTP map.

B.5 Time-resolved GST simulation

In the main text, time-resolved GST is demonstrated on simulated data. Here we detail the error model that was simulated, and use this example to explain how we apply our methods to implement time-resolved GST. The simulation is for GST circuits out to $L=128$ containing the gates $\mathsf{G}_{\mathsf{i}}$, $\mathsf{G}_{\mathsf{x}}$ and $\mathsf{G}_{\mathsf{y}}$ with an error model that, when expressed in the parameterization above, consists of:

• (i)

Setting the $s$ matrix for all three gates so that if there are no coherent errors $\Lambda$ is a uniform depolarizing channel with a 1% error rate (this is $s\propto\mathds{1}$).

• (ii)

Setting all the coherent error parameters to zero, except

[TABLE]

where $\theta_{k}(\mathsf{G})$ is the $H_{k}$ error rate for gate $\mathsf{G}$, and $f_{k}(t)$ is the $k^{\rm th}$ basis function of the Type-II DCT.

Time starts at 0 and increments by 1 after each raster. The trajectories of the time-dependent parameters over the 1000 simulated rasters are shown in Fig. 2D (and are also included in Fig 6C).

To implement time-resolved GST on this simulated data we first apply the general stability analysis. The resulting

[TABLE]

statistic for every circuit is shown in Fig. 6A, where $\mathsf{p}$ is the p-value of the largest peak in the power spectrum for that circuit. (Due to numerical accuracy limitations, the minimal p-value is $\mathsf{p}=10^{-16}$). The value of $\lambda_{\mathsf{p}}$ varies strongly between circuits of the same approximate length $L$. The strong variations in $\lambda_{\mathsf{p}}$ imply that the time-varying errors are coherent, as stochastic errors have a fairly uniform impact across the GST circuits. Moreover, by inspecting the values of $\lambda_{\mathsf{p}}$ for the simple germs $\mathsf{G}_{\mathsf{i}}$, $\mathsf{G}_{\mathsf{x}}$ and $\mathsf{G}_{\mathsf{y}}$ it becomes clear that $\theta_{z}(\mathsf{G}_{\mathsf{i}})$, $\theta_{x}(\mathsf{G}_{\mathsf{x}})$ and $\theta_{y}(\mathsf{G}_{\mathsf{y}})$ are varying. For example, with the $\mathsf{G}_{\mathsf{i}}$ germ it is the generalized Ramsey circuits for which $\lambda_{\mathsf{p}}$ is statistically significant, as highlighted in Fig. 6B. These circuits amplify $\theta_{z}(\mathsf{G}_{\mathsf{i}})$, so this parameter is likely varying over time. Similarly, it is generalized Rabi sequences that exhibit statistically significant $\lambda_{\mathsf{p}}$ for the $\mathsf{G}_{\mathsf{x}}$ and $\mathsf{G}_{\mathsf{y}}$ germs, and these circuits amplify over/under-rotation errors in these gates, i.e., the parameters $\theta_{x}(\mathsf{G}_{\mathsf{x}})$ and $\theta_{y}(\mathsf{G}_{\mathsf{y}})$. So we conclude that $\theta_{z}(\mathsf{G}_{\mathsf{i}})$, $\theta_{x}(\mathsf{G}_{\mathsf{x}})$ and $\theta_{y}(\mathsf{G}_{\mathsf{y}})$ should be expanded into a summation of Fourier coefficients. Having decided on the set of parameters that are to be expanded into a time-dependent functions, we now use the short-hand $\theta_{\mathsf{i}}\equiv\theta_{z}(\mathsf{G}_{\mathsf{i}})$, $\theta_{\mathsf{x}}\equiv\theta_{x}(\mathsf{G}_{\mathsf{x}})$ and $\theta_{\mathsf{y}}\equiv\theta_{y}(\mathsf{G}_{\mathsf{y}})$.

Deciding on the frequencies to include in the expansion of each of these parameters is more involved. One simple option is to use the frequencies that are statistically significant in the global power spectrum — the set $\mathbb{W}=\{1,2,3,4,97,98,99,100,101,102\}$ — as the frequency set for all three parameters. This parameterized model $(\mathcal{M}_{1})$ includes all of the frequencies in the true model, but also sums and differences of these frequencies, as well as high-frequencies in the expansion of $\mathsf{G}_{\mathsf{i}}$ that are not present in the true model. Yet, implementing MLE over this parameterized model still provides an accurate time-resolved estimate of $\theta_{\mathsf{i}}$, $\theta_{\mathsf{x}}$ and $\theta_{\mathsf{y}}$, as shown in Fig. 6C.

A more sophisticated model selection strategy is to select the frequencies to include in the expansion for $\theta_{\mathsf{i}}$ (the set $\mathbb{W}_{\mathsf{i}}$) to be those frequencies that are statistically significant in the circuits in which $\mathsf{G}_{\mathsf{i}}$ is a germ, and similarly for $\theta_{\mathsf{x}}$ and $\theta_{\mathsf{y}}$. This results in the sets $\mathbb{W}_{\mathsf{i}}=\{1,2,3\}$, $\mathbb{W}_{\mathsf{x}}=\{1,2,3,100,103\}$ and $\mathbb{W}_{\mathsf{y}}=\{1,2,3,5,97,100,101,103\}$, and constitutes a smaller parameterized model ($\mathcal{M}_{2}$). Even this smaller model still contains frequencies that are not in the true model, but each frequency set now only contains sums and differences of the frequencies that appear in the corresponding parameter. The MLE over this parameterized model is what is reported in the main text. The difference in the AIC scores of these two candidate parameterized models is $\mathsf{AIC}(\mathcal{M}_{1})-\mathsf{AIC}(\mathcal{M}_{2})\gg 1$. So the AIC favors discarding the erroneous parameters in $\mathcal{M}_{1}$.

B.6 Time-resolved GST on experimental data

We now detail the exact approach used to implement time-resolved GST on the data from the first of our two GST experiments (as discussed in the main text, we implemented standard time-independent GST on the data from the second experiment). As discussed in the main text, we use the $\lambda_{\mathsf{p}}$ statistics to motivate only expanding the $\hat{z}$-axis rotation angle of the $\mathsf{G}_{\mathsf{i}}$ gate — i.e., using the same notation as above, the parameter $\theta_{z}\equiv\theta_{z}(\mathsf{G}_{\mathsf{i}})$ — into a sum of Fourier components. Specifically, we expand $\theta_{z}$ into

[TABLE]

where $\mathbb{W}=\{1,2,5,7,15\}$ is the set of significant frequencies in the average power spectrum (shown in Fig. 3D). This method for the frequency selection is motivated by noting that, as we have decided that only $\theta_{z}$ exhibits time-dependence, all the frequencies in the averaged power spectrum must be attributable to variations in $\theta_{z}$. So far this is exactly the same type of analysis as with the time-resolved GST on simulated data. However, fitting the model requires more care with experimental data.

First we implemented standard time-independent GST on all the data. The predictions of this GST estimate are compared to the aggregated counts data in Fig. 7. This summarizes the discrepancy between the data and the predictions of the model estimated by standard GST using (A) the log-likelihood ratio (LLR) and (B) the total variational distance (TVD). The LLR test statistic for a circuit is

[TABLE]

Here $\mathcal{L}_{\mathsf{gst}}$ is the likelihood, given the data, of the outcome probabilities for that circuit that are predicted by the GST estimate. Conversely, $\mathcal{L}_{\mathsf{freq}}$ is the likelihood of the outcome probabilities that are equal to the observed frequencies (i.e., for each outcome, set the corresponding probability equal to the number of times that particular outcome was observed divided by the total number of outcome counts for the circuit). Under the null hypothesis that the model of standard GST is true, Wilks’ theorem tells us that the LLR is asymptotically $\chi^{2}_{1}$ distributed. This is because, since we have many circuits, the degrees of freedom per circuit is approximately zero for the GST model (used to compute $\mathcal{L}_{\mathsf{gst}}$) and one for the model maximized over to compute $\mathcal{L}_{\mathsf{freq}}$.The TVD for a circuit is

[TABLE]

where $p_{\mathsf{model}}$ is the probability for the circuit to output 1 predicted by the estimated model and $\mathbf{x}$ is the clickstream for the circuit. We do not use the TVD as a test static (it has less convenient properties for hypothesis testing than the LLR), but the TVD is useful as — unlike the LLR — it quantifies the distance between the observed frequencies and the predicted probabilities. There are statistically significant LLRs in Fig. 7A, so the data are inconsistent with the model of standard GST.

From the general instability analysis, we know that there is some instability in the qubit (see Fig. 3A). So we did not expect the time-independent model of standard GST to be consistent with the data. But, if drift was the only effect outside of GST’s model that is needed to explain the data, we would expect the predictions of the estimated GST model to be inconsistent with the data only for those circuits that exhibit drift. This is not the case; comparing Fig. 7 to Fig. 3A shows that the inconsistency between the estimate of standard GST and the data cannot all be attributed to detectable drift. In particular, the data from many of the long circuits containing germs other than $\mathsf{G}_{\mathsf{i}}$ are irreconcilable with the model of GST, but there is no evidence for drift in these circuits. So there is some other non-Markovian effect that our time-resolved parameterized model does not take into account.

This sort of “model violation” interferes with the model selection and data fitting techniques that we rely on. If we fit our time-resolved parameterized model to all the data, the data that are inconsistent with any set of parameter values can pollute the estimates and render them meaningless. So we first estimate most of the time-independent parameters in the model — specifically, the state preparation, measurement, $\mathsf{G}_{\mathsf{x}}$, $\mathsf{G}_{\mathsf{y}}$ and the non-unitary errors matrix $s$ for $\mathsf{G}_{\mathsf{i}}$ — by implementing ordinary GST on the data from circuits with $L\leq 64$ (these circuits are short enough for the model of GST to be reasonably consistent with the data). Only the two time-independent coherent error parameters in $\mathsf{G}_{\mathsf{i}}$ and $\theta_{z}$ — which is a summation of 6 Fourier coefficients — are not fixed by this initial fit. Using MLE, we then fit this 8-parameter time-resolved model to the time-series data for all circuits with $L\leq 64$ and all the circuits where $\mathsf{G}_{\mathsf{i}}$ is a germ with $L\leq 1024$. This allows us to extract information about the time-dependent parameter without polluting the results with data that does not fit this model.

B.7 Estimated gate error rates

In the main text, we report the diamond distance error rate ($\epsilon_{\diamond}$) of the estimated (time-resolved) gate process matrices in the two GST experiments. Table 1 includes additional information on the gate error rates. In addition to $\epsilon_{\diamond}$, this table reports the entanglement infidelity ($\epsilon_{\mathsf{f}}$) and spectral entanglement infidelity ($\epsilon_{\mathsf{spec},\mathsf{f}}$) for each gate. The diamond distance error rate for process matrix $\mathsf{G}$ is

[TABLE]

where $\|\cdot\|_{\diamond}$ is the diamond norm Watrous (2005) and $\mathsf{G}_{\text{target}}$ is the target process matrix that $\mathsf{G}$ ideally should be. The entanglement infidelity is

[TABLE]

where $\varphi$ is any maximally entangled state Nielsen (2002). Although the entanglement infidelity is defined with respect to a maximally entangled state, note that it is well-defined for one-qubit gates and it does not imply that any entangling operations are being performed. Another commonly used gate infidelity is the average gate infidelity ($\bar{\epsilon}$). The entanglement infidelity is related to the average gate infidelity via a linear rescaling. Specifically, $\epsilon_{\textsf{f}}=(d+1)\bar{\epsilon}/2$ Nielsen (2002) where $d$ is the dimension, so here $d=2$. We use the entanglement infidelity because $\epsilon_{\diamond}=\epsilon_{\textsf{f}}$ for depolarizing errors, so the difference between $\epsilon_{\diamond}$ and $\epsilon_{\textsf{f}}$ is a rough measure of how structured the errors are, e.g., for coherent errors $\epsilon_{\diamond}>\epsilon_{\textsf{f}}$.

The spectral entanglement infidelity ($\epsilon_{\mathsf{spec},\mathsf{f}}$) is defined to be

[TABLE]

where $\mathsf{Spec}(\mathsf{A})$ is the diagonal matrix that has the eigenvalues of the matrix $\mathsf{A}$ on the diagonal (and the relative ordering of the eigenvalues in the two matrices is chosen to minimize $\epsilon_{\mathsf{spec},\mathsf{f}}$). The spectral entanglement infidelity is an alternative to the standard entanglement infidelity that is invariant under transformation of the “GST gauge” (see Ref. Blume-Kohout et al. (2017) for discussions of this gauge).

Looking at Table 1 we see that all three gates have a smaller diamond distance error rate in the second experiment, with a decrease from $\epsilon_{\diamond}\sim 0.5\%$ to $\epsilon_{\diamond}\sim 0.05\%$. The infidelity also decreases for $\mathsf{G}_{\mathsf{y}}$, but, in contrast, for $\mathsf{G}_{\mathsf{i}}$ and $\mathsf{G}_{\mathsf{x}}$ the infidelity increases from $\epsilon_{\textsf{f}}\sim 0.003\%$ to $\epsilon_{\textsf{f}}\sim 0.006\%$. The diamond distance error rate is sensitive at order $\theta$ to a coherent error with a rotation angle of $\theta$, whereas the infidelity is sensitive to coherent errors only at order $\theta^{2}$. So, a decrease in $\epsilon_{\diamond}$ and an increase in $\epsilon_{\textsf{f}}$ is a sign that the magnitude of the coherent errors (i.e., calibration errors) has been reduced in the second experiment, but that the rate of stochastic errors has slightly increased.

Appendix C Experiment Design

In this section we discuss how to design experiments to efficiently detect and characterize instability. This includes details on the relationship between the amount of data taken and the sensitivity of our methods.

C.1 How much data do we need?

Consider looking for instability in a given set of circuits. For example, we could be running an algorithm, or syndrome extraction circuits, or the circuits from some benchmarking routine. Furthermore, assume that we are going to collect data by rastering through this set of circuits $N$ times. It then only remains to choose $N$, the length of the clickstream for each circuit. In this circumstance, it is important to understand how the sensitivity of our methods varies with $N$.

Each circuit has some discrete-time probability trajectory $\mathbf{p}=(p_{0},p_{1},\dots,p_{N-1})$. We can always express this as a summation over the basis vectors of $F$. As we use the Type-II DCT for all of our explicit data analysis, we will use this basis for clarity. In this basis, we write

[TABLE]

for some $\{\gamma_{\omega}\}$ amplitudes. We can use our statistical model for the frequency-domain data (see earlier) to approximately calculate the probability of successfully detecting that some $\gamma_{\omega}$ is non-zero, as a function of $\gamma_{\omega}$, $N$ and the test significance $\alpha$. In particular, the probability of detecting that $\mathbf{p}$ has a nonzero component at $\omega$, from a clickstream of length $N$ and with tests at a statistical significance of $\alpha$, is approximately

[TABLE]

Here, $\text{erf}(\cdot)$ is the error function,

[TABLE]

and $T(\alpha,N)$ is the power significance threshold for the hypothesis test at that frequency. In the case of testing a single clickstream, and splitting the test significance equally over all $N-1$ non-zero frequencies (as we do for all our data analysis, but which is not necessarily always desirable), this is given by

[TABLE]

where $\mathsf{CDF}_{k}$ is the $\chi^{2}_{k}$ cumulative distribution function. If testing $C$ clickstreams, we let $\alpha\to\alpha/C$ in this equation.

Fig. 8 shows a heatmap of $P_{\text{det}}(N,\gamma,0.05)$ as $N$ and $\gamma$ are varied. From this plot we can read off the clickstream length required to reliably detect an instability of a given amplitude. It shows that only moderate amounts of data are required to detect relatively small instabilities, e.g., an instability of amplitude $\gamma\approx 0.1$ can be reliably detected with $N\approx 1000$. This quantity of data is routinely taken to estimate static circuit outcome distributions.

As well as detecting instabilities in probability trajectories, our methods also estimate $p(t)$. To quantify the accuracy of our probability trajectory estimation we use simulations. We simulated clickstream data sampled from a pure-tone probability trajectory $\mathbf{p}(\gamma)$:

[TABLE]

Fig. 9 shows both the instability detection probability and the error in our estimator. Fig. 9B shows the root mean squared error (RMSE) of our estimator for three values of $\gamma$ (unbroken lines), alongside the corresponding baselines set by estimating $p(t)$ as a constant given by the mean of the clickstream (dotted lines). This is the mean RMSE, obtained in $10^{5}$ simulations, and note that this simulation used the “Fourier filter” variant of our estimator (as it is faster to calculate than the MLE version). In the regime where $N$ is so low that instability is almost never detected, the RMSE of our estimator is essentially equal to the RMSE estimator of the baseline clickstream-mean estimate. This is because our estimate is the mean of the clickstream unless instability is detected. In fact, when $N$ is so low that instability is almost never detected (as shown in Fig. 9A) our estimator actually performs worse than the baseline, on average. This is because in approximately $5\%$ of instances (for testing at 5% statistical significance) our detection routine will spuriously detect that there is a component of $\mathbf{p}$ at some random frequency, and so we will spuriously add a component at that frequency into our estimator. Comparing Fig. 9A and 9B, we see that the RMSE of our estimator rapidly improves once $N$ moves into the regime where instability at frequency $\omega$ is successfully detected with high probability. Once the detection probability saturates at essentially 1, the RMSE then follows $O(\nicefrac{{1}}{{\sqrt{N}}}$) scaling. This is because the RMSE in the estimate of each amplitude in our selected model is essentially shot-noise limited.

So far we have considered how long a single clickstream needs to be to reliably detect instabilities of a given magnitude. But our analysis can also jointly test the data from multiple circuits, by calculating and then testing an averaged power spectrum. This allows for much greater sensitivity for a fixed $N$, where $N$ is the clickstream length for each circuit. Fig. 10 demonstrates this enhancement, by plotting the minimum $N$ required for a 50% detection probability ($N_{\min}$) when averaging spectra from 1, 10, 100 and 1000 circuits. This is for the case of clickstreams that are all sampled from the same pure-tone $\mathbf{p}(\gamma)$. This case is chosen as it provides the maximal enhancement from averaging spectra, and — although it is unlikely that multiple circuits will have identical $p(t)$ — there are many circumstances under which many or all of a set of circuits would contain drift at the same frequencies.

Fig. 10 shows that increasing the number of circuits provides a substantial reduction in $N_{\min}$. For example, if $\gamma=0.1$ then $N_{\min}\approx 100$ for 100 clickstreams whereas $N_{\min}\approx 1000$ for a single clickstream. But note that the sensitivity enhancement from increasing the number of circuits is smaller than that from increasing the clickstream length. As illustrated in the above example, to obtain a 50% detection probability at $\gamma=0.1$ we need 10 times more bits in total in the 100 clickstream case than in the single clickstream case. Essentially, this is due to the difference between the $O(\nicefrac{{1}}{{\sqrt{C}}})$ shot-noise suppression obtained when averaging $C$ spectra with independent statistical fluctuations verses the linear growth with $N$ of the power in $\mathbf{p}$ (for $\mathbf{p}$ with an $N$-independent variance).

C.2 Choosing the circuit repetition rate

Our methods can be used for both dedicated drift characterization and for auxiliary data analysis. For dedicated drift characterization we need to carefully chose a set of quantum circuits that are as sensitive as possible to the specific parameters under study, as we discussed in the main text. Greater sensitivity results in large fluctuations in the circuits’ outcome probabilities, for a given instability on the parameters of interest, meaning that fewer data are needed to reliably detect and characterize the drift — as quantified above (see Fig. 8-10). But there is also another important aspect for designing dedicated drift characterization experiments: the rate that the circuits need to be repeated. If the dedicated drift characterization circuits can be run sufficiently infrequently then it becomes possible to interleave them with application circuits. Furthermore, understanding the effect of reducing the repetition rate is important even in experiments run only for the purpose of drift characterization — as we will often have to decide between running more circuits, or repeating fewer circuits at a faster rate. For example, in the case of both GST and Ramsey experiments, we must choose the maximum circuit length, and this directly impacts the rate at which we can raster through the full set of circuits.

Reducing the repetition rate of the drift characterization circuits — i.e., increasing the time between each run of a circuit, $t_{\rm gap}$ — will have two effects. Firstly, for a fixed experiment duration we will obtain fewer data points, i.e., the clickstreams will be shorter. The impact of the clickstream length on the sensitivity of our methods has already been discussed in detail earlier (see Fig. 8-10). But note that, for dedicated drift characterization, we are free to choose the circuits that maximally amplify the parameters we expect might be drifting — so that the oscillations in the $p(t)$ need only be small if they are small in every possible circuit (under the assumption that we have correctly identified all the parameters that might be drifting). Therefore, relatively low $N$ will typically still allow for sensitive drift characterization. In the main text we demonstrated this, by discarding 80% of our Ramsey data (giving $N=1200$) and showing that we still obtain a high-precision time-resolved estimate of the qubit detuning (see Fig. 1E). Moreover, much more of that data could be discarded if we only wanted to detect drift and obtain a rough time-resolved estimate of the detuning.

The second impact of increasing $t_{\rm gap}$ is that the sensitivity to high-frequency drift is reduced. As with all signal analysis techniques that use regularly sampled data, there is a Nyquist limit — i.e., oscillations at a frequency of $\nicefrac{{1}}{{2t_{\rm gap}}}$ are undetectable, and oscillations at higher frequencies are either aliased, or are undetectable if they are integer multiples of this frequency. The repetition rate should be chosen to maintain sensitivity to the frequency range of interest. As illustrated by our analysis of the Ramsey data, this means that we can have a slow repetition rate if there is only low-frequency drift.