# Dynamics of superconducting qubit relaxation times

M. Carroll * * These authors contributed equally<EMAIL_ADDRESS>IBM
Quantum, IBM T.J. Watson Research Center,Yorktown Heights, NY, 10598, USA  S.
Rosenblatt* * These authors contributed equally<EMAIL_ADDRESS>IBM
Quantum, IBM T.J. Watson Research Center,Yorktown Heights, NY, 10598, USA  P.
Jurcevic  IBM Quantum, IBM T.J. Watson Research Center,Yorktown Heights, NY,
10598, USA  I. Lauer  IBM Quantum, IBM T.J. Watson Research Center,Yorktown
Heights, NY, 10598, USA  A. Kandala<EMAIL_ADDRESS>IBM Quantum, IBM T.J.
Watson Research Center,Yorktown Heights, NY, 10598, USA

###### Abstract

Superconducting qubits are a leading candidate for quantum computing but
display temporal fluctuations in their energy relaxation times $T_{1}$. This
introduces instabilities in multi-qubit device performance. Furthermore,
autocorrelation in these time fluctuations introduces challenges for obtaining
representative measures of $T_{1}$ for process optimization and device
screening. These $T_{1}$ fluctuations are often attributed to time varying
coupling of the qubit to defects, putative two level systems (TLSs). In this
work, we develop a technique to probe the spectral and temporal dynamics of
$T_{1}$ in single junction transmons by repeated $T_{1}$ measurements in the
frequency vicinity of the bare qubit transition, via the AC-Stark effect.
Across 10 qubits, we observe strong correlations between the mean $T_{1}$
averaged over approximately nine months and a snapshot of an equally weighted
$T_{1}$ average over the Stark shifted frequency range. These observations are
suggestive of an ergodic-like spectral diffusion of TLSs dominating $T_{1}$,
and offer a promising path to more rapid $T_{1}$ characterization for device
screening and process optimization.

##  I Introduction

Superconducting qubits are a leading platform for quantum computing Zhang  _et
al._ (2020); Arute  _et al._ (2019). This has been driven, in part, by
improvements in coherence times over five orders of magnitude since the
realization of coherent dynamics in a cooper pair box Nakamura, Pashkin, and
Tsai (1999). However, further improving coherence times remains crucial for
enhancing the scope of noisy superconducting quantum processors as well as the
long term challenge of building a fault tolerant quantum computer. Recent
advances Hong  _et al._ (2020); Kandala  _et al._ (2020); Hashim  _et al._
(2020); Foxen  _et al._ (2020) in two-qubit gate control have placed their
fidelities at the cusp of their coherence limit, implying that improvements in
coherence could directly drive gate fidelities past the fault tolerant
threshold. In this context, coherence stability and its impact on multi-qubit
device performance is also an important theme, since superconducting qubits
have been shown to display large and correlated temporal fluctuations (i.e.,
$1/f^{\alpha}$) in their energy relaxation times $T_{1}$Müller  _et al._
(2015); Paladino  _et al._ (2014); Weissman (1988); Klimov  _et al._ (2018);
Kandala  _et al._ (2019); Burnett  _et al._ (2019); Schlör  _et al._ (2019).
This places additional challenges for benchmarking the coherence properties of
these devices Burnett  _et al._ (2019), and also for error mitigation
strategies such as zero noise extrapolation Kandala  _et al._ (2019).

The fluctuations of qubit $T_{1}$ are often attributed to resonant couplings
with two level systems (TLSs) that have been historically studied in the
context of amorphous solids Phillips (1972); Müller, Cole, and Lisenfeld
(2019) and their low temperature properties. More recently, TLSs have
attracted renewed interest due to their effect on the coherence properties of
superconducting quantum circuitsMartinis  _et al._ (2005); Grabovskij  _et
al._ (2012); Barends  _et al._ (2013); Klimov  _et al._ (2018); Lisenfeld  _et
al._ (2019); Burnett  _et al._ (2019); Schlör  _et al._ (2019), and are
attributed to defects in amorphous materials at surfaces, interfaces, and the
Josephson junction tunnel barrier. Frequency resolved measurements of $T_{1}$
in flux and stress tunable devicesBarends  _et al._ (2013); Klimov  _et al._
(2018); Lisenfeld  _et al._ (2019) have also displayed fluctuations,
suggesting an environment of TLSs with varying coupling strengths around the
qubit frequency. The variability of $T_{1}$ over time is explainedMüller,
Cole, and Lisenfeld (2019); Klimov  _et al._ (2018), at least in part, by
temporal fluctuations in this frequency environment, associated with the
spectral diffusion of the TLSs Black and Halperin ; Phillips (1972).

Furthermore, two-qubit gates that involve frequency excursions Kandala  _et
al._ (2020); Stehlik  _et al._ (2021); Klimov  _et al._ (2018) can also
interact with TLS near the qubit frequency leading to additional incoherent
error. The fluctuations in TLS peak positions, therefore, can also introduce
fluctuations in two qubit fidelity. Spectroscopy of defect TLS is, therefore,
central to understanding the short and long time $T_{1}$ and gate fidelity of
qubits.

Figure 1: (a) Map of a 20 qubit device including the qubit number, two-qubit
connectivity and the qubit frequency in (GHz). (b) Box and whisker plot of
$T_{1}$’s measured daily over approximately nine months. The box covers the
first to third quartile, the vertical line indicates the median and the
whiskers are drawn to the maximum or minimum values that fall within 1.5 times
the interquartile range. All other points are outliers.

Single Josephson junction transmons with fixed frequency couplings represent a
successful device architecture achieving networks of over 60 qubits Zhang  _et
al._ (2020) with all microwave control and state of the art device coherence.
The single junction configuration offers advantages such as reduced
sensitivity to flux noise, while preserving the transmon charge insensitivity
and also reducing system complexity with fewer control inputs. However, there
is little TLS spectroscopy of single junction transmons because of the limited
tunability, despite the central importance of understanding the TLS
environment both for device and process characterization.

Figure 2: (a) Schematic of qubit $T_{1}$ response to shift in qubit frequency
(red traces), at different times, $t$. The qubit frequency is tuned by an
amount $\Delta\omega_{q}$ by an off-resonant tone placed $|\Delta_{qs}|$ above
the qubit frequency $\omega_{q}$. The dependence of the frequency shift on the
off-resonant tone amplitude $\Omega_{s}$ is depicted by the blue trace. The
$T_{1}$ dips are indicative of the qubit coming into resonance with a strongly
coupled TLS, at the frequency $\omega_{TLS}$. The TLS frequency shifts in time
due to spectral diffusion, schematically indicated by changes in the $T_{1}$
dips at different time snapshots. (b) Schematic of a Ramsey pulse sequence
used to calibrate $\Delta\omega_{q}$ as a function of $\Omega_{s}$ (top); and
(bottom) schematic pulse sequence for the relaxation time spectroscopy. For
each $\Omega_{s}$ (i.e., $\Delta\omega_{q}$), the $|1\rangle$ occupation is
measured at a fixed time (i.e., $\tau$ = 50 $\mu$s in this work). (c) An
illustrative case of the 01 transition dependence on $\omega_{s}$ for constant
$\Omega_{s}$. The 01 qubit frequency, $\omega_{q}$ \+ $\Delta\omega_{q}$, uses
an unperturbed frequency of $\omega_{q}$ = 5.0 GHz and an anharmonicity of
$\delta_{q}$ = -340 MHz. The locations of the unperturbed 01 and 12
transitions are shown as vertical lines overlaid with 5 MHz offsets to make
their locations more visible on the figure. Negative and positive qubit shifts
can be produced and large shifts can be induced depending on $\Delta_{qs}$.
(d) Measured $\Delta\omega_{q}$ as a function of normalized DAC amplitude,
$\Omega_{s}$ using the AC Stark shifted Ramsey technique. Solid line is a
quadratic fit functionally consistent with a perturbative model.

In this work, we introduce an all-microwave technique for the fast
spectroscopy of TLSs in single junction transmon qubits that requires no
additional hardware resources. In contrast to flux based approaches to TLS
spectroscopy, we employ off-resonant microwave tones to drive AC-Stark shifts
of the fundamental qubit transition and spectrally resolve qubit relaxation
times. Dips in relaxation times serve as a probe of the frequency location of
a strongly coupled TLS. We use repeated frequency sweeps to probe the time
dynamics of the relaxation probabilities including tracking the spectral
diffusion of strongly coupled TLS. Across 10 qubits, we observe strong
correlations between the long time mean, averaged over several months
$\langle{T_{1}}\rangle_{T}$, and the short time mean, averaged around the
local qubit frequency $\langle{T_{1}}\rangle_{\omega,t}$.

This strong correlation suggests a quasi-ergodic behavior of the TLS spectral
diffusion in the nearby frequency neighborhood of the qubit. In contrast,
there is lower correlation between $\langle{T_{1}}\rangle_{T}$ and $T_{1}$
measured over a single day. The $\langle{T_{1}}\rangle_{\omega,t}$ can
provide, therefore, a more rapid estimate of long time behavior.

##  II Device and spectroscopy technique

The experiments reported in this letter were performed on ibmq_almaden, a 20
qubit processor based off single junction transmons and fixed couplings. The
device topology is shown in Fig. 1 (a), and qubit frequencies are around
$\sim$5 GHz. Fig. 1 (b) depicts the characteristic spread of the qubit
$T_{1}$s and their mean, from $\sim$250 measurements over 9 months. The base
plate (to which the device was mounted) temperature of the dilution
refrigerator was typically $\sim$13 mK excepting several temperature
excursions to $\sim$1 K, which were not observed to have any significant
effects on the  long time $T_{1}$ time series or distributions of $T_{1}$
values discussed in this work, discussed later. Several qubits on the device
display mean $T_{1}$s exceeding 100 $\mu$s. However, the large spread in
individual qubit $T_{1}$s highlights the challenge for rapid benchmarking of
device coherence, since any single $T_{1}$ measurement can disagree
substantially from its longtime mean.

We study the spectral dynamics of these $T_{1}$ times by employing off-
resonant microwave tones Gambetta  _et al._ (2006) to induce an effective
frequency shift $\Delta\omega_{q}$ in single junction transmons by the AC
Stark effect. This has been employed previously for coherent state transfer
between coupled qubits that are Stark shifted into resonance Majer  _et al._
(2007). In this work, shifting the qubit frequency into resonance with a
defect TLS induces a faster relaxation time, which in turn is used to detect
the frequency location of the TLS Simmonds  _et al._ (2004), as depicted in
Fig. 2 (a). The Stark shift can be described analytically by a Duffing
oscillator model Magesan and Gambetta (2020); Schneider  _et al._ (2018)

$\Delta\omega_{q}=\frac{\delta_{q}\Omega_{s}^{2}}{2\Delta_{qs}(\delta_{q}+\Delta_{qs})}$
(1)

where $\delta_{q}$ is the qubit anharmonicity, $\Omega_{s}$ is the drive
amplitude and $\Delta_{qs}=\omega_{q}-\omega_{s}$ is the detuning between the
qubit frequency and the Stark tone.

As seen from the expression above, the magnitude and sign of the Stark shift
can be manipulated by the detuning and the drive amplitude of the Stark tone,
Fig. 2 (c). Very large frequency shifts can be obtained by driving close to
the transmon transitions, but this typically leads to undesired
excitations/leakage out the two-state manifold. In this work, we obtain Stark
shifts of 10’s of MHz, with modest drive amplitudes and a fixed detuning
$\Delta_{qs}$ of $\pm$50 MHz. The frequency shifts are experimentally measured
using a modified Ramsey sequence Ramsey (1950), schematically shown in Fig. 2
(b), and display good agreement with the quadratic dependence of the
perturbative model in the low-drive limit. A representative case is shown in
Fig. 2 (d).

Figure 3:  (a) Measured probability of being in the $|1\rangle$ state,
$P_{1}$, at 50 $\mu$s wait time with varying $\Delta\omega_{q}$ and tone
detuned 50 MHz above ${\omega_{q}/2\pi}$ for $Q_{19}$. (b) Example of Ramsey
measurements used to extract frequency shifts ${\Delta\omega_{q}/2\pi}$ from
pulse amplitudes. The two curves result from starting the Ramsey oscillations
with a $X_{\pi/2}$ or $Y_{\pi/2}$. (c) and (d) are $T_{1}$ measurement with
Stark shifts $\Delta{\omega_{q}/2\pi}$ = -18.9 and -5.3 MHz, respectively. (e)
$T_{1}$ measurement with no Stark shift (i.e., no Stark tone).

We focus on the spectrally resolved $T_{1}$ measurements in Fig. 3 that we use
as a probe of defect TLS transition frequencies. However, instead of measuring
the entire $T_{1}$ decay, we use the excited state probability, $P_{1}$, after
a fixed delay time as a measure of $T_{1}$. This speeds up the spectral scans
significantly. Our experiments are performed at a repetition rate of 1 kHz,
but our scheme can be further accelerated with reset techniques Egger  _et
al._ (2018), which can be crucial for probing faster TLS dynamics. For an
effective frequency sweep, we run an amplitude sweep with off-resonant pulses
at fixed detuning ($\pm$50 MHz) and duration (delay time of 50 $\mu$s), after
exciting the qubit with an initial $\pi$ pulse. The pulsed Stark sequence
enables faster spectroscopy by circumventing the need to re-calibrate the
$\pi$, $\pi/2$ pulses at every frequency. The off-resonant pulses have
Gaussian-square envelopes with a 2$\sigma$ rise-fall profile, where $\sigma$ =
10 ns. This pulse sequence is shown in Fig. 2 (b). The amplitude points in the
sweep are then related to Stark shifts by Ramsey sequences. Fig. 3 shows
representative data of such a sweep on qubit 19 ($Q_{19}$) with distinctive
dips in $P_{1}$ that we attribute to strongly coupled TLS at their transition
frequencies. $T_{1}$ measurements at Stark tone amplitudes corresponding to
high/low $P_{1}$ points, as seen in the bottom panel of Fig. 3, explicitly
show the substantial variation in $T_{1}$ as a function of frequency and the
consistent tracking of $T_{1}$ with $P_{1}$.

Variations in $P_{1}$ can potentially be caused by sources other than TLS. In
our experiments, $P_{1}$ is spectrally resolved to $\sim\pm$25 MHz around the
individual qubit frequencies. The narrow frequency range combined with
measuring non-neighbor sets of qubits simultaneously avoids strong $P_{1}$
suppression from resonances with neighboring qubits, the coupling bus or
common low-Q parasitic microwave modes. Control experiments show that time
insensitive features in the $P_{1}$ fingerprint are robust to choice of the
Stark tone detuning, ruling out a power dependence for the power range used in
this work. Finally, while a recent report Yan  _et al._ (2016) modeled their
broadband $T_{1}$ scatter as arising from quasi-particle fluctuations, this is
not sufficient to explain the sharp frequency dependent $P_{1}$ features
depicted, for instance in Fig. 3. Furthermore, recent experiments on our
qubits suggest a quasi-particle limit to $T_{1}$ that exceeds several
milliseconds.Kurter  _et al._ (2020)

##  III TLS dynamics and correlations between $P_{1}(\omega,t)$ and $\langle
T_{1}\rangle_{T}$

We repeat the line traces of Fig. 3 for both positive and negative 50 MHz
detuning, approximately once every 3-4 hours, extended over hundreds of hours
for all the qubits. A representative example of the cumulative scans is shown
in Fig. 4 for $Q_{15}$. Spectroscopy of the other qubits is shown in the
supplemental information, appendix A. The TLS dynamics around the qubit
frequency are qualitatively similar to previous TLS spectroscopy using flux or
stress tunable devices Müller, Cole, and Lisenfeld (2019).

In the case of $Q_{15}$, Fig. 4, there are prominent dips in relaxation
probability around positive 1 MHz, negative 5-10 MHz, and negative 15-20 MHz.
The spectral diffusion of the positions of the $T_{1}$ dips can vary between
order of 1 MHz to 10 MHz over the 272 hours of measurement providing a
qualitative measure of linewidths. A more quantitative discussion of
linewidths can be found in appendix I. The background is covered by an
ensemble of smaller dips of relaxation, Fig. 3, that also dynamically evolve,
with features that are larger than the sampling noise in the measurement.

Figure 4: Time dependence of the energy relaxation spectroscopy for $Q_{15}$
using $\Delta\omega_{s}=\pm$50 MHz and varying $\Omega_{s}$ over 501 different
j amplitude points, $\Omega_{s,j}$, to sweep $\Delta\omega_{q}(\Omega_{s,j})$
in a positive or negative direction (see eqn. 1). The $P_{1}$ is measured at
50 $\mu$s.

As discussed previously, $T_{1}$ fluctuations introduce uncertainty in the
coherence benchmarking, stability of multi-qubit circuit performance and
process optimization of superconducting qubit devices. In this context of
better estimator, we examine if the long time averages ($T\sim$ 9 months)
$\langle T_{1}\rangle_{T}$ and $\langle P_{1}\rangle_{T}$ are correlated with
the frequency neighborhood of the qubit $\langle T_{1}\rangle_{\omega,t}$ and
$\langle P_{1}\rangle_{\omega,t}$, respectively. The averaged relaxation
probabilities and $T_{1}$s are defined as

$\langle
P_{1}\rangle_{T}=\frac{1}{N}\sum_{i=1}^{N}P_{1}(\omega_{q},\tau,T_{i})$ (2)

$\langle T_{1}\rangle_{T}=\frac{1}{N}\sum_{i=1}^{N}T_{1}(\omega_{q},T_{i})$
(3)

$\langle
P_{1}\rangle_{\omega,t}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{2\Delta\omega}\sum_{\omega_{j}=-\Delta\omega}^{\Delta\omega}P_{1}(\omega_{q}+\omega_{j},\tau,t_{i})d\omega_{j},$
(4)

$\langle
T_{1}\rangle_{\omega,t}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{2\Delta\omega}\sum_{-\Delta\omega}^{\Delta\omega}\frac{-\tau}{\ln(P_{1}(\omega_{q}+\omega_{j},\tau,t_{i})d\omega_{j})},$
(5)

where definitions of variables can be found in table 1.

We compare $\langle P_{1}\rangle_{\omega,t}$ to $\langle P_{1}\rangle_{T}$
from the daily $T_{1}$ measurements over $T_{\mathrm{max}}\sim$ 9 months
evaluated at $\tau=53~{}\mu s$, shown in Fig. 1. The $\langle
P_{1}\rangle_{\omega,t}$ are calculated for a $T_{1}$ delay time of $\tau$ =
50 $\mu$s for 10 qubits in the device for the first time slice and a cutoff
frequency $\Delta\omega/2\pi$ = 5 MHz. A qualitatively close agreement for all
10 qubits is observed, see Fig. 5 (a).

Symbol |  Definition
---|---
$\omega_{q}$ | Qubit frequency
$\Omega_{s,j}$ |  $j^{th}$ Stark drive amplitude point
| in a frequency scan (see eqn. 1)
$\omega_{j}$ | $j^{th}$ Qubit Stark shift,
| $\Delta\omega_{q}(\Omega_{s,j})$ (see eqn. 1)
$d\omega_{j}$ | Frequency bin size at the
| Stark shifted frequency location
$\Delta\omega$ | Maximum qubit Stark shift
$N$ | Number of $T_{1}$ measurements
| for $\sim$9 month time series in an average
$n$ | Total number of spectroscopy
| time slices in a moving average
$T_{i}$ | Time of $i^{th}$ time bin for the
|  $\sim$9 month $T_{1}$ time series
$t_{i}$ | Time of $i^{th}$ time bin
| of the spectroscopy time series
$T_{1}(T_{i})$ | $T_{1}$ measured at $i^{th}$ time
$\tau$ | Decay time at
| which $P_{1}$ was evaluated
$P_{1}(\omega_{q}+\omega_{j},\tau,t_{i})$ | Probability of $|1\rangle$ at $\tau$ delay for
| time slice $t_{i}$ and frequency shift $\omega_{j}$
$\langle P_{1}\rangle_{T}$ | Probability of $|1\rangle$ at $\tau$
| averaged over $\sim$9 month time series
$\langle P_{1}\rangle_{\omega,t}$ | Probability of $|1\rangle$ at $\tau$ averaged over
| frequency and spectroscopy time
$\langle T_{1}\rangle_{T_{0\rightarrow N}}$ | Moving average of $T_{i}$
| measurements from $T_{0}$ to $T_{n}$
$\langle T_{1}\rangle_{T}$ |  $T_{1}$ averaged over entire $\sim$9 months
$\langle T_{1}\rangle_{\omega,t}$ |  $T_{1}$ average from $\langle P_{1}\rangle_{\omega,t}$
$Q_{k}$ | k’th qubit
$R(t_{i})$ |  $\langle T_{1}\rangle_{T}$ correlation to $\langle T_{1}\rangle_{\omega}$ at $t_{i}$ for {$Q_{k}$}
$R(T_{i})$ |  $\langle T_{1}\rangle_{T}$ correlation to $T_{1}$ at a single time
|  $T_{i}$ of the $\sim$9 month series for {$Q_{k}$}
$\langle R\rangle_{t_{0\rightarrow n}}$ | $\langle T_{1}\rangle_{T}$ correlation to $\langle T_{1}\rangle_{\omega,t_{0\rightarrow n}}$ for {$Q_{k}$}
$\langle R\rangle_{T_{0\rightarrow N}}$ | $\langle T_{1}\rangle_{T}$ correlation to $\langle T_{1}\rangle_{T_{0\rightarrow N}}$ for {$Q_{k}$}
Table 1: List of symbols

A $\langle T_{1}\rangle_{\omega,t}$ can also be estimated for each $\langle
P_{1}\rangle_{\omega,t}$ at $\tau$ = 50 $\mu$s by assuming an exponential
decay. The approximate equivalence of $\langle T_{1}\rangle_{\omega,t}$ and
$\langle T_{1}\rangle_{T}$ is seen in the scatter plot of Fig 5 (a) inset. A
near 1:1 relationship is observed when this approach is applied more broadly
across many IBM devices, see Appendix K. Furthermore, the poorer correlation
between $\langle T_{1}\rangle_{T}$ and a single instance of $T_{1}$
measurements, is also shown by larger scatter, as seen in Fig 5 (a) inset.

To quantify with a single value the correlation between $\langle
T_{1}\rangle_{T}$ or $\langle P_{1}\rangle_{T}$ and their estimators for many
qubits, we use a Pearson $R$ measure across the ten odd-labeled qubits,

$R=\frac{\sum\limits_{k=0}^{d-1}(\langle X\rangle_{T,Q_{i}}-\overline{\langle
X\rangle_{T}})(\langle X\rangle_{\omega,t,Q_{i}}-\overline{\langle
X\rangle}_{\omega,t})}{\sqrt{\sum\limits_{k=0}^{d-1}(\langle
X\rangle_{T,Q_{k}}-\overline{\langle
X\rangle_{T}})^{2}\sum\limits_{k=0}^{d-1}(\langle
X\rangle_{\omega,t,Q_{i}}-\overline{\langle X\rangle}_{\omega,t})^{2}}}$ (6)

where d is the number of qubits in the device or analysis, 10 in this case,
and X is the observable $P_{1}$ or $T_{1}$. The Pearson correlation is a
normalized covariance between two variables reflecting a linear correlation
from 1 to -1, where $R$ = 1 (-1) represents a $100\%$ positive (negative)
correlation and $R$ = 0 indicates no correlation. Strong R correlation can
therefore signal the existence of a potential linear mapping between the
estimator and $\langle T_{1}\rangle_{T}$, in particular, possibly one that is
1:1 or a scaling factor that will reliably estimate $\langle
T_{1}\rangle_{T}$.

For a single frequency sweep that takes $\sim$20 minutes, we obtain 0.76 $<$
$R(t_{i})$ $<$ 0.84 correlation between $\langle T_{1}\rangle_{T}$ and
$\langle T_{1}\rangle_{\omega,t}$ for 0.5 MHz $<$ $\Delta\omega$ $<$ 5 MHz.
Using the $P_{1}$ values without assuming an exponential dependence leads to
stronger correlations of 0.87 $<$ $R(t_{i})$ $<$ 0.91. Both of these are
substantially stronger than the correlation found between the representative
instance of $T_{1}$ and $\langle T_{1}\rangle_{T}$, which was $R$ = 0.29. We
note this instance of $R$ can have a large spread, as seen by simulations of
Gaussian distributed fluctuations in Appendix C.

A better estimate of the $\langle T_{1}\rangle_{T}$ for each qubit, $Q_{k}$,
in the device can be obtained from a moving average of multiple, $N$,
measurements. We show the evolution of $\langle R\rangle_{T_{0\rightarrow N}}$
using a moving average of the $T_{1}(T_{i})$ measurements, $\langle
T_{1}\rangle_{T_{0\rightarrow N}}$, for each qubit, Fig. 5 (b). The $\langle
R\rangle_{T_{0\rightarrow N}}$ exceeds $R$ $\sim$ 0.8 (i.e., strong
correlation) after $\sim$10 measurements, corresponding to a time exceeding
100 hours.  Approximately 10 independent measurements is sufficient for
fluctuations with magnitude $\sim$0.2$\langle T_{1}\rangle_{T}$ to obtain a
strong correlation, R $\sim$ 0.8, between an estimator (e.g., $\langle
T_{1}\rangle_{T_{0}\rightarrow N}$) and $\langle T_{1}\rangle_{T}$. The
details of R dependence on fluctuation magnitude and number of measurements in
the moving average are discussed more completely in Appendix C.

Autocorrelation between $T_{1}(T_{i})$ and $T_{1}(T_{i-1})$ measurements is an
underlying challenge to fast estimation of $\langle T_{1}\rangle_{T}$.
Evidence of autocorrelation can be seen for example in long term drifts in the
average and short term correlations between $T_{1}$, inset of Fig. 5 (b). On
shorter time scales, our experimental data shows evidence of stronger
autocorrelation frustrating faster accurate estimation of $\langle
T_{1}\rangle_{T}$ and that the fastest $R\sim$ 0.8 can be obtained on order of
1-2 days, see appendices D and H. We conclude that $\langle
T_{1}\rangle_{\omega,t}$ shows promise as a method for faster estimation of
$\langle T_{1}\rangle_{T}$ than repeated $T_{1}(\omega=\omega_{q})$
measurements at only the qubit frequency. Extending the $\langle
T_{1}\rangle_{\omega,t}$ estimator to a set of many qubits, {$Q_{k}$}, in a
device result in larger R, in the same time, compared to relying only on
$T_{1}(\omega_{q})$ measurements for each qubit. The R value simply being a
quantitative single value expression of the high correlation between each
$\langle T_{1}\rangle_{\omega,t}$ and $\langle T_{1}\rangle_{T}$ across the
entire set of qubits.

Figure 5: (a) Comparison of $\langle P_{1}\rangle_{\omega,t}$ at $50\mu$s and
$\langle P_{1}\rangle_{T}$ averaged for $\sim$9 months and evaluated at a
$\tau$ of 53 $\mu$s decay time. $\langle P_{1}\rangle_{\omega,t}$, is averaged
over $\Delta\omega/2\pi$ = 5 MHz after a single measurement that took $\sim$20
minutes. (inset) A scatter plot using $\langle T_{1}\rangle_{T}$’s averaged
over 9 months of measurement as the dependent variable and $\langle
T_{1}\rangle_{\omega,t}$ or $T_{1}$s from a single day. The line is a guide to
the eye showing a 1:1 correlation. (b) The Pearson $R$ dependence on time
averaging of the $T_{1}$’s of the odd numbered qubits up to time, T,  for
three cases: (i) the entire time series (dash), (ii) the time series between
temperature excursion B and D (dash-dot), lettered locations indicated in the
inset, and (iii) the time series between C and D for which no temperature
excursions were recorded (solid). The intermediate time series are shifted in
time index to compare more directly at short times with the full time series.
The differences in $R$ are within the standard deviation calculated for
sampling $T_{1}$ time series with a Gaussian distributed range of values, see
appendix C. (inset) The $T_{1}$ (black, left) and mixing chamber temperature
(blue, right) time series for $Q_{13}$. Spacing of measurements is non-
uniform. The minimum spacing is $\sim$24 hours apart. Each temperature
excursion is labeled with a letter. (c) Pearson correlation, $R$, dependence
on time slice averaging and frequency range, $\Delta\omega$, of the odd
numbered qubits.

It is important to note that our calculations of $\langle
T_{1}\rangle_{\omega,t}$ employ an equal weighting of $P_{1}$ associated with
every frequency bin  and the same choice of $\Delta\omega$ for every qubit.
However, it is not a priori clear that equal weighting is a representative
choice over the $\Delta\omega$ range. For example, how evenly does the
spectral diffusion of each TLS contribute to the $T_{1}$ of the qubit? The
strong correlation of $\langle T_{1}\rangle_{\omega,t}$ with $\langle
T_{1}\rangle_{T}$ with equal weighting suggests that an ergodic-like sampling
of the TLSs near the qubit frequency is a reasonable first approximation.  The
ergodic behavior of the $T_{1}$ estimators is examined more completely in
Appendix F and Appendix G. Central to the question of assigning a $T_{1}$
estimate to any qubit, we observe that $\langle T_{1}\rangle_{T}$ behaves
ergodically for all the qubits despite short term $1/f^{\alpha}$ correlated
behavior (i.e., a constant mean $\langle T_{1}\rangle_{T}$ can be identified).
Assignment of any $T_{1}$ estimate could alternatively be made impossible in
the presence of drift, which is not observed in these qubits, see Appendix B
and Appendix F for further details about weak stationarity and ergodicity.
Furthermore, the strong correlation of $\langle T_{1}\rangle_{T}$ to $\langle
T_{1}\rangle_{\omega,t}$ using only the $P_{1}(\omega,\tau,t)$ spectrum around
the qubit is consistent with a leading hypothesis that the $\langle
T_{1}\rangle_{T}$ is dominated by TLS behavior rather than other stochastic or
static contributions.

##  IV Correlation dependence on frequency and measurement time

A natural question about the estimator $\langle T_{1}\rangle_{\omega,t}$ is,
what are the optimal parameter choices for frequency range $\Delta\omega$, $n$
autocorrelated samples and the spacing in time, $\Delta t=t_{i}-t_{i-1}$, to
obtain sufficiently weakly autocorrelated measurements and a fast, accurate
measure of $\langle T_{1}\rangle_{T}$. Since the optimum choices are presently
not known a priori, we evaluate and plot  $\langle R\rangle_{t_{0\rightarrow
n}}$ versus $\Delta\omega$ and $t_{i}$ in Fig. 5 (c) to guide future
application of this approach. Equal frequency bin weighting of $P_{1}$ is
used. While this order of magnitude choice of $\Delta\omega$ produces a
reasonably good first approximation for correlation across the entire range,
the plot displays several unexplained features (e.g., non-monotonic dependence
on $\Delta\omega$) indicating the unsurprising insufficiency of these two
globally applied parameters (i.e., $\Delta\omega$ and $t$) alone to weight the
frequency contribution of all the qubits and approach $R\sim$ 1\.  Additional
sensitivity analysis in Appendix G also examines correlation between
frequencies and highlights that individual qubits have different sensitivity
to the range sampled, $\Delta\omega$. We see that a wide span of
$\Delta\omega$ produces high $\langle R\rangle_{t_{0\rightarrow n}}$,
comparable or better than $R(T_{i})$ from a single $T_{1}(\omega_{q})$
measurement. We further show that not only is there a strong R correlation
(e.g., linear dependence) but that $\langle T_{1}\rangle_{\omega,t}$
approaches 1:1 quantitative agreement with $\langle T_{1}\rangle_{T}$. The
degree to which a $T_{1}$ estimator, from sampling the nearby frequency space,
is quasi-ergodic and would converge to 1:1 agreement is addressed in much more
detail in Appendix G and Appendix K.

##  V Discussion: implications for process characterization

The strong correlation between $\langle T_{1}\rangle_{\omega,t}$ and $\langle
T_{1}\rangle_{T}$ suggests that long time $T_{1}$ averages might be estimated
relatively rapidly using spectroscopy. This is in contrast to overcoming
correlation times in $T_{1}$ at a single $\omega_{q}$ to obtain a
representative $\langle T_{1}\rangle_{T}$ for the qubit.

Identification of better choices of $\Delta\omega$ and $n$ in this study were
made with pre-knowledge of what $\langle T_{1}\rangle_{T}$ was. These
parameters will have to be chosen without this pre-characterization for future
implementation of this method. Encouragingly, the R dependence on both these
parameters appears to be relatively weak suggesting that a heuristic choice
for a single $\Delta\omega$ and $n$ might be sufficient to obtain useful
estimates (i.e., $R>0.8$) of $\langle T_{1}\rangle_{T}$ for new processes when
using this simple equal weighting approach until improved choices can be
formulated (i.e., different frequency spans for each qubit and or weighted
averaging over frequency).

More specifically we observe that $\mathcal{O}$(10) independent measurements
is sufficient to obtain an R $\sim$ 0.8 or higher, see Appendix C. We
conjecture that one can obtain 10 approximately independent samples, $S$, in a
single scan by sampling at frequency spacings, $\chi$, that are greater than
the autocorrelation frequency width (i.e., a frequency spacing where
correlation drops below $\sim$0.2). In this work, we found the correlation to
become weak for $\mathcal{O}$(1 MHz), see Appendix G. Then by this heuristic,
a single spectroscopy scan would require a $\Delta\omega$ =
$\frac{(S-1)}{2}\chi$, where $S$ = 10 for the target of $R\sim$0.8. We assume
one of the measurements is done at the qubit frequency, $T_{1}(\omega_{q})$,
so for a $\chi\sim$ 1 MHz, a scan from $\pm$4.5 MHz would be suggested by such
a heuristic. Extra $n$ measurements can be obtained by waiting longer than the
autocorrelation time. The autocorrelation width, furthermore, can be evaluated
in the same scan as that used for the $\langle T_{1}\rangle_{T}$ estimate as
long as a sufficiently wide range is sampled. Alternatively a second scan can
be taken if the initial $\Delta\omega$ guess was too small.

Empirically we see diminishing gains in using ever larger $\Delta\omega$.
Further research is needed to guide better limits on $\Delta\omega$ beyond the
operational observation that $S\sim\mathcal{O}(10)$ produces a quasi-ergodic
result for qubits with $\langle T_{1}\rangle_{T}$ in the range of 10-200
$\mu$s, see Appendix G for more details on quasi-ergodicity. Since we do find
$\sim$1:1 agreement using a relatively small $\Delta\omega\sim$ 10 MHz for the
$\sim$9 month time series and we observe that the distribution of
$T_{1}(\omega_{q},T_{i})$ produces a constant standard deviation, see Appendix
B, rather than growing (e.g., proportional to a random walk
$\propto\sqrt{t}$), we speculate that optimal $\Delta\omega$ is bounded rather
than growing indefinitely from spectral diffusion processes. Notably Klauder
et al. calculate that dipole coupled ensembles that are proposed for TLS
spectral diffusion Black and Halperin , will produce a truncated linewidth
Klauder and Anderson (1962).

##  VI Conclusion

In this work, we probe the temporal and spectral dynamics of superconducting
qubit relaxation times. We study these dynamics in high coherence, single-
junction transmons by developing a technique for energy relaxation
spectroscopy of defect TLSs via the AC Stark effect. Our technique requires no
additional hardware resources and can be easily sped up further by integration
with reset schemes. Autocorrelation of $T_{1}$ frustrates rapid
characterization of the long-time average $\langle T_{1}\rangle_{T}$ and
therefore accurate characterization of devices. Our analysis of the dynamics
identifies a strong correlation between $\langle T_{1}\rangle_{T}$ and its
short time average over the local frequency span, $\langle
T_{1}\rangle_{\omega,t}$. The strong correlation of $\langle T_{1}\rangle_{T}$
with $\langle T_{1}\rangle_{\omega,t}$ is also consistent with a TLS dominated
$T_{1}$ that quasi-ergodically samples the qubit local frequency neighborhood
in contrast to static or uncorrelated stochastic processes. This work opens up
several new promising directions for rapid process characterization and
evaluation of device stability.

##  VII Data Availability

The data that support the findings of this study are available from the
corresponding author on reasonable request.

##  VIII References

## References

  * Zhang  _et al._ (2020) E. J. Zhang, S. Srinivasan, N. Sundaresan, D. F. Bogorin, Y. Martin, J. B. Hertzberg, J. Timmerwilke, E. J. Pritchett, J.-B. Yau, C. Wang, _et al._ , “High-fidelity superconducting quantum processors via laser-annealing of transmon qubits,” arXiv preprint arXiv:2012.08475 (2020).
  * Arute  _et al._ (2019) F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , “Quantum supremacy using a programmable superconducting processor,” Nature 574, 505–510 (2019).
  * Nakamura, Pashkin, and Tsai (1999) Y. Nakamura, Y. A. Pashkin, and J. S. Tsai, “Coherent control of macroscopic quantum states in a single-cooper-pair box,” nature 398, 786–788 (1999).
  * Hong  _et al._ (2020) S. S. Hong, A. T. Papageorge, P. Sivarajah, G. Crossman, N. Didier, A. M. Polloreno, E. A. Sete, S. W. Turkowski, M. P. da Silva, and B. R. Johnson, “Demonstration of a parametrically activated entangling gate protected from flux noise,” Physical Review A 101, 012302 (2020).
  * Kandala  _et al._ (2020) A. Kandala, K. Wei, S. Srinivasan, E. Magesan, S. Carnevale, G. Keefe, D. Klaus, O. Dial, and D. McKay, “Demonstration of a high-fidelity cnot for fixed-frequency transmons with engineered zz suppression,” arXiv preprint arXiv:2011.07050 (2020).
  * Hashim  _et al._ (2020) A. Hashim, R. K. Naik, A. Morvan, J.-L. Ville, B. Mitchell, J. M. Kreikebaum, M. Davis, E. Smith, C. Iancu, K. P. O’Brien, I. Hincks, J. J. Wallman, J. Emerson, and I. Siddiqi, “Randomized compiling for scalable quantum computing on a noisy superconducting quantum processor,” (2020), arXiv:2010.00215 [quant-ph] .
  * Foxen  _et al._ (2020) B. Foxen, C. Neill, A. Dunsworth, P. Roushan, B. Chiaro, A. Megrant, J. Kelly, Z. Chen, K. Satzinger, R. Barends, F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, S. Boixo, D. Buell, B. Burkett, Y. Chen, R. Collins, E. Farhi, A. Fowler, C. Gidney, M. Giustina, R. Graff, M. Harrigan, T. Huang, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, P. Klimov, A. Korotkov, F. Kostritsa, D. Landhuis, E. Lucero, J. McClean, M. McEwen, X. Mi, M. Mohseni, J. Y. Mutus, O. Naaman, M. Neeley, M. Niu, A. Petukhov, C. Quintana, N. Rubin, D. Sank, V. Smelyanskiy, A. Vainsencher, T. C. White, Z. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis (Google AI Quantum), “Demonstrating a continuous set of two-qubit gates for near-term quantum algorithms,” Phys. Rev. Lett. 125, 120504 (2020).
  * Müller  _et al._ (2015) C. Müller, J. Lisenfeld, A. Shnirman, and S. Poletto, “Interacting two-level defects as sources of fluctuating high-frequency noise in superconducting circuits,” Physical Review B 92, 035442 (2015).
  * Paladino  _et al._ (2014) E. Paladino, Y. M. Galperin, G. Falci, and B. L. Altshuler, “ noise: Implications for solid-state quantum information” Rev. Mod. Phys. 86, 361–418 (2014).
  * Weissman (1988) M. B. Weissman, “$\frac{1}{f}$ noise and other slow, nonexponential kinetics in condensed matter,” Rev. Mod. Phys. 60, 537–571 (1988).
  * Klimov  _et al._ (2018) P. Klimov, J. Kelly, Z. Chen, M. Neeley, A. Megrant, B. Burkett, R. Barends, K. Arya, B. Chiaro, Y. Chen, A. Dunsworth, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, T. Huang, E. Jeffrey, E. Lucero, J. Mutus, O. Naaman, C. Neill, C. Quintana, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. White, S. Boixo, R. Babbush, V. Smelyanskiy, H. Neven, and J. Martinis, “Fluctuations of Energy-Relaxation Times in Superconducting Qubits,” Phys. Rev. Lett. 121, 090502 (2018).
  * Kandala  _et al._ (2019) A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, “Error mitigation extends the computational reach of a noisy quantum processor,” Nature 567, 491–495 (2019).
  * Burnett  _et al._ (2019) J. J. Burnett, A. Bengtsson, M. Scigliuzzo, D. Niepce, M. Kudra, P. Delsing, and J. Bylander, “Decoherence benchmarking of superconducting qubits,” npj Quantum Information 5, 1 (2019).
  * Schlör  _et al._ (2019) S. Schlör, J. Lisenfeld, C. Müller, A. Bilmes, A. Schneider, D. P. Pappas, A. V. Ustinov, and M. Weides, “Correlating decoherence in transmon qubits: Low frequency noise by single fluctuators,” Physical review letters 123, 190502 (2019).
  * Phillips (1972) W. Phillips, “Tunneling states in amorphous solids,” Journal of Low Temperature Physics 7, 351–360 (1972).
  * Müller, Cole, and Lisenfeld (2019) C. Müller, J. H. Cole, and J. Lisenfeld, “Towards understanding two-level-systems in amorphous solids – Insights from quantum circuits,” Rep. Prog. Phys. 82, 124501 (2019), arXiv: 1705.01108.
  * Martinis  _et al._ (2005) J. M. Martinis, K. B. Cooper, R. McDermott, M. Steffen, M. Ansmann, K. Osborn, K. Cicak, S. Oh, D. P. Pappas, R. W. Simmonds, _et al._ , “Decoherence in josephson qubits from dielectric loss,” Physical review letters 95, 210503 (2005).
  * Grabovskij  _et al._ (2012) G. J. Grabovskij, T. Peichl, J. Lisenfeld, G. Weiss, and A. V. Ustinov, “Strain tuning of individual atomic tunneling systems detected by a superconducting qubit,” Science 338, 232–234 (2012).
  * Barends  _et al._ (2013) R. Barends, J. Kelly, A. Megrant, D. Sank, E. Jeffrey, Y. Chen, Y. Yin, B. Chiaro, J. Mutus, C. Neill, P. O’Malley, P. Roushan, J. Wenner, T. C. White, A. N. Cleland, and J. M. Martinis, “Coherent josephson qubit suitable for scalable quantum integrated circuits,” Phys. Rev. Lett. 111, 080502 (2013).
  * Lisenfeld  _et al._ (2019) J. Lisenfeld, A. Bilmes, A. Megrant, R. Barends, J. Kelly, P. Klimov, G. Weiss, J. M. Martinis, and A. V. Ustinov, “Electric field spectroscopy of material defects in transmon qubits,” npj Quantum Information 5, 1–6 (2019).
  * (21) J. L. Black and B. I. Halperin, “Spectral diffusion, phonon echoes, and saturation recovery in glasses at low temperatures,” Physical Review B 16, 2879–2895.
  * Stehlik  _et al._ (2021) J. Stehlik, D. Zajac, D. Underwood, T. Phung, J. Blair, S. Carnevale, D. Klaus, G. Keefe, A. Carniol, M. Kumph, _et al._ , “Tunable coupling architecture for fixed-frequency transmons,” arXiv preprint arXiv:2101.07746 (2021).
  * Gambetta  _et al._ (2006) J. Gambetta, A. Blais, D. I. Schuster, A. Wallraff, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, “Qubit-photon interactions in a cavity: Measurement-induced dephasing and number splitting,” Phys. Rev. A 74, 042318 (2006).
  * Majer  _et al._ (2007) J. Majer, J. Chow, J. Gambetta, J. Koch, B. Johnson, J. Schreier, L. Frunzio, D. Schuster, A. A. Houck, A. Wallraff, _et al._ , “Coupling superconducting qubits via a cavity bus,” Nature 449, 443–447 (2007).
  * Simmonds  _et al._ (2004) R. W. Simmonds, K. M. Lang, D. A. Hite, D. P. Pappas, and J. M. Martinis, “Decoherence in Josephson Qubits from Junction Resonances,” Physical Review Letters 93, 077003 (2004), arXiv: cond-mat/0402470.
  * Magesan and Gambetta (2020) E. Magesan and J. M. Gambetta, “Effective hamiltonian models of the cross-resonance gate,” Physical Review A 101, 052308 (2020).
  * Schneider  _et al._ (2018) A. Schneider, J. Braumüller, L. Guo, P. Stehle, H. Rotzinger, M. Marthaler, A. V. Ustinov, and M. Weides, “Local sensing with the multilevel ac stark effect,” Phys. Rev. A 97, 062334 (2018).
  * Ramsey (1950) N. F. Ramsey, “A Molecular Beam Resonance Method with Separated Oscillating Fields,” Phys. Rev. 78, 695–699 (1950).
  * Egger  _et al._ (2018) D. J. Egger, M. Werninghaus, M. Ganzhorn, G. Salis, A. Fuhrer, P. Mueller, and S. Filipp, “Pulsed reset protocol for fixed-frequency superconducting qubits,” Physical Review Applied 10, 044030 (2018).
  * Yan  _et al._ (2016) F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, _et al._ , “The flux qubit revisited to enhance coherence and reproducibility,” Nature communications 7, 1–9 (2016).
  * Kurter  _et al._ (2020) C. Kurter, M. Sandberg, V. Adiga, Z. Minev, A. Finck, A. Eddins, R. Shelby, M. Brink, and J. Chow, “Detection of non-equilibrium quasiparticles in charge sensitive transmons,” Bulletin of the American Physical Society 65 (2020).
  * Klauder and Anderson (1962) J. R. Klauder and P. W. Anderson, “Spectral diffusion decay in spin resonance experiments,” Phys. Rev. 125, 912–932 (1962).
  * Krishnan (2016) V. Krishnan,  _Probability and Random Processes_ , 2nd ed. (Wiley, 2016).
  * Milotti (2002) E. Milotti, “1/f noise: a pedagogical review,” (2002), arXiv:0204033 [physics.class-ph] .
  * ANSCOMBE and GLYNN (1983) F. J. ANSCOMBE and W. J. GLYNN, “Distribution of the kurtosis statistic b2 for normal samples,” Biometrika 70, 227–234 (1983), https://academic.oup.com/biomet/article-pdf/70/1/227/661551/70-1-227.pdf .
  * D’agostino, Belanger, and Jr. (1990) R. B. D’agostino, A. Belanger, and R. B. D. Jr., “A suggestion for using powerful and informative tests of normality,” The American Statistician 44, 316–321 (1990), https://www.tandfonline.com/doi/pdf/10.1080/00031305.1990.10475751 .
  * (37) J. G. MacKinnon, “Approximate asymptotic distribution functions for unit-root and cointegration tests,”  12, 167–176, publisher: [American Statistical Association, Taylor & Francis, Ltd.].
  * MacKinnon (2010) J. G. MacKinnon, “Critical Values For Cointegration Tests,” Working Paper 1227 (Economics Department, Queen’s University, 2010).
  * Harris (1992) R. Harris, “Testing for unit roots using the augmented dickey-fuller test: Some issues relating to the size, power and the lag structure of the test,” Economics Letters 38, 381–386 (1992).
  * Singh, Ghosh, and Adhikari (2018) R. Singh, D. Ghosh, and R. Adhikari, “Fast bayesian inference of the multivariate ornstein-uhlenbeck process,” Phys. Rev. E 98, 012136 (2018).
  * (41) M. T. Madzik, T. D. Ladd, F. E. Hudson, K. M. Itoh, A. M. Jakob, B. C. Johnson, J. C. McCallum, D. N. Jamieson, A. S. Dzurak, A. Laucht, and A. Morello, “Controllable freezing of the nuclear spin bath in a single-atom spin qubit,”  6, eaba3442, publisher: American Association for the Advancement of Science, 32937454 .
  * (42) X. Brokmann, J.-P. Hermier, G. Messin, P. Desbiolles, J.-P. Bouchaud, and M. Dahan, “Statistical aging and nonergodicity in the fluorescence of single nanocrystals,”  90, 120601.
  * Reif (1965) F. Reif,  _Fundamentals of statistical and thermal physics_ , McGraw-Hill series in fundamentals of physics (McGraw-Hill, New York, 1965).
  * Wald and Wolfowitz (1943) A. Wald and J. Wolfowitz, “An exact test for randomness in the non-parametric case based on serial correlation,” The Annals of Mathematical Statistics 14, 378–388 (1943).
  * de Graaf  _et al._ (2020) S. E. de Graaf, L. Faoro, L. B. Ioffe, S. Mahashabde, J. J. Burnett, T. Lindström, S. E. Kubatkin, A. V. Danilov, and A. Y. Tzalenchuk, “Two-level systems in superconducting quantum devices due to trapped quasiparticles,” Science Advances 6 (2020), 10.1126/sciadv.abc5055, https://advances.sciencemag.org/content/6/51/eabc5055.full.pdf .
  * Herzog and Hahn (1956) B. Herzog and E. L. Hahn, “Transient Nuclear Induction and Double Nuclear Resonance in Solids,” Physical Review 103, 148–166 (1956).
  * Witzel  _et al._ (2012) W. M. Witzel, M. S. Carroll, L. Cywiński, and S. Das Sarma, “Quantum decoherence of the central spin in a sparse system of dipolar coupled spins,” Phys. Rev. B 86, 035452 (2012).

##  IX Acknowledgements

We acknowledge technical support on the ibmq_almaden device from the IBM
Quantum deployment team. Additional insightful discussions, suggestions and
assistance came from Nick Bronn, Andrew Cross, Oliver Dial, Doug McClure,
Easwar Magesan, Hasan Nayfeh, James Raferty, Martin Sandberg, Srikanth
Srinivasan, Neereja Sundaresan, Jerry Tersoff, Ben Fearon, Karthik
Balakrishnan, James Hannon and Jerry Chow.

MC also acknowledges support from Princeton Plasma Physics Laboratory through
the Department of Energy Laboratory Directed Research and Development program
and contract number DE-AC02-09CH11466 to complete parts of the analysis and
manuscript.

##  X Author Contributions

S. R. and A. K. developed the technique with contributions from M. C., I. L.
and P. J. M. C. and S. R. performed the experiments. M.C., S.R and A. K.
analyzed the data. M. C., S. R. and A.K, wrote the manuscript with feedback
from the other authors.

##  XI Competing Interests

The authors declare that elements of this work will be included  in patents
filed by the International Business Machines Corporation with the US Patent
and Trademark office. The authors declare no other financial or non/financial
competing interests in relation ot this published work.

##  Appendix A Spectroscopy of odd numbered qubits in device

The spectrally and temporally resolved dynamics of $T_{1}$ for all the odd
numbered qubits in the device are provided for reference in Fig. 6. The data
was taken under the same conditions and at the same time as Fig. 4 shown in
the main body of the paper. Care was taken to avoid frequency collisions both
over the range of qubit frequency shift and the placement of $\omega_{s}$. We
note that the odd qubits do not have direct connectivity with each other,
making cross talk effects negligible for these measurements.

Figure 6: (a)-(i) Spectroscopy for odd qubits 1-19 excluding qubit 15  in
respective order. Explicitly: (a) $Q_{1}$, (b) $Q_{3}$, (c) $Q_{5}$, (d)
$Q_{7}$, (e) $Q_{9}$, (f) $Q_{11}$, (g) $Q_{13}$, (h) $Q_{17}$, (i) $Q_{19}$.
$Q_{15}$ is shown in the main text.

##  Appendix B Stationarity of $T_{1}$ time series

An important question for $T_{1}$ analysis of qubits is whether there is in
fact a representative $T_{1}$ mean and variance that can be assigned to each
qubit. In time series analysis, this concept is described as weak stationarity
Krishnan (2016). Weak stationarity is also a necessary condition for
ergodicity Milotti (2002).

The mean and standard deviations qualitatively appear relatively constant. A
histogram of the $T_{1}$s over 9 months for $Q_{13}$ is plotted for a
representative case along with a fit to a normal distribution Fig. 7 (a).
Visually, the distributions appear relatively normal indicating that the mean
is not drifting substantially relative to the variance, while the skew and
kurtosis are also relatively small. In other qubits, similar near normally
distributed $T_{1}$ fluctuations are observed, see table 2. The skew and
kurtosis are frequently not discernible statistically from a normal
distribution ANSCOMBE and GLYNN (1983); D’agostino, Belanger, and Jr. (1990)
with notable exceptions such as $Q_{7}$, which has a tight distribution around
the mean and some instances for which the skew was distinguishably larger than
normal in $Q_{3}$ and $Q_{17}$.

Qubit |  $\mu_{(1)}$ ($\mu$s) |  $\mu_{(2)}$ ($\mu$s) |  $\mu_{(3)}$ |  $\mu_{(4)}$ |  p-value $\mu_{(3)}$ test |  p-value $\mu_{(4)}$ test
---|---|---|---|---|---|---
1 | 86.1 | 28.3 | -0.12 | -0.32 | 0.37 | 0.23
3 | 97.2 | 26.4 | -0.44 | 0.30 | 0.002 | 0.25
5 | 107.0 | 28.5 | -0.25 | 0.23 | 0.07 | 0.32
7 | 31.4 | 5.6 | -0.83 | 1.92 |  $\sim$0 |  $\sim$0
9 | 107.9 | 29.3 | -0.18 | -0.03 | 0.20 | 0.94
11 | 68.4 | 17.7 | -0.11 | -0.09 | 0.42 | 0.88
13 | 106 | 32.7 | -0.09 | 0.22 | 0.50 | 0.34
15 | 98.8 | 28.4 | -0.19 | -0.32 | 0.17 | 0.23
17 | 104.1 | 28.3 | -0.28 | -0.32 | 0.04 | 0.22
19 | 113.4 | 35.1 | -0.14 | -0.24 | 0.29 | 0.44
Table 2:  Calculated first four moments of the $T_{1}$ distributions for all
the odd qubits. The mean, standard deviation, skew and kurtosis are indicated
as $\mu_{(1)}$, $\mu_{(2)}$, $\mu_{(3)}$, and $\mu_{(4)}$, respectively, in
the table. The table also contains the p-value results of skew and kurtosis
tests for normality ANSCOMBE and GLYNN (1983); D’agostino, Belanger, and Jr.
(1990). The null hypothesis is that the test distribution comes from a normal
distribution. A value $\geq$ 0.05 is often used as the threshold for the null
hypothesis to be accepted (i.e., the distribution is normal-like for that
tested moment).

We also show the moving average of the standard deviation of the $T_{1}$
distributions for each qubit in Fig. 7 (b). The standard deviation of each
qubit is normalized to its respective mean, $\mu$. We see that the general
trend is for the standard deviations to converge towards their mean,
$\sigma\sim\mathcal{O}$($\mu$). Drift in $\sigma$ is small relative to $\mu$.
The $\sigma$ behaves weakly stationary rather than, for example, random walk
like (i.e., $\sigma\propto\sqrt{t}$).

Figure 7:  (a) Probability density of $T_{1}$ fluctuations for a
representative qubit, $Q_{13}$. The normalized histograms are derived from the
$T_{1}$ time series. (b) Normalized cumulative moving average of the standard
deviation of the $T_{1}$ fluctuations as a function of time index spanning
$\sim$9 months.

To more rigorously test the weak stationarity, we apply an augmented Dickey
Fuller (ADF) test MacKinnon ; MacKinnon (2010); Harris (1992) to the
timeseries. The ADF test is a commonly used test of weak stationarity, testing
both drift and constant variance. In our case we are most concerned with how
to test whether the variance is stationary after observing that the mean is
stationary, at least within $\mathcal{O}(\sigma_{T1})$ over 9 months. Random
walks are the canonical case for which the mean is constant but the timeseries
will be non-stationary because the variance grows with time. In particular,
the ADF tests the likelihood of a unit root difference equation regression
with the timeseries in question. Unit root is synonomous with random walk
behavior. ADF uses the following parameterized model:

$y_{i}=\alpha y_{i-1}+\sum_{j=1}^{\rho}\beta\Delta y_{i-j}+\delta+\gamma
i+\epsilon_{i}$ (7)

where i is the time step index, $\alpha$ is the root, the sum is over
additional lag terms, $\delta$ is a constant offset, $\gamma$ is the slope of
a linear trend and $\epsilon$ is a random error term that is normally
distributed with a standard deviation of $\sigma_{\epsilon}$. The lag terms
importantly account for effects of serial correlation (e.g., non-Markovian
behavior expected in 1/f noise), while the drift term can be used to establish
’trend stationary’ behavior. The null hypothesis, $H_{0}$, is that there is a
unit root. If the time series is stationary, the ADF $H_{0}$ must be rejected.

We test using $\gamma$ = 0, which will accept the null hypothesis, $H_{0}$,
for either the case of a unit root (e.g., random walk) or for an
$\mathcal{O}(\sigma_{\epsilon})$ non-stationary trend summed with a random
error term. Table 3 shows the results of the tests. All $T_{1}$ timeseries
reject the ADF $H_{0}$. All $T_{1}$ timeseries are therefore consistent with
being weakly stationary.

Qubit |  t-stat |  p-value
---|---|---
1 | -9.3 | $9.1\times 10^{-14}$
3 | -6.14 | $1.0\times 10^{-6}$
5 | -10.1 | $1.6\times 10^{-15}$
7 | -11.3 | $5.2\times 10^{-18}$
9 | -12.1 | $2.5\times 10^{-19}$
11 | -8.5 | $5.8\times 10^{-12}$
13 | -6.0 | $1.9\times 10^{-6}$
15 | -9.9 | $3.0\times 10^{-15}$
17 | -9.6 | $1.4\times 10^{-14}$
19 | -4.1 | $7.3\times 10^{-3}$
Table 3: ADF tests for odd qubit $T_{1}$ time series using the python
statsmodel library. The tests use a Bayesian information criteria to determine
the number of lags. P-values of greater than 0.05 are the typical statistical
threshold for acceptance of the $H_{0}$, null hypothesis. The table values are
for $\gamma$ = 0.

##  Appendix C Pearson correlation dependence on standard deviation and
sampling

We calculate an expectation of how many uncorrelated measurements are
necessary to achieve strong $R$ correlation (i.e., $R\sim$ 0.8). We calculate
a Pearson $\langle R(n_{q},\sigma_{mk})\rangle_{T_{0\rightarrow N}}$ for
$n_{q}$ being a ten qubit device, $\sigma_{mk}$ being the $T_{1}$ standard
deviation for the $k^{th}$ qubit and where $N$ is the number of uncorrelated
$T_{1}$ measurements for each of the ten qubits. The simulated $\langle
R\rangle_{T_{0\rightarrow N}}$ calculation has a long-time average $T_{1}$
assigned to each of the 10 qubits, $\langle T_{1}\rangle_{T,k}$. For every
$k$-th qubit, $\langle T_{1}\rangle_{T,k}$ is chosen randomly from a normal
distribution with a mean of 100 $\mu$s and a standard deviation of 10 $\mu$s.
$\langle T_{1}\rangle_{T,k}$ simulates a long time stationary $\langle
T_{1}\rangle_{T}$ for each qubit. This standard deviation is representative of
process variation in the qubits of each simulated device.

We then simulate a sequence of measurements. Each measurement obtains an
instantaneous $T_{1}$ for each qubit in the device. The $T_{1}$ measurement is
chosen from a normal distribution with a  standard deviation of $\sigma_{mk}$
= 0.2$\langle T_{1}\rangle_{T,k}$, which is of similar magnitude to what is
observed in the device. The $\sigma_{mk}$ is a measure of the time fluctuating
$T_{1}$ centered around $\langle T_{1}\rangle_{T,k}$, in contrast to the 10
$\mu$s above, which is the the variability of the static $\langle
T_{1}\rangle_{T,k}$ centered around $\langle T_{1}\rangle_{T,n_{q}}$ =
$\frac{1}{n_{q}}\sum_{k}^{n_{q}}\langle T_{1}\rangle_{T,k}$.

We can obtain an estimate of the Pearson $R$ correlation between $\langle
T_{1}\rangle_{T,k}$ and a single measurement instance $T_{1,k}$ for all the
qubits. The effect of multiple uncorrelated device measurements is then
simulated by repeating the device measurement and updating the average of the
$T_{1,k}$ with all previous $T_{1}$ measurements. In order to simulate the
dependence of $\langle R\rangle_{T_{0\rightarrow N}}$ on the number of
uncorrelated measurements $N$, we first define an expectation value,
$\langle\langle R\rangle_{T_{0\rightarrow N}}\rangle_{Ndev}$, averaged over
$N_{dev}$ = 200 devices with 4, 10, 20 and 50 qubits.

$\langle\langle R\rangle_{T_{0\rightarrow N}}\rangle_{Ndev}$ approaches unity
with increasing uncorrelated measurements, as seen in Fig. 8. This indicates
that the averaging of uncorrelated measurements increasingly produces an
accurate estimate of $\langle T_{1}\rangle_{T,k}$. The standard deviation,
$\sigma_{R}$, of $\langle\langle R\rangle_{T_{0\rightarrow N}}\rangle_{Ndev}$
is also shown in the inset, Fig. 8. We can see that initial $\langle
R\rangle_{T_{0\rightarrow N}}$ values can be very low, and around $N$ = 10
uncorrelated measurements are required to obtain $R\sim$ 0.8, the correlation
obtained from the frequency averaging of $T_{1}$ discussed in the main text.

To provide additional insight into the $\langle R\rangle_{T_{0\rightarrow N}}$
dependence on $n_{q}$, $\sigma_{mk}$ and $N$ measurements, we derive an
analytic expression for $\langle R(n_{q},\sigma_{mk})\rangle_{T_{0\rightarrow
N}}$. We express R in terms of differences, $\delta X$, from $\langle
T_{1}\rangle_{T,n_{q}}$. Explicitly writing out the first few terms of
$\langle R(n_{q},\sigma_{mk})\rangle$’s sums for a multiqubit device:

$\langle R(n_{q},\sigma_{mk})\rangle=\frac{\delta X_{1}(\delta
X_{1}+\sigma_{m1})+\delta X_{2}(\delta X_{2}+\sigma_{m2})+...}{\sqrt{[\delta
X_{1}^{2}+\delta X_{2}^{2}+...][(\delta X_{1}+\sigma_{m1})^{2}+...]}}$ (8)

where we define $\delta X_{k}=\langle T_{1}\rangle_{T,k}-\langle
T_{1}\rangle_{T,n_{q}}$ and $\sigma_{mk}$ is the standard deviation for
measurements of $T_{1,k}$ for the $k^{th}$ qubit in the device.

We parameterize the standard deviation of the $T_{1}$ measurement with a
scaling constant $\alpha$ for the $k^{th}$ qubit as
$\sigma_{mk}=\alpha_{i}\langle T_{1}\rangle_{T,k}$. Likewise, we assume that
the $\langle T_{1}\rangle_{T,k}$ in a device are normally distributed and may
be parameterized with $\beta_{i}$ as $\sigma_{di}=\beta_{i}\langle
T_{1}\rangle_{n_{q}}$ (i.e., $\langle T_{1}\rangle_{T,k}=\langle
T_{1}\rangle_{n_{q}}+\sigma_{dk}$).

We solve for $\langle R\rangle_{T_{0\rightarrow N}}$ for a particular device
instance with $\langle T_{1}\rangle_{T,k}$, in the simpler case that
$\alpha_{k}$ = $\alpha$ (i.e., $\sigma_{mk}$ = $\sigma_{m}$), and with the
device defined by $\vec{\beta}$ = $\\{\beta_{1},\beta_{2},...\\}$. The
$\vec{\beta}$ defines the $\langle T_{1}\rangle_{T}$s of a device instance.
Multiplying and reorganizing the first few terms:

$\langle R(n_{q},\sigma_{m})\rangle=\frac{(\delta X_{1}^{2}+\delta
X_{2}^{2}+...)+\sigma_{m}(\delta X_{1}+\delta X_{2}+...)}{\sqrt{[\delta
X_{1}^{2}+\delta X_{2}^{2}+...][(\delta X_{1}+\sigma_{m})^{2}+...]}}$ (9)

Substituting $\delta X_{k}+\sigma_{m}\rightarrow\delta
X_{k}+\frac{\sigma}{\sqrt{N}}$ to express the moving average dependence of
$T_{1}$ on N measurements, we obtain:

$\langle
R(n_{q},\sigma_{m},N)\rangle=\frac{(\sum_{i}^{n_{q}}\beta_{i}^{2})+\frac{\alpha}{\sqrt{N}}(\sum_{i}^{n_{q}}\beta_{i})}{\sqrt{(\sum_{i}^{n_{q}}\beta_{i}^{2})(\sum_{i}^{n_{q}}(\beta_{i}+\frac{\alpha}{\sqrt{N}})^{2})}}$
(10)

To compare to the simulations in Fig. 8, we numerically sample many instances
of eqn. 10, $N_{dev}$ instances of $\vec{\beta}$. The $\vec{\beta}$ are
assumed to have a normal distribution. The expression shows good agreement
with the full numerical simulation described above, Fig. 8. The expression 10
provides the quantitative dependence of how increasing $N$ reduces the
uncertainty in R through reducing the uncertainty in each $T_{1,k}$ estimator.
The uncertainty is reduced by the familiar $\sqrt{N}$ dependence resulting in
$\lim_{N\rightarrow\infty}\langle R\rangle_{T_{0\rightarrow N}}\rightarrow 1$.
Increasing the number of qubits in the device also can be used to reduce
$\sigma_{R}$ for a single device measurement, playing a similar role of
averaging $\vec{\beta}$ instead of alternatively measuring many devices.

Figure 8: Simulated Pearson $R$ dependence  on number of qubits in the device
and on $N$ uncorrelated device measurements of $T_{1}$ for each qubit in the
device. The simulation is done for 200 simulated devices to find an
expectation value $\langle R\rangle$. The semi-analytic expression in the text
agrees well to the 4 qubit simulation. The inset shows the dependence of the
standard deviation of $R$ on $N$ and qubit number.

##  Appendix D Autocorrelation of $T_{1}$($\omega_{q},T_{i})$

We now ask the question: How long does it take to obtain an uncorrelated
measurement of $T_{1}$? To address this, we calculate the autocorrelation for
the $\sim$9 month $T_{1}$ time series, depicted in Fig 9. The autocorrelation
is detrended using the mean and normalized using the estimated variance.  The
period between sampling is approximately 24 hours with some variability (i.e.,
hours). All qubits have some time correlation in the first few measurements
(i.e., 1-2 days). For short time behavior see appendix H. Some qubits also
show weak and decaying autocorrelation at longer times. For example, $Q_{13}$
and $Q_{1}$’s time series appear to have long time trends, Fig. 9. The
decaying oscillatory autocorrelation is qualitatively consistent with a mean
reverting time series, for example, an Ornstein-Uhlenbeck process commonly
used to simulate pink noise Singh, Ghosh, and Adhikari (2018); and is also
consistent with the weak stationarity of the same time series found in
Appendix B. In general, these longer timescale autocorrelations highlight a
challenge in extrapolating $\langle T_{1}\rangle_{T}$ from sampling $T_{1}$ at
a single frequency at early times in the time evolution.

Figure 9: (a), (b), (c) and (d) Autocorrelation as a function of measurement
lag time for qubits 1, 7, 13 and 15, respectively. The autocorrelation is
detrended using the mean of the time series and normalized using the estimated
variance. The insets show the $T_{1}$ time series from which the
autocorrelation was calculated and the corresponding temperature of the mixing
chamber plate.

##  Appendix E Thermal excursion effects on autocorrelation

There were mixing chamber plate temperature excursions during the duration of
the $\sim$9 month $T_{1}(T_{i})$ measurement time series, Fig. 9. To address
doubts about the impact of the temperature excursions, we show autocorrelation
for representative qubits for the longest time series in which there are no
thermal excursions for comparison, Fig 10. The qualitative behavior and
magnitudes of the long time correlations are similar to those cases with the
temperature excursions. Exact quantitative agreement does change but this
would also be expected from truncating the time series at any different
starting time.

In Appendix D we show that the autocorrelation time is order of days. This is
consistent with observing no strong effects of the temperature excursions on
the weak stationarity of the time series or long term behavior of the
distributions. That is, a temperature excursion would be expected to have a
very localized effect on a $\sim$9 month time series behavior. The
measurements using spectroscopy were done when the temperature was stable.

We conclude that the temperature excursions don’t effect the conclusions of
the paper.

Figure 10: (a) and (b) autocorrelation as a function of measurement lag time
for qubits 1 and 13, respectively, for an $\sim$2890 hour time period between
temperature excursions. (c) and (d) are the autocorrelations using the entire
time series, $\sim$2890 hours, for the corresponding qubits 1 and 13. The
autocorrelation is detrended using the mean of the time series and normalized
using the estimated variance.

##  Appendix F Ergodicity of $T_{1}(\omega_{q},T_{i})$ time series

A foundational question is whether a $T_{1}$ time series behaves ergodically.
That is, the time series $T_{1}(T_{i})$ both converges to a $\langle
T_{1}\rangle$ (e.g., after sufficient time lag) and the pair correlations are
well behaved (e.g., decay at long lag). Ergodicity is not guaranteed, drift or
$1/f^{\alpha}$ behavior being illustrative reasons to doubt whether a reliable
long time estimator of a physical property can be obtained. Since breaking of
ergodicity can signal physical phenomena of interest such as switching between
phases (i.e., isolated systems of an ensemble) with distinct mean values
including special cases of spectral diffusion Madzik  _et al._ ; Brokmann  _et
al._ , the establishment of whether the $T_{1}$ time series behaves
ergodically represents an important step in clarifying the dynamics of the
$T_{1}$ fluctuations.

The ergodic assumption is that given sufficient time a system will visit
through all the accessible states (i.e., values) available to it. Such a
sufficiently long time series trajectory can then be divided into $k$
independent subsets to form an ensemble of $k$ new systems that should
represent the statistical behavior of the original system at any given time
index, $m$ of the $k$ systems, Fig. 11 (a) Reif (1965). The new time series of
the systems are re-indexed, $m$, with equal lengths. An ensemble average is
defined by selecting from the same time index for all equally sized subsets,

$\\{T_{1}(T_{m})\\}=\frac{1}{k}\sum_{k}T_{1}(T_{m})^{(k)}\mid m=constant$ (11)

and for an ergodic system,

$\\{T_{1}(T_{m})\\}=\langle T_{1}\rangle\simeq\langle T_{1}\rangle_{T}$ (12)

when averaging across the ensemble of newly defined systems for any time index
$i$. We test the $\sim$9 month $T_{1}$ time series for ergodic behavior.
Autocorrelation of the different qubit time series shows some correlation over
1-2 days (i.e., first several measurements), see Appendix D. We examine a
range of ensembles of size, $2\leq k_{max}\leq 40$ containing $160\geq
m_{max}\geq 5$ time points, respectively.

Ensemble averages are distributed around the time average $\langle
T_{1}\rangle_{T}$, Fig. 11 (c). We test the likelihood that the
$\\{T_{1}(T_{m})\\}$ is statistically indistinguishable from $\langle
T_{1}\rangle_{T}$ using a t-test comparison of the mean values of the two
$T_{1}$ estimators. We show $Q_{13}$ results as an example of the analysis,
Fig. 11 (d). In general the ensemble means are statistically indistinguishable
from the $\sim$9 month means, $\\{T_{1}\\}\sim\langle T_{1}\rangle_{T}$.
Similar results are observed for the remaining odd qubits.

Figure 11: (a) The $k$ = 5 systems (i.e., subsets) formed from the $Q_{13}$
time series. The $T_{1}(T_{i})$ measurements are re-indexed to form equal
sized subsets with a new maximum index, $m$, of $m_{max}$ = 61. An ensemble, a
set of 5 points in this case, is formed from the subsets for each $m$ time
index. The vertical line is a guide to the eye representing one ensemble
{$T_{1,k}(T_{m}$ = $20)$}. (b) Ensemble average dependence on time index for
the 5 subset and 32 subset partitions. (c) violin plot results of t-tests for
different system choices, $k_{max}$. More systems in the ensemble, $k_{max}$,
corresponds to shorter time series. The dotted line marks 0.05 as a guide to
the eye.

We note that we also investigated statistical tests of independence of the
subsets. The overall conclusions are not changed when subsets are rejected
from the ensemble based on failing statistical independence tests (i.e., the
non-parametric Wald-Wolfowitz runs testWald and Wolfowitz (1943)) as the
majority of subsets are found independent to the limits of the sensitivity of
the runs test.

##  Appendix G Ergodicity of ensemble averaging of $T_{1}(\omega,t_{i})$

We found in Appendix B that the $T_{1}$ time series of the qubits was
generally weakly stationary as well as behaving ergodically, Appendix F. The
time series $T_{1}(\omega_{q}+\omega_{j},t_{i})$ can also be expected to
behave as a weakly stationary ergodic time series as there is nothing uniquely
distinctive about the frequency $\omega_{q}$. We further explicitly note that
the sums and averages of neighboring time series will be weakly stationary and
ergodic.

For frequencies close to $\omega_{q}$, we may assume $\lim_{\omega_{j}\to
0}\langle T_{1}(\omega_{q}+\omega_{j})\rangle\rightarrow\langle
T_{1}\rangle_{T}$. We then ask whether an accurate $\langle T_{1}\rangle_{T}$
estimator may be formed from an ensemble of neighboring
$T_{1}(\omega_{q}+\omega_{j},t_{i})$ values. That is:

$\\{T_{1}(t_{i})\\}=\frac{1}{S}\sum_{j=0}^{S}T_{1}(\omega_{q}+(j\chi-\Delta\omega),t_{i})\simeq\langle
T_{1}(\omega_{q})\rangle_{T}$ (13)

where $S$ is the number of frequencies at which $T_{1}$ is sampled to form an
ensemble estimator; $\chi$ is the frequency spacing of sampling in a single
scan (see main text section V or Fig. 12 inset); and $\Delta\omega$ is the
maximum frequency range sampled, parameterized by $\chi$ as
$\Delta\omega=\frac{(S-1)}{2}\chi$ (see section V in main text).

To develop a heuristic for $\chi$, we choose a minimum frequency spacing
$\omega_{j}-\omega_{j+1}$ that produces approximately independent
$T_{1}(\omega_{q}+\omega_{j},t_{i})$s in the ensemble average of
$\\{T_{1}(t_{i})\\}$ (i.e., weak correlation with neighboring
$T_{1}(\omega_{q}+\omega_{j\pm 1},t_{i})$). We calculate the correlation
between frequencies to identify a $\chi$ that reduces the correlation below
$\sim$0.2. We show an illustrative example of the autocorrelation of Stark
shifted frequencies, $\langle
T_{1}(\omega,0)T_{1}(\omega+\omega_{lag},0)\rangle$), for $Q_{3}$, Fig. 12(a).
The autocorrelation shows a fall off to weak correlation over $\sim$1-2 MHz
centered around $\omega_{lag}\sim$ 0 for a single time index. Similar
magnitude fall offs were observed for all qubits examined. This suggests a
heuristic spacing of $\chi\simeq$ 1-2 MHz to obtain weakly correlated
$T_{1}(\omega_{q}+\omega_{j})$ with $T_{1}(\omega_{q})$.

Figure 12:  (a) Autocorrelation of $Q_{3}$ at $t$ = 0. The autocorrelation is
for negative Stark shift only. The autocorrelation is normalized and mean
value detrended. (inset) Schematic of the formation of an ensemble average,
$\\{T_{1}\\}$, from a set of $T_{1}(\omega_{q}+\omega_{j},t_{i})$ spaced by a
frequency $\chi$ to minimize correlation between the $T_{1}$ values and
sampling a range, $\Delta\omega$.

We empirically test the sensitivity of the error in the $\\{T_{1}\\}$
estimator to changing $\Delta\omega$ using $\chi$ = 2 MHz. We define
$\\{T_{1}\\}$ estimator error for each qubit as,
$\frac{\mid\\{T_{1}(t_{i})\\}-\langle T_{1}\rangle_{T}\mid}{\langle
T_{1}\rangle_{T}}$. The error dependence, averaged over the 75 time steps for
different $\Delta\omega$ (i.e., different $\chi$), is shown in Fig. 13 (a).
Notably the qubits have different dependencies on increasing $\Delta\omega$.
There is no single $\Delta\omega$ that is optimal for every qubit. The
differences between qubits helps explain the complex $\langle
R(\Delta\omega)\rangle_{t_{0\rightarrow n}}$ dependence shown in the main
text, Fig. 5 (c). The underlying cause for the differences is perhaps related
to the details of the local spectral diffusion for each qubit. The nearly 1:1
relationship is further supported by measurements on other IBM devices, see
Appendix K.

We now ask whether $\\{T_{1}\\}\simeq\langle T_{1}\rangle_{T}$. The 75
$\\{T_{1}\\}$ distributions for each qubit are compared to the $\langle
T_{1}\rangle_{T}$ time series distributions using a t-test. Many of the 75
time indexed $\\{T_{1}(t_{i})\\}$s are statistically indistinguishable but
many are not. The probability of $\\{T_{1}\\}$ being indistinguishable from
$\langle T_{1}\rangle_{T}$ is shown in Fig. 13 (b).

We illustrate the proximity of the estimator by showing the combined
distribution of the 75 means for $Q_{7}$, Fig. 13 (c) over the $\sim$270
hours. The means, $\\{T_{1}(t_{i})\\}$ and $\langle T_{1}\rangle_{T}$, are
within a few $\mu$s, less than a $\sigma$, but the null hypothesis (i.e.,
indistinguishable) is rejected because the distributions are sufficiently
different. It is possible that the overall $\\{T_{1}\\}$ distribution would
converge accumulating over a longer time. The $\sim$270 hours is insufficient
for all the $T_{1}(\omega_{q}+\omega_{j},t_{i})$ time series to unambiguously
converge around a mean value, see for example Fig. 14 (a) in Appendix H. Using
longer periods of time to measure $\langle
T_{1}(\omega_{q}+\omega_{j})\rangle$ could therefore potentially lead
$\\{T_{1}\\}$ to be more strictly ergodic.

Empirically, $\Delta\omega\leq$ 10 MHz produces $\\{T_{1}(t_{i})\\}$
distributions that are more likely to pass the t-test. The R values are also
slightly better ranging from 0.89 to 0.91, reflecting that regardless of
choice we find $\\{T_{1}\\}\simeq\langle T_{1}\rangle_{T}\pm\sigma$. For
illustration, we show the $\langle T_{1}\rangle_{T}$, $\\{T_{1}(t_{i})\\}$
pairs for $\Delta\omega$ = $\pm$6 MHz and $\chi$ = 2 MHz in Fig. 13 (d), which
has a Pearson R = 0.91. A very reliable 1:1 relationship is observed in many
other qubits measured on other IBM devices, see Appendix K.

Figure 13:  (a) Dependence of mean estimator error on $\Delta\omega$ using a
frequency spacing of $\chi$ = 2 MHz. The $\Delta\omega$ shown are $\pm$2 MHz
($\bullet$), $\pm$4 MHz ($\blacksquare$), $\pm$6 MHz ($\times$), $\pm$10 MHz
($\blacktriangle$) and $\pm$15 MHz ($\bigstar$). (b) Probability, out of 75
t-tests, that $\\{T_{1}\\}$ and $\langle T_{1}\rangle_{T}$ are found
indistinguishable as a function of range $\Delta\omega$ using a frequency
step, $\chi$ = 2 MHz. The $\Delta$s are $\pm$2 MHz ($\bullet$), $\pm$4 MHz
($\blacksquare$), $\pm$6 MHz ($\times$), $\pm$10 MHz ($\blacktriangle$) and
$\pm$15 MHz ($\bigstar$).(c) Histogram of $T_{1}$ measured at frequency steps
of 2 MHz over $\Delta\omega$ = $\pm$15 MHz for $Q_{7}$. Overlaid are the
$T_{1}(T_{i})$ measurements from the $\sim$9 month time series. Normal fits
for both distributions are shown in solid lines. (d) $\\{T_{1}\\}$ as a
function of $<T_{1}>_{T}$ for each qubit. One sigma standard deviations of the
distributions are shown as error bars for each $\\{T_{1}\\}$. $\Delta\omega$ =
$\pm$6 MHz and $\chi$ = 2 MHz.

To summarize, $\\{T_{1}\\}\simeq\langle T_{1}\rangle_{T}\pm\sigma$. More
statistically rigorous comparison by t-test indicates that
$\\{T_{1}(\Delta\omega,\chi)\\}$ is more quasi-ergodic than strictly ergodic
as it produces indistinguishable estimates of $\langle T_{1}\rangle_{T}$ for
less than 100% of the ensembles. The quasi-ergodic results were found for
$\chi$ and $\Delta\omega$ that were heuristically chosen and applied to all
qubits. Individually optimized $\chi$ and $\Delta\omega$ will reduce
disagreement and should become fully ergodic, certainly in the limit of
$\Delta\omega,\chi\rightarrow 0$ converging trivially on the single time
series $T_{1}(\omega_{q},t_{i})$. Optimal choices to achieve full ergodocity,
while minimizing the number of measurements (e.g., total time to obtain the
$T_{1}$ estimator), are left for future work. We speculate that this might
include forming an ensemble average with a physical model guided weighting of
the $j^{th}$ elements of $\\{T_{1}(\omega_{q}+\omega_{j})\\}$. For immediate
application of this approach, similar magnitude $\Delta\omega$ and $\chi$
values could be applied to other devices with the expectation that similar
magnitude R values will be obtained as the R value is not strongly sensitive
to the detailed choice of $\Delta\omega$ and $\chi$, see Appendix K.

##  Appendix H $R_{T_{1}}$ time dependence estimate from spectroscopy

We examine autocorrelation of $T_{1}$ from the data of Fig. 4 and 6 at
$\Delta\omega_{q}$ = 0 MHz. The $T_{1}$ time series (see Fig. 14(a)) for each
of the qubits is deduced assuming,

$T_{1}=\frac{-\tau}{\ln(P_{1}(\Delta\omega_{q}=0))},$ (14)

where $\tau=50\mu$s is the $T_{1}$ delay time, and $P_{1}(\Delta\omega_{q}=0)$
is the measured probability of being in the $|1\rangle$ state at time $\tau$
at the bare qubit frequency.

Figure 14: (a) $T_{1}$ time series calculated from $P_{1}$ spectroscopy for
$Q_{15}$ at zero frequency shift, (b) autocorrelation of the $Q_{15}$ time
series shown in (a) detrended with the mean and using padded zeros at the end
of the time series to maintain a constant number of elements for all negative
and positive lags (i.e., see online documentation for numpy.correlate). (c)
Pearson $R$ correlation between the odd qubits $T_{1}$ from spectroscopy and
$\langle T_{1}\rangle_{T}$. The autocorrelation is done for every time points
or intervals as indicated in the legend.

Autocorrelation of the time series in many of the qubits show decaying
correlation over the first 10-30 hours followed by weaker autocorrelation at
longer times consistent with the longer $T_{1}$ series from the $\sim$9 month
data set. This is seen in Fig. 14(b). The autocorrelation is detrended with
the mean and normalized to its estimated variance. Some qubits appear to have
longer term drifts, see appendix D.

We calculate the Pearson $R$ correlations to the odd qubit $\langle
T_{1}\rangle_{T}$ in Fig. 14 (c). The correlations are shown for time series
for different intervals of time available in the data set. This is to
highlight the effects of correlation. For example, if the interval in time was
doubled in an uncorrelated time series, it will double the time to achieve the
same $R$ value on average. However, if there are strong correlations,
increasing the interval times can decrease the time needed to achieve the same
$R$ value. Achieving $R$ $\sim$ 0.8 requires order of days or longer if the
best interval is unknown. We note that the dependence on interval does vary
some depending on what frequency is used for the $P_{1}$ data. The
autocorrelation hampers achieving strong correlation (i.e., $R$ $\sim$ 0.8) in
times shorter than 1-2 days.

##  Appendix I Time dependence of TLS linewidths

The time dependence of the TLS position in frequency is a topic of great
interest Müller, Cole, and Lisenfeld (2019); Klimov  _et al._ (2018); Schlör
_et al._ (2019); de Graaf  _et al._ (2020). The linewidth (i.e., the width of
the distribution of TLS frequency positions over time) is, for example,
suggested to depend on the volume density of thermal fluctuators (TF) Black
and Halperin ; Klimov  _et al._ (2018). TFs are low energy TLSs (i.e.,
$E_{TLS}\lesssim kT$). They are called thermal fluctuators (TF) because their
configurations dynamically change in time due to thermal excitation. This
results in a bath coupling to the higher energy TLSs (e.g.,
$E_{TLS}\sim\hbar\omega_{q}$) that produces the TLS spectral diffusion. A
linewidth time dependence is, therefore, expected to depend on the TF density
and coupling strength. An example of this dependence has been simulated in
previous work Klimov  _et al._ (2018). This appendix provides a reference and
comparison to previous TLS linewidth characterization. We present TLS
linewidths for $Q_{15}$.

The spectroscopy produces a discrete function of $P_{1}$ bins,
$P_{1}(\omega,\tau$,t). To track the time evolution of the TLS position, we
putatively associate each minima of a significant dip in $P_{1}$ with a TLS.
We find the position of each min($P_{1}(\omega,\tau$,t)) and record the
location for each time slice. We only record the TLS position if min($P_{1})$
is below a threshold, $P_{th}$, of 0.315 to remove spurious markers of TLS
location due to smaller fluctuations in $P_{1}$ (i.e., we focus on strongly
coupled TLSs). The threshold of $P_{\textrm{th}}=$0.315 corresponds to a
$T_{1}<43\mu$s. The threshold was chosen by visual inspection to best minimize
spurious points. The resulting TLS tracks are shown in Fig. 15. White
indicates a min($P_{1}(\omega,\tau$,t)). As mentioned earlier in the paper,
similar qualitative behavior is observed in the other qubits, see for example
spectroscopy of odd qubits in appendix A as well as in other works Müller,
Cole, and Lisenfeld (2019); Klimov  _et al._ (2018).

Figure 15: Peak locations of TLSs using a binary filtered replot of $Q_{15}$
with a $P_{1}$ = 0.315 threshold at a delay of $50\mu$s.

A linewidth can be crudely estimated by a cumulative histogram of the TLS
positions as a function of time. The resulting histogram for $Q_{15}$ is shown
in Fig. 16. Previous reports have fit the time dependence of spectral
diffusion linewidths with a normal distribution Klimov  _et al._ (2018);
Herzog and Hahn (1956) although we note that other functional forms are also
predicted Black and Halperin . Nevertheless, to obtain a simple quantitative
description of the linewidth and time dependence, we follow the
phenomenological model of a normally distributed peak position, fitting each
linewidth with a Gaussian and reporting the standard deviations, $\sigma$, in
Table 4. We report diffusivities $D_{1d}$ assuming a standard one-dimensional
model of the following form:

$C(t)=Ne^{{((\Delta\omega_{q}-\mu)^{2}}/{4\sigma(t)^{2}})}$ (15)

where $C$ is counts, $t$ is time, $\sigma$ is the standard deviation of the
distribution, $\mu$ is the center frequency, $N$ is a fitting constant and for
one-dimensional diffusion $\sigma(t)=\sqrt{2D_{1d}t}$. The spectral diffusion
visually appears to be of the same order of magnitude in $Q_{15}$ as the other
qubits so we present the values in Table 4 as representative of the order of
magnitude of Ds for the device. We also calculate $D_{K}$ with the units used
by the recent work of Klimov et al. Klimov  _et al._ (2018) (i.e.,
$\sigma(t)=2D_{K}t^{0.5}$) to provide rapid comparison to the extracted
diffusivity and modeling in that work. However, we note that while we
associate diffusivities to individual features that we track over time, the
estimate provided in ref. Klimov  _et al._ (2018) is a single value fit to the
consolidated linewidths of thirteen TLSs, an effective ensemble diffusivity.

Figure 16: (a) histogram of the TLS positions from Fig. 15. The distribution
of TLS positions are fit to Gaussians. The range over which the fit was done
is indicated by the truncation points of the Guassians.

Simulations of TLS spectral diffusion from Klimov et al. offer a suggestive
and appealing link between TF volume densities and $D_{K}$ Klimov  _et al._
(2018). Despite the differences in analysis, the $D_{K}$s from this work for
the most diffusive TLSs are over an order of magnitude smaller than reported
by Klimov et al., 2.5 MHz (hr)0.5, which might be interpreted as lower thermal
fluctuator densities. However, the authors emphasize that this is a dubious
conjecture. Doubts include: the model assumption of a time dependence of
$\propto\sqrt{t}$; the related model assumption of an unbounded random walk of
the TLS; and differences in the details of the temperature, between the two
experimental setups, which might lead to differing numbers of active
fluctuators. Different time dependencies might indeed be expected (e.g., for
different time regimes or dominant bath couplings Herzog and Hahn (1956);
Black and Halperin ; Witzel  _et al._ (2012); Madzik  _et al._ ). Perhaps more
importantly, the linewidths are likely truncated due to distance attenuated
coupling mechanisms in the bath Klauder and Anderson (1962) (e.g., dipole
coupling to thermal fluctuators). The assumption of ever increasing
$\sigma(t)$ can lead to significant disagreement (i.e., for cases of longer
time intervals of collection such as done in this work). We therefore do not
presently put substantial weight on the comparison of extracted diffusivities
until a more complete understanding of the linewidth time dependence is
established.

Position $(\textrm{MHz})$ |  $\sigma$ $(\textrm{MHz}$) |  $D_{K}$ $(\textrm{MHz~{}hr}^{-0.5}$) |  $D_{1d}$ ($\textrm{MHz${}^{2}$~{}hr}^{-1}$)
---|---|---|---
-19.3 | 0.65 | $2.4\times 10^{-2}$ | $1.1\times 10^{-3}$
-16.9 | 0.55 | $2.6\times 10^{-2}$ | $0.8\times 10^{-4}$
-7.9 | 1.06 | $3.9\times 10^{-2}$ | $3.1\times 10^{-3}$
-0.97 | 0.16 | $6.1\times 10^{-3}$ | $7.3\times 10^{-5}$
1.2 | 0.23 | $8.6\times 10^{-3}$ | $1.4\times 10^{-4}$
9.9 | 1.03 | $3.1\times 10^{-3}$ | $3.0\times 10^{-3}$
22.6 | 0.78 | $3.0\times 10^{-2}$ | $1.7\times 10^{-3}$
Table 4: Table 4 Estimated standard deviations for each of the TLS peak
distributions after 175 hours and diffusivities following ref. Klimov  _et
al._ (2018) $D_{K}$ or a one-dimensional diffusivity $D_{1d}$.

We also note two additional challenges to the accuracy of linewidth analysis
of individual TLSs, beyond the limits of validity of the one-dimensional
model, those being: (1) dips can overlap in ambiguous ways and (2) there can
be uncertainty in assignment of positions related to other TLS-like features
that migrate through the frequency region of interest. In Fig. 15 one can see
at least one other weaker dip that weaves between the two prominent ones in
the -15 to -20 MHz range, for example.

##  Appendix J Almaden measurement details

The Almaden device was a deployed system with cloud access. A daily
calibration was done, which included $T_{1}$ measurement. A database recorded
calibration measurements and the measurement times. Some additional
measurements were added to the database due to custom checks and recalibration
of qubits that were outside of the daily cals. The database was queried from
2019-09-13 to 2020-07-15.

The $T_{1}$ measurement was done for 41 time points logarithmically spaced up
to 500 $\mu$s using 300 shots per time point. A simple exponential fit was
made to the decay. The TLS spectroscopy was done using 501 frequency points
per direction of Stark shift with 1000 repeated rounds for each point. The
repetition time was 1 ms. This time can be substantially reduced with faster
reset of the initial state, for example.

##  Appendix K Supporting evidence for 1:1 correlation between $\langle
T_{1}\rangle_{\omega,t}$ and $\langle T_{1}\rangle_{T}$

We provide additional evidence that $\langle T_{1}\rangle_{\omega,t}$ is a 1:1
estimator of $\langle T_{1}\rangle_{T}$. We compile long time averages
$\langle T_{1}\rangle_{T}$ for 458 qubits and compare them to a single $T_{1}$
measurement in Figure 17 (a) to provide an illustrative example of statistical
spread and resulting R value. A 1:1 guide to the eye is overlaid and a Pearson
R of 0.72 is measured for the single $T_{1}$ estimator of $\langle
T_{1}\rangle_{T}$ for this instance. The length of the time series of daily
$T_{1}$ measurements depends on the amount of time the device was deployed.
The total time duration over which the $T_{1}$ measurements were done are
indicated in Table 5.

We contrast this with a $\langle T_{1}(t_{i},n=1)\rangle_{\omega,t}$ measured
from a single spectroscopy scan of each of the qubits, Fig. 17 (b) randomly
selected in the spectroscopy time series of measurements taken approximately
every 6 hours. For the spectroscopy measurements, the $\omega_{s}=\pm$80 MHz
and drive amplitude was swept to a fixed amplitude in both the negative and
positive shifts. This results in a total $\Delta\omega\sim$ 25 MHz. Each qubit
shifts slightly differently due to differences such as line attenuation. The R
value between $\langle T_{1}(t_{i})\rangle_{\omega}$ and $\langle
T_{1}\rangle_{T}$ was 0.82. Visually, a tighter concentration around a one-to-
one correspondence with $\langle T_{1}\rangle_{T}$ is observed from the
$\langle T_{1}\rangle_{\omega,t}$ estimator than relying on single $T_{1}$
measurements, consistent with observations made in the main text.

We also examine how $\langle T_{1}(n)\rangle_{\omega,t}$ converges with
$\langle T_{1}\rangle_{T}$ averaging over the duration of available $\sim$6
hour repeated spectroscopy measurements for each device. The time series
durations for the spectroscopy, $t_{max}$, are indicated for each device in
table 5. The R value improves to 0.91. The source of residual disagreement
likely comes in part from the lack of custom optimization of $\Delta\omega$
and best choice of weighting of $T_{1}(\omega_{q}+\omega_{j},t_{i})$.

Figure 17: Scatter plots of $\langle T_{1}\rangle_{T}$ compared to different estimators for 458 qubits from devices listed in table 5. (a) Comparison of $\langle T_{1}\rangle_{T}$ to a single $T_{1}(T_{i})$ measurement for each qubit randomly selected from the $T_{1}$ time series. (b) $\langle T_{1}(t_{i})\rangle_{\omega}$ estimated from $P_{1}(\omega_{q}+\omega_{j},t_{i})$ measurements for a spectroscopy scan at a randomly selected time, $t_{i}$, in the spectroscopy time series. (c) $\langle T_{1}\rangle_{\omega,t}$ averaging over all scans available, see the duration of the time series for each device in table 5. The Pearson Rs for the three cases are 0.72, 0.82 and 0.91, respectively Device |  Qubits |  $T_{max}$ (days) |  $t_{max}$ (days)
---|---|---|---
Santiago | 5 | 453 | 169
Bogota | 5 | 437 | 169
Montreal | 27 | 452 | 220
Toronto | 27 | 449 | 169
Belem | 5 | 220 | 192
Brooklyn | 65 | 145 | 108
Casablana | 7 | 362 | 169
Dublin | 27 | 350 | 168
Kolkata | 27 | 236 | 180
Lima | 5 | 216 | 192
Quito | 5 | 222 | 192
Sydney | 27 | 366 | 168
Manila | 5 | 123 | 101
Jakarta | 7 | 138 | 118
Hanoi | 27 | 132 | 107
Bangkok | 27 | 130 | 78
Lagos | 7 | 68 | 45
Nairobi | 7 | 76 | 55
Mumbai | 27 | 278 | 198
Cairo | 27 | 55 | 30
Algiers | 27 | 73 | 51
Manhattan | 65 | 402 | 191
Table 5:  Device names, number of qubits and the length of time series for the
$T_{1}$ measurements, $T_{max}$, and Stark spectroscopy, $t_{max}$. The
$T_{1}$ measurements are taken every 24 hours and the spectroscopy
measurements are taken approximately every 6 hours.