# Stability of Quantum Computers

Samudra Dasgupta

###### keywords:

Quantum Computing, Probabilistic error cancellation, Non-stationary quantum
channels, Bayesian statistics, Stability, Reproducibility, Reliability,
Accuracy

2024 May Dr. Travis Humble Dixie L. Thompson Data Science and Engineering
Doctor of Philosophy Engineering Bredesen Center for Interdisciplinary
Research The University of Tennessee, Knoxville 4 Dr. Travis Humble Dr.
Russell Lee Zaretzki Dr. Himanshu Thapliyal Dr. Rebekah Herrman

##

Dedicated to my parents.

## Acknowledgments

I am deeply thankful to Travis for his invaluable guidance, support, and
expertise during my entire PhD journey. His steadfast dedication and
insightful mentorship has significantly influenced my research and personal
development. I also want to express my gratitude for the opportunity to
utilize the resources of the Oak Ridge Leadership Computing Facility (a United
States Department of Energy, Office of Science User Facility).

## Abstract

Quantum computing’s potential is immense, promising super-polynomial
reductions in execution time, energy use, and memory requirements compared to
classical computers. This technology has the power to revolutionize scientific
applications such as simulating many-body quantum systems for molecular
structure understanding, factorization of large integers, enhance machine
learning, and in the process, disrupt industries like telecommunications,
material science, pharmaceuticals and artificial intelligence. However,
quantum computing’s potential is curtailed by noise, further complicated by
non-stationary noise parameter distributions across time and qubits. This
dissertation focuses on the persistent issue of noise in quantum computing,
particularly non-stationarity of noise parameters in transmon processors. It
establishes a framework comprising computational accuracy, device reliability,
outcome stability, and result reproducibility for assessing noisy outcomes
amidst time-varying quantum noise. It further aims to determine the upper and
lower bounds for this framework using available noise characterization data,
in terms of the distance between time-varying noise densities. Using real data
from a transmon processor, it validates the bounds on a test quantum circuit.
It also demonstrates that if the physical platform’s noise stays within the
bounds determined by the analysis, experimental reproducibility can be
guaranteed with a high degree of certainty. Furthermore, it develops a
Bayesian algorithm to enhance outcome stability and accuracy for probabilistic
error cancellation (PEC) in presence of time-varying quantum noise. The
results obtained from experiments using a 5-qubit implementation of the
Bernstein-Vazirani algorithm conducted on the ibm_kolkata device, underscore
the effectiveness of the adaptive algorithm, showing a 42% improvement in
accuracy over non-adaptive PEC and a 60% improvement in stability. Considering
the time-varying stochastic nature of quantum noise, integrating adaptive
estimation in error mitigation is crucial. In summary, by delving into the
complexities of non-stationary noise in quantum computing, this dissertation
provides valuable insights into quantifying and enhancing stability of
outcomes from noisy quantum computers.

## Table of Contents

toc

###### List of Tables

  1. 1 *
  2. 2 *
  3. 3.1 Moment-based distance by Distribution
  4. 3.2 Transition Probabilities for Uncorrelated Noise
  5. 4.1 The 16-parameter model derived from the ibm_washington data set has four types of quantum noise processes: (i) ‘ABC’: asymmetric binary channel model, (ii) ‘CP’: coherent phase error model, (iii) ‘DP’: depolarizing noise model and (iv) ‘TR’: thermal relaxation noise model. Note that the two-qubit model ‘DP$\otimes$DP’ is a tensor product of depolarizing noise.
  6. 4.2 Hellinger distance values for the device parameters.
  7. 4.3 Gate lookup table
  8. 5.1 Register values for $\gamma_{D}(\tau)$ and $\gamma_{\textrm{max}}$
  9. 5.2 Device parameters
  10. 6.1 True means of the time-varying noise Pauli channel coefficients
  11. 6.2 Estimated Pauli coefficients for peiod 1
  12. 6.3 Estimated coefficients for period 2

###### List of Figures

  1. 1.1  Impact of the non-stationary noise, leading to irreproducible outcomes. The results show histograms (after state preparation and measurement (SPAM) noise mitigation) upon executing the Bernstein-Vazirani circuit on ibm_mumbai.
  2. 1.2 Evidence of spatial non-stationarity in the mean values for quibt decoherence times $T_{1}$ and $T_{2}$ for the 127-qubits of the ibm_washington device, generated on 14 Jan 2023 10:20 PM UTC.
  3. 1.3 state preparation and measurement (SPAM) fidelity distributions on ibm_toronto for qubits $0-8$ as measured on 8 April 8 2021, between 8:00-10:00pm (UTC-05:00).
  4. 1.4  state preparation and measurement (SPAM) fidelity distributions on ibm_toronto for qubits $9-17$ as measured on 8 April 8 2021, between 8:00-10:00pm (UTC-05:00).
  5. 1.5  state preparation and measurement (SPAM) fidelity distributions on ibm_toronto for qubits $18-26$ as measured on 8 April 8 2021, between 8:00-10:00pm (UTC-05:00).
  6. 2.1 Interleaved sub-circuit post transpilation.
  7. 2.2  Spatial non-stationarity of decoherence times. Indivudal error-bars do not capture the variation across qubits. (a) In terms of the $T_{1}$ parameter, qubit 15 performed the best, exhibiting a median $T_{1}$ time of 184 microseconds, while qubit 4 performed the worst with a median $T_{1}$ time of 74 microseconds. Across all 27 qubits, the median $T_{1}$ time was 116 microseconds, with a standard deviation of 1.5 microseconds, and the range of $T_{1}$ times spanned from 32 to 297 microseconds. (b) For the $T_{2}$ parameter, qubit 3 demonstrated the highest performance, displaying a median $T_{2}$ time of 76 microseconds, while qubit 19 had the poorest performance with a median $T_{2}$ time of 15 microseconds. Across all 27 qubits, the median $T_{2}$ time averaged 30 microseconds, with a standard deviation of 1.1 microseconds, and the range of $T_{2}$ times varied from 3 to 191 microseconds.
  8. 2.3 This figure presents non-stationary temporal dynamics of decoherence for a qubit on IBM’s belem device. The top figure shows $T_{1}$ relaxation time series for qubit 0, where two datasets were collected for 5 ms each on the same day, separated by a vertical line. The blue dataset varies between 116-126$\mu$s with a mean of 122 $\mu$s and a standard deviation of 2 $\mu$s, while the green dataset varies between 104-114$\mu$s with a mean of 108 $\mu$s and a standard deviation of 2 $\mu$s. The bottom figure displays Ramsey dephasing time ($T_{2}$) series for qubit 0, where two datasets were collected for 5 ms each on the same day, separated by a vertical line. The blue dataset varies between 67-73$\mu$s with a mean of 70 $\mu$s and a standard deviation of 1 $\mu$s, while the green dataset varies between 60-64$\mu$s with a mean of 62 $\mu$s and a standard deviation of 1 $\mu$s, collected around different times on the same day. The data shows significant non-stationarity in decoherence values over a 30-minute interval.
  9. 2.4  (a) Schematic of the 27-qubit device ibm_kolkata. (b) Estimated $T_{1}$ and $T_{2}$ time-series for qubit 0 as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023.
  10. 2.5  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 1 and (b) qubit 2.
  11. 2.6  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 3 and (b) qubit 4.
  12. 2.7  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 5 and (b) qubit 6.
  13. 2.8  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 7 and (b) qubit 8.
  14. 2.9  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 9 and (b) qubit 10.
  15. 2.10  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 11 and (b) qubit 12.
  16. 2.11  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 13 and (b) qubit 14.
  17. 2.12  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 15 and (b) qubit 16.
  18. 2.13  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 17 and (b) qubit 18.
  19. 2.14  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 19 and (b) qubit 20.
  20. 2.15  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 21 and (b) qubit 22.
  21. 2.16  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 23 and (b) qubit 24.
  22. 2.17  Estimated $T_{1}$ and $T_{2}$ time-series as collected between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for (a) qubit 25 and (b) qubit 26.
  23. 3.1 A classical Markov model for the correlated error process.
  24. 3.2 Simulation result for addressability $F_{A}$ showing sensitivity to intra-register correlation.
  25. 3.3 Addressability: pairwise comparison when correlation parameter $u=0.12$. Ideally, there should have been zero spatial variation (as idealized simulations do not differentiate between qubits at different physical locations). However, small fluctuations are seen about a mean value of 95.80 with a standard deviation of 0.02. This arises due to the readout error fluctuations in the ibm_yorktown noise model.
  26. 3.4 The Bernstein-Vazirani circuit for a 4-bit secret string $r=1000$. $H$ represents the Hadamard gate and $Z$ represents the $Z$-gate. The meter symbols are measurement operations that project to the computational basis states. The measurement results are recorded in a classical register $c$.
  27. 3.5 Contribution to moment based distance ($d$) from increasing moment orders. The graph shows the results of comparing two normal distributions: $\mathcal{N}_{1}(\mu=\mu_{0},\sigma=\sigma_{0})$ and $\mathcal{N}_{2}(\mu=2\mu_{0},\sigma=2\sigma_{0})$ where $\mu_{0}=40$ and $\sigma_{0}=4$.
  28. 3.6 Simulated signal-to-noise ratios of the moment-based distance and total variation distance for two normal distributions of varying width.
  29. 3.7 When the two distributions are similar, then we expect each $S_{m}$ to be zero. Empirically, that happens as we increase the sample size. The lower order moments take longer to go to zero.
  30. 4.1 Schematic layout of the ibm_yorktown device produced by IBM. Circles denote register elements and edges denote connectivity of 2-qubit operations.
  31. 4.2 Schematic of the ibm_toronto device, produced by IBM. Circles represent register elements, while edges denote the connectivity for performing two qubit operations.
  32. 4.3 Experimentally observed probability density for SPAM (state preparation and measurement) fidelity of one of the register elements of the IBM transmon device named ibm_washington.
  33. 4.4  (a)-(e) Temporal stability of the SPAM fidelity $F_{\text{SPAM}}$ of each register element in the ibm_yorktown device. The top panel shows the average $F_{\text{SPAM}}$ of the register with associated variance, and the bottom panel shows a running calculation of the Hellinger distance using a one-month window. The dashed red line is the median value.
  34. 4.5 Spatial stability of the SPAM fidelity for the ibm_yorktown device from May 2019 to December 2020, where the inset highlights the registers with the maximum distance.
  35. 4.6  (a)-(e) Temporal stability of the gate fidelity $F_{G}$ for the CNOT gate for sequential register pairs in the ibm_yorktown device from March 2019 to December 2020. The top panel shows the average $F_{G}$ of the register pair and the bottom panel shows a running calculation of the Hellinger distance with respect to May 2019. The dashed red line is the median value.
  36. 4.7  (a)-(d) Temporal stability of the CNOT duty cycle for sequential register pairs in the yorktown device. The top panel shows the harmonic mean of the register decoherence time $T_{2}$ for the elements, the upper-middle panel shows the gate duration $T_{G}$, the lower-middle panel plots the corresponding duty cycle $\tau$, and the bottom panel presents the Hellinger distance for the duty cycle averaged over a one-month window. The dashed red line is the median value.
  37. 4.8 The spatial stability of the duty cycle $\tau$ for ibm_yorktown. The inset shows the experimental histograms for register pairs (0,1) and (2,3) which are separated by the largest Helligner distance of 0.789.
  38. 4.9 Spatial stability of the gate fidelity $F_{G}$ for the CNOT gates of the ibm_yorktown device from March 2019 to December 2020. The inset shows the distribution of gate fidelities for pairs (1,2) and (3,2), which yield a Hellinger distance 0.467.
  39. 4.10 The figure presents results of reliability testing on a transmon platform over 16 months. Plots of the three Hellinger distance measures, with the unmodified measure $H$ (blue line) being insensitive due to dimensionality issues. The blue and orange lines both capture the correlation structure of the joint distribution (with the orange line being normalized to enhance discrimination power), while the green line lacks correlation capture. The normalized measure $H_{\textrm{normalized}}$ ranges between 0.41 and 0.92, while the average measure $H_{\textrm{avg}}$ varies between 0.431 and 0.51. The latter captures monthly variations in the marginal Hellinger distance for each of the 16 error parameters but fails to account for correlations.
  40. 4.11  Decomposition of the sources of quantum noise non-stationarity. The degree of non-stationarity varies amongst the sources. The plot shows the contributions made by each noise type to the composite Hellinger distance (a measure of the degree of non-stationarity). The Hellinger distance measures the statistical distance between the joint distribution of the noise observed in Apr-2023 to the joint distribution of the noise observed in Jan-2022. The various noise types contribute varying percentages but no single term dominates the sum.
  41. 4.12  The Pearson correlation coefficients between the 16 characterization parameters as observed in April-2023. Dark blue and dark red colors represent the Pearson coefficients $1$ and $-1$, respectively. The axes labels $(x0,x1,\cdots x15)$ correspond to various noise sources as listed in Table 5.2.
  42. 4.13 Normalized mutual information of register pairs in the toronto device sampled 11:00-11:30 PM (UTC-5) on 11 December 2020. Data corresponds to register pairs prepared in the Bell state and the inset shows the range of the lowest and highest values for the normalized mutual information.
  43. 4.14 Experimentally observed duty cycle ($\tau$) for the CNOT gate for washington.
  44. 4.15 The 95% temporal confidence interval (vertical orange lines) for CNOT duty cycle ($\tau$) for the physical nearest-neighbor connections of ibm_washington. The data-set contains the values for all the 144 physical CNOT gates of ibm_washington between 1-Dec-2021 to 31-May-2022. The dashed black line represents the mean.
  45. 4.16 Temporal reliability of the SPAM fidelity $F_{\text{SPAM}}$ of register element $q=37$
  46. 4.17 Temporal reliability of the gate fidelity $F_{G}$ for the CNOT gate for register pairs $0$ and $14$.
  47. 4.18 Spatial reliability of the gate fidelity $F_{G}$ for the CNOT gates of the ibm_washington device. The heat map shows the Hellinger distance between the distributions for the gate fidelities of nearest-neighbor connections. Only the lower triangular matrix is shown to avoid redundancy. The upper triangular matrix as well as any data gaps in the lower triangular matrix are colored white. The inset shows the distribution of gate fidelities for pairs $(11,12)$ and $(19,20)$, which yield a Hellinger distance 0.99. This metric captures the probability that the quantum processor register is spatially dissimilar. The estimation of the distribution for CNOT-$i$ utilizes data for Sep-2022.
  48. 4.19 Spatial reliability of the duty cycle $\tau$ for the CNOT gates of the ibm_washington device. The heat map shows the Hellinger distance between the distributions for the duty cycles of nearest-neighbor connections. The inset shows the distribution of duty cycle for pairs (46, 47) and (96, 109), which yield a Hellinger distance 0.99. This metric captures the probability that the quantum processor register is spatially dissimilar. The estimation of the distribution for CNOT-$i$ utilizes duty-cycle data for Sep-2022.
  49. 4.20 Spatial reliability of the SPAM fidelity $F_{\text{SPAM}}$ for the ibm_washington device. The heat map shows the Hellinger distance between the distributions for each register pair. This metric captures the probability that the quantum processor register is spatially dissimilar. The inset shows the distributions of $F_{\text{SPAM}}$ for registers 12 and 37, which represent distance=0.99 due to minimal overlap. The estimation of the distribution for qubit-$i$ utilizes data for Sep-2022.
  50. 4.21  Spatial reliability of the gate fidelity $F_{G}$ for the CNOT gates of the ibm_washington device. The heat map shows the Hellinger distance between the distributions for the gate fidelities of nearest-neighbor connections. Only the lower triangular matrix is shown to avoid redundancy. The upper triangular matrix as well as any data gaps in the lower triangular matrix are colored white. The inset shows the distribution of gate fidelities for pairs $(11,12)$ and $(19,20)$, which yield a Hellinger distance 0.99. This metric captures the probability that the quantum processor register is spatially dissimilar. The estimation of the distribution for CNOT-$i$ utilizes data for Sep-2022.
  51. 4.22 Addressability of register pairs in the ibm_toronto device sampled 08:00-08:30 AM (UTC-5) on 11 December 2020. This data corresponds to the register prepared in the separable fiducial state. The inset shows the range i.e. the lowest and highest values for addressability. The average of (23,21) is the lowest value at 0.887 while all other values lie in the range [0.992, 1). The outlier is the only value that does not appear in the plot.
  52. 5.1 Circuit used for our experiment. In this figure, $H$ represents the Hadamard gate. The meter symbol denotes measurement gate.
  53. 5.2 Plot illustrating significant spatial non-stationarity in the register-wise variation of the SPAM asymmetry for ibm_toronto.
  54. 5.3 The plot depicts distance variation across the $27$ register elements of ibm_toronto device, illustrating the impact of spatially non-stationary noise on program outcomes and the dependence of output on register choice.
  55. 5.4 The plot illustrates the register-wise Hadamard gate angle error (in degrees) within ibm_toronto’s quantum register of $27$ elements, emphasizing spatial non-stationarity and cautioning against averaging qubit values for coherent noise mitigation.
  56. 5.5 This plot shows the register-wise variation of the composite accuracy metric $\gamma_{D}$ for ibm_toronto device with 27 qubits, where a higher value significantly impacts program precision. The graph highlights the necessity of re-estimating the metric due to temporal non-stationary noise in unreliable devices, emphasizing the crucial role of analyzing noise parameter interactions for desired accuracy.
  57. 5.6 Characterizing circuit accuracy on toronto. Plot of $\gamma_{\textrm{max}}$ (dashed line) and $\gamma_{D}$ (blue dots) for ibm_toronto on 8 April 2021. The blue dots are experimentally-observed data plotted using the characterization data versus the actual observed distance $(d)$ for each register element. Only a subset of qubits are shown.
  58. 5.7 Simulation demonstrating that when $H\leq H_{\textrm{max}}$ (i.e. a reliable, slowly varying noise platform), then $s\leq s_{\textrm{max}}$ (i.e. the temporal difference of the observable stays within the predicted upper bound).
  59. 5.8 A quantum circuit implementation of the Bernstein-Vazirani algorithm that employs 5 qubits, denoted $q_{0}$ to $q_{4}$. The first four qubits are used to compute the 4-bit secret string, while the fifth qubit serves as an ancilla and initially resides in the $\ket{-}$ superposition state. The symbol $H$ denotes the Hadamard gate while the oracle unitary ($U_{r}$) implements the secret string ($r$). The depolarizing noise channel is denoted by $\mathcal{E}_{\textrm{x}}(\cdot)$. A quantum measurement operation is represented by the meter box symbol at the circuit’s end.
  60. 5.9 The ratio $s/s_{\textrm{max}}$ for a simulated time-varying noisy circuit is plotted with respect to the increasing noise variance across 15 months. The results from noisy simulations of the Bernstein-Vazirani circuit with a secret string of 4, 8, or 12 bits validate the bound expected.
  61. 5.10 Contour plots to compare the probability densities of a two-dimensional subset of Hadamard gate errors for qubit 0 and 3 in Apr-2023, (a) with and (b) without correlation modeling using a copula function.
  62. 5.11  Simulations of the stability ratio $s/s_{max}$ times 100 for a 4-qubit Bernstein-Vazirani circuit using the noise characterization from the ibm_washington platform. The box-and-whisker plot of the monthly statistics are based on noisy circuit simulations using the joint probability distribution derived from data from 1-Jan-2022 to 30-Apr-2023. Ratio values below unity confirm that the upper bound is never exceeded.
  63. 5.12 Schematic layout of the $127$-qubit washington device produced by IBM. Circles denote register elements and edges denote connectivity of 2-qubit operations. The register elements $0$, $1$, $2$, $3$, and $4$ are mapped to the physical qubits $4$, $3$, $2$, $1$, and $0$, respectively, in the diagram above. The CNOT gates used in the circuit connect the physical qubits $(0,1)$ and $(2,1)$, in the diagram above, where the first number represents the control qubit and the second one represents the target qubit.
  64. 6.1  Time-varying density of SPAM fidelity. Data shown for register element #26.
  65. 6.2  SPAM fidelity time-series for qubit #26 for Dec-May 2022.
  66. 6.3 We simulated SPAM and gate error channels for a quantum circuit with 4 qubits that creates a uniform superposition across all the computational basis states using Hadamard gates. The blue bars represent the probability distribution across the measurement outcomes for the ideal, noiseless circuit while the orange bars represent the same for a realization of an execution on an unstable device.
  67. 6.4 The plot shows improved accuracy achieved using Bayesian optimization in presence of time-varying noise. Mitigation using average parameters deteriorates accuracy. The number of samples used in the Bayesian estimation process is held constant at $10^{4}$.
  68. 6.5 The 4-qubit quantum circuit used for the simulation experiment. Each qubit is assumed to have a different, independent SPAM error process. Each Hadamard gate is similarly assumed to have a different, independent gate error process.
  69. 6.6  The figure displays the noisy basis circuits for the linear combination step (the first step) in PEC, where each figure represents the noisy operations $\mathcal{E}_{\textrm{x}}\circ\mathcal{P}_{\sigma,\sigma^{\prime}}\circ\mathcal{G}$ and $G=H\otimes H$ is the desired, ideal operation. Note that different basis set choices are possible, depending on the hardware.
  70. 6.7 We model the degradation of a non-stationary Pauli noise channel by assuming the coherence time steadily decreases over time. Using a non-stationary Dirichlet distribution, we model the joint distribution of coefficients for a two-qubit circuit, which fluctuate as coherence times deteriorate. The y-axis represents the degree of non-stationarity, and the x-axis shows four time-periods. The Hellinger distance between the Dirichlet distributions at time $\tau$=0 and a later time is a measure of non-stationarity, increasing from 0 to 57%. This model is based on transmon platforms and is used as an experimental setup for our simulation experiments.
  71. 6.8 This graph compares adaptive and non-adaptive PEC for the $H\otimes H$ gate under time-varying noise. The y-axis shows the probability of observing a basis state. The black bars represent the ideal histogram for the test input, while red/orange bars are non-adaptive PEC results. Adaptive PEC improves accuracy. For example, the $\ket{00}$ probability increases from 57% to 72%.
  72. 6.9 This graph compares adaptive and non-adaptive PEC implementations over four time-periods with time-varying noise, using Hellinger distance to measure the difference from the ideal distribution. The black and orange bars represent the observed distributions with non-adaptive and adaptive PEC, respectively. Adaptive PEC significantly outperforms non-adaptive PEC in reducing Hellinger distance to 1.1%, and 3.2% compared to 7%, and 15% for non-adaptive PEC across time-periods. Adaptive PEC uses adaptive estimation of noise super-operators to improve accuracy compared to non-adaptive PEC. Time-varying noise underscores the need for adaptive PEC.
  73. 6.10 The likelihood function used in the Bayesian inference procedure for the second period.
  74. 6.11 The cost function used in the Bayesian inference procedure for the second period.
  75. 7.1 Non-stationary distribution functions of the state preparation and measurement (SPAM) fidelity for qubit 2 on ibm_kolkata superconducting device collected on Jan 15, 2024.
  76. 7.2 A 5-qubit implementation of the Bernstein-Vazirani algorithm with secret bit-string $r$.
  77. 7.3 Circuit diagram for implementing one of the 16 noisy super-operators for CNOT noise mitigation using PEC. The initial noisy CNOT gate is followed up by $\mathds{Y}$ and $\mathds{Z}$ gates on the control and target qubits respectively.
  78. 7.4 One of the 512 noisy basis circuits for mitigation using PEC in the 5-qubit implementation of the Bernstein-Vazirani circuit. The initial noisy CNOT gate is followed up with $\mathds{Y}$ and $\mathds{Z}$ gates on the control and target qubits respectively. The readout lines for qubit $0,1,2,3,4$, and $5$ are subjected to the Pauli gates $\mathds{X},\mathds{I},\mathds{X},\mathds{X},\mathds{I}$ respetively, prior to measurement, in this specific noisy basis circuit for PEC.
  79. 7.5  Comparison of the adaptive PEC algorithm with a non-adaptive approach for a density matrix simulation that implements a noisy 5-qubit circuit for solving the Bernstein Vazirani problem. (a) Average accuracy across the 5-qubits and (b) Average stability across the 5-qubits.
  80. 7.6 Qubit layout of the 27-qubit superconducting device ibm_kolkata. Circles represent superconducting transmons and lines indicate possible gate operations between sites.
  81. 7.7  This figure depicts the non-stationary noise on the experimental device ibm_kolkata. In plot (a), the blue line represents the depolarizing parameter for the target qubit, while the black line denotes the depolarizing parameter for the control qubit in the CNOT gate. Plot (b) illustrates five lines, each indicating the SPAM fidelity for the register elements. The x-axis corresponds to intra-calibration timestamps for January 15. The shaded regions denote the time-varying standard deviations.
  82. 7.8  These plots depict the non-stationary nature of the quasi-probability distribution. To maintain clarity, we opted not to plot all 512 bins of the histogram in a single plot. In (a), we observe the time-varying weight $\mathds{Q}_{0}$ for the raw Bernstein-Vazirani circuit which carries the most substantial weight, accounting for almost 80% of the distribution. In (b), we show the values of the quasi-probability bins for the subsequent 10 basis circuits, which are more than 10 times lower in magnitude compared to the first circuit. The non-stationary nature of the quasi-probability distribution becomes crucial given the lengthy data collection process required for PEC mitigation because the noise estimation becomes inaccurate in these time-frames. Our experiment took approximately 2 hours for each dataset comprising 512 circuits.
  83. 7.9  The figure comprises two plots, each depicting four graphs: (1) "No mitigation" presents raw Bernstein-Vazirani metrics without error mitigation, (2) "ROEM" shows metrics after readout error mitigation, with constant SPAM noise parameters, (3) "Non-adaptive PEC" displays metrics for PEC, and (4) "Adaptive PEC" exhibits metrics for PEC with adaptive noise mitigation. The x-axis denotes intra-calibration time-stamps (UTC) for Jan 15. Plot (a) illustrates the time-varying accuracy metric from Eqn. 7.16. It demonstrates the adaptive PEC’s 42% accuracy improvement over non-adaptive PEC. Plot (b) shows the time-varying stability metric from Eqn. 7.19. It illustrates the adaptive PEC’s 60% stability enhancement compared to non-adaptive PEC. These plots underscore the significant impact of non-stationary noise on PEC resilience. Due to PEC’s lengthy completion time (a couple of hours), adaptive methods are able to handle non-stationary noise conditions better.
  84. 7.10 Non-stationary distribution function of SPAM fidelity for qubit 0 of ibm_kolkata as observed on Jan 15, 2024.
  85. 7.11 Non-stationary distribution function of SPAM fidelity for qubit 1 of ibm_kolkata as observed on Jan 15, 2024.
  86. 7.12 Time-varying density for probability of observing ‘00’ on qubits 3 and 4 of ibm_kolkata, for the Bernstein-Vazirani circuit, as observed on Jan 15, 2024.
  87. 7.13 Time-varying density for probability of observing ‘01’ on qubits 3 and 4 of ibm_kolkata, for the Bernstein-Vazirani circuit, as observed on Jan 15, 2024.
  88. 7.14 Time-varying density for probability of observing ‘10’ on qubits 3 and 4 of ibm_kolkata, for the Bernstein-Vazirani circuit, as observed on Jan 15, 2024.
  89. 7.15 Time-varying density for probability of observing ‘11’ on qubits 3 and 4 of ibm_kolkata, for the Bernstein-Vazirani circuit, as observed on Jan 15, 2024.

## List of abbreviations

Table 1: * NSN |  Non-stationary noise
---|---
PEC |  probabilistic error cancellation
BC |  Bhattacharyya coefficient
BV |  Bernstein-Vazirani
SPAM |  State preparation and measurement
Tr |  Trace
WSS |  Wide-sense stationary
MCMC |  Markov chain Monte-Carlo

## Notations

Table 2: * $\ket{\psi}$ |  a pure quantum state
---|---
$\otimes$ |  tensor operator
Tr$(\cdot)$, Pr$(\cdot)$ |  Trace and probability operators respectively
$\rho$ |  density matrix
$\Omega,\lambda$ |  Hamiltonian and its eigenvalue
$T_{1},T_{2}$ |  Qubit relaxation and dephasing time respectively
$N_{s}$ |  Number of shots per circuit execution
$L$ |  Number of times a circuit is executed for collecting statistics
$n$ |  quantum register size
$\ket{v},\ket{v_{n-1}\cdots v_{0}},v_{i},v$ |  $\ket{v}$ is the short form for $\ket{v_{n-1}\cdots v_{0}}$. The latter denotes an n-bit string with each $v_{i}\in\\{0,1\\}$. Lastly, $v$ denotes the decimal form for the bit-string e.g. if $\ket{v}=\ket{0011}$, then $v=3$.
$\mathcal{U},\tilde{\mathcal{U}}$ |  a unitary matrix and its noisy version (could be non-unitary)
$R_{Y}(\theta)$ |  operator for rotation by an angle ($\theta$) about the Y-axis on the Bloch sphere
$\mathds{I},\mathds{X},\mathds{Y},\mathds{Z}$ |  2x2 identity and the Pauli X, Y, Z matrices
$\mathds{H}$, CNOT, SX, RZ |  Various quantum gates (Hadamard gate, Controlled-NOT, Square root of NOT and Unitary gate for rotation about Z-axis
$T_{D}$ |  delay gate duration
$\tilde{\mathds{H}}$ |  noisy Hadamard gate
$\mathcal{I}(X,Y)$ |  Mutual Information between two random variables $X$ and $Y$
$\mathcal{H}(X)$ |  Entropy of $X$
$\mathcal{O}$ |  A quantum observable
$\mathcal{E}_{\textrm{x}},\mathcal{E}_{\textrm{x}_{i}}(\rho)$ |  Multi-qubit quantum noise channel parameterized by the vector x; $\mathcal{E}_{\textrm{x}_{i}}(\rho)$ is a single-qubit quantum noise channel
$M_{k}$ |  Kraus operators characterizing the error channel $\mathcal{E}(\cdot)$
$F_{\text{SPAM}},\epsilon_{\text{SPAM}}$ |  SPAM Fidelity and SPAM error respectively
$F_{\text{G}},e_{G}$ |  Gate Fidelity and error per Clifford gate respectively
$T_{G}$ |  CNOT Gate Length
$\tau_{G}$ |  Duty Cycle
$\text{Beta}(\alpha,\beta)$ |  Beta function $(=\int_{0}^{1}t^{\alpha-1}(1-t)^{\beta-1}dt)$
$\Gamma(z)$ |  Gamma function $(=\int_{0}^{\infty}t^{z-1}\exp^{-t}dt)$
$\mathds{R},\mathds{C}$ |  Set of real and complex numbers respectively
$P_{2}(\mathds{R})$ |  Set of real polynomials of degree not more than 2
$\textrm{x},\textrm{x}_{i}$ |  $x=(\textrm{x}_{1},\cdots,\textrm{x}_{d})$ is a vector of parameters characterizing the noise of the quantum circuit
$X,X_{i}$ |  $X$ is a vector of random variables (stochastic device parameters) corresponding to x above; $X_{i}$ is a scalar random variable
$f_{X}(\textrm{x};t)$ |  multi-variate continuous probability distribution function for x; same as $f_{X}(\textrm{x};t)$
$F_{X}(\textrm{x};t)$ |  multi-variate continuous cumulative distribution function for x; same as $F_{X}(\textrm{x};t)$
$\eta(X,Y)$ |  Normalized Mutual Information between two random variables $X$ and $Y$
$F_{\text{A}}$ |  Addressability
$BC$ |  Bhattacharyya coefficient
$H_{X}(t_{1},t_{2})$ |  Hellinger distance between the multi-variate joint distributions of $X$ at times $t_{1}$ and $t_{2}$; if $t_{1}=0$ then we sometimes simply write it as $H_{X}(t)$
$H_{\text{avg}}(t_{1},t_{2})$ |  average of the Hellinger distances over the $d$ univariate marginal distributions $\left(=\frac{1}{d}\sum\limits_{k=1}^{d}H_{X_{k}}(t_{1},t_{2})\right)$
$H_{\text{normalized}}(t_{1},t_{2})$ |  Hellinger distance normalized with respect to the dimension $(=\sqrt{1-BC^{1/d}})$
$\Theta(\cdot)$ |  copula function
$f^{\text{BV}}(\cdot)$ |  Oracle function in the Bernstein-Vazirani problem
$r$ |  secret $n$-bit string in the Bernstein-Vazirani problem
$U_{f}$ |  Unitary for the Oracle function in the Bernstein-Vazirani problem
$P_{a},P_{C},P_{T}$ |  Pauli gates
$\mathcal{Q}$ |  quasi probability distribution
$\\{\eta_{{}_{k}}\\}$ |  generally refers to a PEC linear combination coefficient
$\tilde{\mathcal{B}}_{i}$ |  noisy basis superoperators in PEC

##  Chapter 1 Introduction

###  1.1 Background

####  1.1.1 Promise of quantum computing

Error-resilient quantum computing holds great promise, offering significant
advancements over conventional computing. Once realized, it is expected to
super-polynomially reduce execution time, energy consumption, and memory
storage needs compared to conventional state-of-the-art computers [1]. The
potential impact of error-resilient quantum computing includes revolutionizing
scientific applications such as simulating many-body quantum systems [2],
solving large-scale optimization problems [3], efficiently sampling high-
dimensional probability distributions [4], factorizing large integers, and
enhancing the security of communication networks [5]. Consequently, this
technology is expected to be disruptive to sectors such as telecommunications,
cyber-security, pharmaceuticals, logistics, supply chain management,
artificial intelligence, and materials science [6].

####  1.1.2 Quantum computing vs classical computing

Classical mechanics, rooted in the laws of Newtonian physics, has served as a
successful framework for understanding the macroscopic world for centuries.
However, when examining the behavior of particles at the atomic and subatomic
scales, classical mechanics began to exhibit limitations and inconsistencies.
For instance, classical mechanics predicted absurd outcomes, like suggesting
that a blackbody emits an infinite amount of energy across all wavelengths.

As scientists delved deeper into the microscopic realm, counter-intuitive
phenomena such as wave-particle duality, quantized energy levels, and non-
locality emerged, challenging the classical paradigm. These challenges
necessitated the development of quantum mechanics, which offered a novel and
revolutionary approach to describe the behavior of particles at the quantum
level [7].

Quantum mechanics introduced probabilistic interpretations, superposition
states, and entanglement, providing a more accurate and comprehensive
understanding of the intricate workings of nature [8]. For example, consider a
particle moving in one dimension under the influence of a conservative force,
such as a harmonic oscillator. In classical mechanics, we can describe the
particle’s motion using Newton’s second law and the equation of motion for a
harmonic oscillator, which yields a sinusoidal trajectory and continuous
energy levels. In contrast, in quantum mechanics, we describe the particle
using the Schrödinger equation for a harmonic oscillator, which results in
quantized energy levels and wave functions corresponding to discrete energy
states. However, as the value of the Planck’s constant (h), representing the
fundamental scale of quantum mechanics, approaches zero ($h\rightarrow 0$),
the quantum system converges towards the classical limit. In this limit, the
quantized energy levels of the quantum harmonic oscillator become densely
spaced and form a continuous energy spectrum, matching the classical behavior.
The quantum wave function also converges to the classical trajectory, and the
classical and quantum results become indistinguishable.

This convergence phenomenon is known as the correspondence principle [9],
where classical mechanics emerges as the limiting case of quantum mechanics at
large scales or when the quantum effects become negligible. Understanding the
system’s scale is crucial in selecting the appropriate framework for a given
physical problem, with classical mechanics suitable for macroscopic objects
with well-defined trajectories and negligible quantum effects, while quantum
mechanics is employed for microscopic particles, providing a more accurate
description of phenomena like wave-particle duality, quantization of energy
levels, and quantum entanglement.

Quantum mechanics is built upon several fundamental postulates [10, 11] that
were formulated to address experimental observations in the early 20th
century. These postulates provide the framework for understanding the behavior
of quantum systems. The first postulate states that every quantum system is
associated with a complex Hilbert space [12]. The quantum state of a system is
described by a density operator [13, 12], often denoted by $\rho$, which
belongs to the set of density operators defined on the Hilbert space. The
second postulate deals with measurements [13]. When a measurement is performed
on a quantum system, it can have random outcomes with finite probabilities.
The third postulate connects the quantum state with the measurement outcomes.
When a measurement $M$ is performed on a quantum system in the state $\rho$,
then the observed outcome is a realization of the random variable $M\rho$
[11]. The fourth postulate addresses composite quantum systems. When we have
two quantum systems, each associated with its own Hilbert space, say
$\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, the combined Hilbert space of the
composite system is given by the tensor product
$\mathcal{H}_{1}\otimes\mathcal{H}_{2}$. This tensor product construction
allows us to represent the joint states of the individual systems [9].

These fundamental postulates provide the tools to analyze quantum systems. The
concepts of quantum state, measurement, and composite systems are the key
building blocks for understanding the intriguing phenomena that occur in the
quantum world.

Classical computing is the conventional form of computing that relies on
classical bits, represented by the binary numbers 0 and 1. Classical computers
crucially rely on components like transistors, which function based on quantum
mechanical principles. However, despite the quantum nature of the transistors,
the interactions between these components within a classical computer follow a
classical framework. This distinction can sometimes lead to confusion, as it
may seem unsatisfactory to say that classical computers, which are built using
components based on quantum principles, operate according to classical laws.
But while specific underlying components leverage quantum phenomena, the
interactions between these components, and the analysis of the data produced
by them, are modeled satisfactorily using classical physics and classical
information theory.

Analogously, the link between quantum computing [14] and quantum mechanics is
fundamental, as the fundamental interaction between the components of quantum
computers depend on quantum mechanical phenomenon [11] such as superposition
and entanglement. Quantum computations use quantum bits or qubits [15] as the
basic unit of information. Superposition enables qubits to exist in multiple
states simultaneously, allowing quantum computers to perform parallel
calculations. Entanglement [16] creates strong correlations between qubits,
even when they are physically separated, potentially leading to increased
computational power. Different qubit technologies, such as superconducting
qubits, trapped ions, and topological qubits, each utilize distinct quantum
phenomena. Researchers continuously draw upon quantum mechanical principles to
optimize performance, addressing challenges like quantum decoherence [17, 18,
19, 20] resulting from interactions with the external environment.

Classical computing remains the practical and efficient choice for many
computing needs. In fact, the computing power of classical computers has been
doubling every two years, as per Moore’s law. But still classical computing
faces limitations with computationally challenging problems that grow
exponentially with problem size such as unstructured large-scale optimization,
factorization of large integers, and simulation of many-body systems. The
quantum equivalent of Moore’s law [11] states that adding just one perfect
qubit to a quantum computer doubles its computational capability. Therefore,
to match the progress of classical computers, a single error-resilient qubit
needs to be integrated into quantum computers every two years.

However, quantum algorithms [10] are more efficient for tackling only a sub-
set of classically computationally challenging problems (not universally).
Example of problems that have an efficient quantum algorithm include prime
factorization and discrete logarithm, both of which were developed by Shor. An
efficient algorithm [21] operates within a time frame that corresponds to a
polynomial function of the problem size, whereas an inefficient algorithm
takes time that corresponds to a super-polynomial function of the problem
size.

The study of algorithm efficiency [21] is a fundamental aspect of complexity
theory [22], a branch of computer science. Problems are categorized based on
their resource requirements, such as time and memory. For instance, problems
solvable in polynomial time by classical computers fall into the class P,
whereas those with solutions verifiable in polynomial time belong to class NP.
While it’s evident that P is a subset of NP, the question of whether there are
problems in NP not in P remains unresolved.

In the realm of quantum computing, problems solvable by quantum algorithms
within polynomial time (with bounded error probability) are classified as
belonging to class BQP [11]. It has not been formally established whether BQP
contains P. So, we are not certain that quantum computing contains classical
computing as a special case, but evidence supports this assertion.

BQP of course contains QP which represents problems that a quantum computer
can solve with a 100% probability of success in polynomial time. Examples that
belong to the BQP class include Deutsch-Jozsa, Bernstein-Vazirani, and Simon’s
algorithm.

Why can we not simulate quantum computation using classical computers? The
reason lies in the exponential space and time complexity involved in storing
quantum gates as classical matrices and tracking entangled qubits after
logical operations. For instance, even a system with just 500 atoms would need
$2^{500}$ complex coefficients for perfect description. Attempting
computations with such requirements would overwhelm classical computers.
However, quantum computers excel in simulating such scenarios by storing,
representing, and evolving states as native quantum states on qubit registers,
bypassing the need for managing $2^{500}$ complex floating-point numbers with
limited precision

Note that, even if we restricted ourselves to a small scale register, then
also a perfect quantum computer can never be built using classical computers
because quantum measurement cannot be perfectly simulated as there is no
perfect random number generator.

In classical computing, logic gates are basic building blocks that manipulate
classical bits (0s and 1s) to perform logical operations. Two examples of
classical logic gates are the AND gate and the NOT gate. The AND gate takes
two input bits, and its output is 1 (true) only when both input bits are 1;
otherwise, the output is 0 (false). The NOT gate takes a single input bit and
produces the opposite value as output.

Input 1 | Input 2 | AND Gate Output
---|---|---
0 | 0 | 0
0 | 1 | 0
1 | 0 | 0
1 | 1 | 1

Input | NOT Gate Output
---|---
0 | 1
1 | 0

In quantum computing, instead of bits, we have qubits, which are represented
as vectors in a complex two-dimensional Hilbert space. The basis vectors
spanning this space are commonly expressed in three ways: the
$\mathds{Z}$-Basis (also called computational basis or standard basis), the
$\mathds{X}$-basis, and the $\mathds{Y}$-basis. Each of these bases is made up
of orthogonal vectors in two dimensions, ideal for representing a two-level
quantum system in a two-dimensional vector space.

The need for these bases becomes clear when studying the underlying physics of
quantum systems. Take, for example, the spin of an electron in a magnetic
field. The electron’s spin resembles a minuscule magnetic moment. Its
orientation relative to the magnetic field influences the electron’s energy.

For the $\mathds{Z}$-Basis, consider a magnetic field applied vertically. An
electron’s spin might align with this field, represented by the lower energy
state $\ket{0}$, or it could oppose the field, corresponding to the higher
energy state $\ket{1}$. This basis provides an intuitive way to think about
qubits, likening the lower and higher energy states to the classical binary
values of $0$ and $1$ respectively. A qubit’s state can be in a superposition
of both $\ket{0}$ and $\ket{1}$, expressed as $c_{1}\ket{0}+c_{2}\ket{1}$,
where $c_{1}$ and $c_{2}$ are complex amplitudes satisfying:

$|c_{1}|^{2}+|c_{2}|^{2}=1.$ (1.1)

However, if we change our perspective and measure the spin horizontally, along
the $\mathds{X}$-axis, the electron’s spin might point left or right. These
orientations, when related back to the $\mathds{Z}$-Basis, are actually
superpositions of the $\ket{0}$ and $\ket{1}$ states. These superpositions,

$\ket{+}=\frac{\ket{0}}{\sqrt{2}}+\frac{\ket{1}}{\sqrt{2}},$ (1.2)

and

$\ket{-}=\frac{\ket{0}}{\sqrt{2}}-\frac{\ket{1}}{\sqrt{2}},$ (1.3)

define the $\mathds{X}$-basis.

The $\mathds{Y}$-basis offers yet another viewpoint. Measuring perpendicular
to both the $\mathds{X}$ and $\mathds{Z}$ axes gives states that are complex
superpositions of the $\mathds{Z}$-Basis:

$\ket{i}=\frac{\ket{0}}{\sqrt{2}}+i\frac{\ket{1}}{\sqrt{2}},$ (1.4)

and

$\ket{-i}=\frac{\ket{0}}{\sqrt{2}}-i\frac{\ket{1}}{\sqrt{2}},$ (1.5)

Quantum operations can be visualized as rotations around these different basis
axes, offering valuable insights for designing quantum algorithms.
Furthermore, using different bases for measurements can be instrumental in
pinpointing various types of errors in quantum systems.

Quantum gates, represented by unitary matrices, manipulate qubits, with
essential examples including the Pauli gates ($\mathds{I}$, $\mathds{X}$,
$\mathds{Y}$, and $\mathds{Z}$), the Hadamard gate ($\mathds{H}$), and the
two-qubit entangling CNOT gate ($\mathds{U}_{\text{CNOT}}$). The
Pauli-$\mathds{X}$ gate acts as the quantum analog of the NOT gate, flipping
the qubit’s state between $\ket{0}$ and $\ket{1}$. The Pauli-$\mathds{Y}$ gate
transforms $\ket{0}$ to $i\ket{1}$ and $\ket{1}$ to $-i\ket{0}$. The
Pauli-$\mathds{Z}$ gate introduces a relative phase shift between the basis
states of a qubit (it leaves $\ket{0}$ unchanged and flips the sign of
$\ket{1}$). The Hadamard gate puts a qubit in an equal superposition of
$\ket{0}$ and $\ket{1}$. The CNOT gate is a two-qubit gate that flips the
target qubit if and only if the control qubit is in state $\ket{1}$.

Pauli-$\mathds{X}$ Gate

Input | Output
---|---
$\ket{0}$ | $\ket{1}$
$\ket{1}$ | $\ket{0}$

Pauli-$\mathds{Y}$ Gate

Input | Output
---|---
$\ket{0}$ | $i\ket{1}$
$\ket{1}$ | $-i\ket{0}$

Pauli-$\mathds{Z}$ Gate

Input | Output
---|---
$\ket{0}$ | $\ket{0}$
$\ket{1}$ | -$\ket{1}$

Hadamard Gate

Input | Output
---|---
$\ket{0}$ | $(\ket{0}+\ket{1})/\sqrt{2}$
$\ket{1}$ | $(\ket{0}-\ket{1})/\sqrt{2}$

CNOT gate
Control qubit Target qubit Output (state of target qubit) $\ket{0}$ $\ket{0}$
$\ket{0}$ $\ket{0}$ $\ket{1}$ $\ket{1}$ $\ket{1}$ $\ket{0}$ $\ket{1}$
$\ket{1}$ $\ket{1}$ $\ket{0}$

In classical circuits, information flows through the movement of electrons
from one transistor to another in a well-defined spatial layout. Classical
logic gates, such as AND, OR, and NOT, manipulate classical bits (0 or 1) and
perform logical operations. In contrast, quantum circuits process information
using qubits. Quantum circuits evolve the quantum state in situ, which modify
the quantum amplitudes and phases of the state of the quantum register.
Quantum gates, like Pauli-$\mathds{X}$, CNOT, and Hadamard, implement quantum
algorithms by performing quantum operations on the qubits.

Designing quantum algorithms [10] is significantly more challenging than
classical algorithms [21] for several reasons. Firstly, quantum computing
requires a departure from classical intuition, as quantum phenomena behave
differently from classical physics. For instance, computer scientists
experienced in conventional parallel programming understand the challenges
associated with designing algorithms that can effectively harness GPU
parallelism. They can thus empathize with the complexity of leveraging
computing power through superposition, a form of parallel computing utilized
in quantum systems. However, phenomena such as entanglement and quantum
interference present unique opportunities that lack analogues in classical
algorithm development. Secondly, to demonstrate the utility of a quantum
algorithm, it must be more efficient than the best-known classical algorithm
for a specific problem. The competitiveness of the latter introduces a moving
target for quantum algorithm developers, where the best classical algorithms
keep evolving, demanding continuous advancements in quantum algorithms to
maintain claims of utility. Lastly, mapping a real-world (often classical)
use-case into a quantum representation is a non-trivial task.

A quantum program is a series of instructions that can be executed by a
quantum device in a specific sequence to perform a specific task. These
instructions are typically written in a high-level programming language like
Qiskit [23], designed to be readable and writable by humans.

A classical computer architecture consists of several key components that work
together to perform various computational tasks. At its core, a classical
computer contains a central processing unit (CPU). The CPU executes
instructions, and coordinates data movement. It is supported by memory units,
including random-access memory (RAM) and cache memory, where data and
instructions are temporarily stored for faster access. The architecture also
includes input and output (I/O) devices, such as keyboards, mice, monitors,
and storage devices like hard drives or solid-state drives. These allow users
to interact with the computer and store data for future use. The CPU
communicates with other components via buses, which are pathways that transfer
data and control signals between different parts of the computer. The system
clock generates regular pulses that synchronize the activities of various
components, ensuring smooth coordination of operations. Moreover, a classical
computer architecture often involves a graphics processing unit (GPU)
dedicated to handling graphics-intensive tasks, such as rendering images and
videos.

The development and elaboration of a quantum computer architecture [24, 25,
26] are still in their early stages, mainly because we have not achieved
fault-tolerant quantum computing yet. The full stack architecture will need
focus beyond physical layer and must include error correction, feedback
stabilization, hardware-aware compilation, logical level compilation, circuit
optimization, application layer, and user interface. While progress is being
made towards this vision, it remains a distant goal.

####  1.1.3 The problem of noise

The behavior of an ideal quantum computer can be modeled as follows. An
$n$-qubit register spans a complex Hilbert Space denoted by
$(\mathbb{C}^{2})^{\otimes n}$. The initial state of the register can be
represented as a tensor:

$\ket{\psi}=\ket{0}^{\otimes n}.$ (1.6)

A logical operation on the register state (also called a quantum gate) can be
represented by a linear, unitary operation U:

$\ket{\psi}\rightarrow U\ket{\psi}.$ (1.7)

A measurement reads out a $n$-bit string v in the computational basis:

$Pr(v)=|\braket{v}{\psi}|^{2},$ (1.8)

where $v\in\\{0,1\\}^{\otimes n}$

Unlike modern classical computers with extremely low failure rates (e.g.,
$10^{-17}$ or less), superconducting quantum computers exhibit higher gate-
level failure rates (nearly $0.01$) [27, 28]. Thus, it is imperative to
investigate the causes of noise and mitigate and correct them so that quantum
computers can provide correct results.

Noise, in the context of this dissertation, refers to deviations from the
ideal description of a quantum computer. Practical efforts to build quantum
computers introduce noise, which affects technologies like superconducting
qubits, trapped ions, and silicon quantum dots. Our focus is primarily on
superconducting quantum computers [29, 30]. The underlying noisy processes
that impact such a computer can be classified into three groups: noise
affecting the the quantum register (such as leakage [31], undesired coupling,
decay processes, non-uniformity, and cross-talk [32]), noise affecting the
quantum operations (such as pulse distortion, attenuation, drift, and mis-
calibration), and noise in the thermodynamic isolation system (due to issues
with dilution refrigerators, vacuum chambers, shields, and vibration
suppression mechanisms [33]).

Quantum computers today are referred to as existing in the NISQ era [34],
which stands for noisy intermediate-scale quantum. The noise threshold for
NISQ is defined by a single-qubit error rate being worse than $10^{-4}$. The
intermediate-scale label is often associated with having fewer than a hundred
thousand qubits. Computing done with NISQ devices is called NISQ computing.
The bare minimum requirements for NISQ computing includes: quantum registers
for storing data, quantum gates to execute logical operations on the
registers, and a measurement interface for extracting the computation
outcomes. The field has witnessed rapid advancements, with NISQ devices now
operating as systems with hundreds of interacting qubits. Remarkably, the
field is already witnessing a transition towards a phase where NISQ devices
are performing scientific computations at a scale that rivals classical
supercomputers in terms of computational power [35]. These experiments apply
error mitigation techniques to the outcomes of the noisy computations
performed on the NISQ devices.

Note that error mitigation and error correction are distinct strategies to
tackle challenges arising from noise: mitigation employs statistical
techniques focused on minimizing noise effects rather than eliminating errors
entirely, such as zero noise extrapolation and probabilistic error
cancellation, while error correction seeks to actively detect and rectify
errors during computations with an objective of fault tolerance, but requires
the device noise to be below a threshold, which has not been achieved yet. One
example of the unique challenges in the area of quantum error-correction is
the no-cloning theorem states that it’s impossible to create an exact
duplicate of an unknown quantum state. This makes it challenging to
incorporate redundancy into quantum computing systems to protect against
information corruption.

###  1.2 Research focus

The fact that noise exhibits non-stationarity, underscores the core motivation
of this dissertation. A quantum noise channel is often used to model how
quantum information becomes distorted during its passage through a physical
system. It describes how interactions with the environment can modify the
quantum state of a system. Examples include depolarizing channel, Pauli noise
channel, amplitude damping channel, phase damping channel, and the SPAM noise
channel.

The term SPAM denotes state preparation and measurement. The statistics of the
SPAM error channel are commonly quantified using the SPAM fidelity, a metric
that evaluates a device’s ability to prepare and measure a qubit.
Specifically, SPAM fidelity quantifies the likelihood that the device readies
the qubit (or a set of qubits) in the desired state and subsequently measures
it in that same state. Fig. 1.3 \- 1.5 depicts the time-varying noise
densities of SPAM fidelity.

Concurrent experimental studies [17, 18, 19, 36, 37] support our concern that
the assumption of fixed and invariant parameters for quantum processors is
flawed. In fact, these studies show that the noise parameters can exhibit time
variations of up to 50% of their mean value within an hour. Spatially varying
noise in quantum devices have also been extensively studied, including the
role of circuit geometry [38, 39] and cross-talk between neighbouring qubits
[40, 41]. Unlike temporal variations in noise, the effects of which are
magnified by the complexity of a quantum circuit, the spatial variations, such
as seen in Fig. 1.2, are dependent on the geometry and scale of the circuit
implementation.

To understand the impact of such non-stationary noise on program outcomes,
consider Fig. 1.1 which depicts the time-varying histogram obtained from IBM’s
ibm_mumbai device for the Bernstein-Vazirani circuit. It is evident that the
associated error bars on a particular day do not provide insights into results
from a different date, highlighting the problem of reproducibility of results
in quantum computing today.

The causes of non-stationary noise are not fully understood but are believed
to stem from TLS defects in transmon registers which might be arising from
deviations in crystalline order. The current consensus attributes these
defects to the presence of certain oxides on the superconductors’ surface [42,
36]. Thus, static quantum channel models do not accurately capture the
dynamics in realistic quantum computations, particularly in superconducting
qubits. Cosmic rays [43, 44] also contribute by ionizing the substrate upon
impact, leading to the emission of high-energy phonons, which in turn triggers
a burst of quasi-particles. These quasi-particles disrupt qubit coherence
across the device. It has been shown that quantum computers can experience
catastrophic errors in multi-qubit registers approximately every 10 seconds
due to cosmic rays originating from outer space [44]. Studies that address
non-stationary noise in superconducting quantum computers include
investigations on output reproducibility [45], noise modeling [46], tracking
the non-stationary profile of quantum noise [47], and quantum error mitigation
using continuous control [48]. Non-stationary quantum channel models have been
proposed [46, 49, 50, 51] that use stochastic processes. Our dissertation
focuses on understanding the effect of non-stationary noise on program
outcomes, as well as devising strategies to address it. Specifically, we model
the noise channel as a random variable and implement adaptive methods to
manage it.

Before we can begin, we need to fix the precise language for performance
assessment. This task is not trivial due to the complexity inherent in quantum
technology, which both distinguishes it from classical computing and hinders
its rigorous checking [52, 53]. Challenges include the inherent randomness in
quantum measurements, error accumulation without clear source attribution, the
curse of dimensionality, and the inability to step-through program execution
in quantum circuits.

In fact, the performance evaluation [54] of noisy quantum computations is a
vast topic that is crucial for several additional reasons, apart from our
motivation of studying the impact of non-stationary noise. Firstly, as quantum
computing is still in its early stages [55, 56, 57], understanding the sources
of errors and noise is vital. Through rigorous evaluation, researchers can
model [58], identify, model and quantify sources of noise, such as
decoherence, gate errors, and readout errors. Secondly, this understanding is
essential for developing error mitigation techniques [59]. Thirdly,
reproducibility of results is critical, and rigorous performance evaluation
ensures experiments can be replicated by other researchers, contributing to
the validation and verification of quantum algorithms. An additional challenge
is the diverse range of terms encountered in quantum computing today which can
blur distinctions between them, making it challenging to appreciate their
nuanced differences. Examples include verification [60, 61, 62] (ensuring
correct transpilation), validation [63, 64] (validating correctness of output
or the quantum nature of a device), benchmarking [65, 66, 67] (assigning a
performance measure to a processor), accreditation [68], and certification
[54]. Thus out first task is to precisely define computational accuracy,
result reproducibility, device reliability, and observable stability in the
presence of non-stationary noise.

Our next objective is to experimentally assess hardware reliability, with a
particular emphasis on analyzing spatial and temporal variations in noise
statistics. IBM [55] has introduced a range of processors in recent years,
each with an expanding register size. These include the Canary processors with
2-16 qubits, Falcon processors with 27 qubits, Egret processors with 33
qubits, Hummingbird processors with 65 qubits, Eagle processors with 127
qubits, and Osprey processors with 433 qubits. Quantifying the spatial and
temporal reliability of these quantum computers is crucial to understand
system-wide performance changes over time. This evaluation should encompass
both component-level metrics, such as individual gates and qubits, and
composite-level metrics, such as circuits, to assess the degree of non-
stationarity in noise and its implications on program outcomes. Holistically
measuring reliability at the circuit level is essential, as examining
thousands or millions of qubits and gates may not provide conclusive insights
at the application level [69].

Our third objective in this dissertation is to establish stability bounds for
error-mitigated outcomes affected. We aim to determine the upper and lower
bounds for our performance evaluation metrics. Our inquiries include
determining the minimum sample size necessary to ensure histogram
reproducibility with a confidence level of $1-\delta$, bounding outcome
stability based on the variation in time-varying noise densities, and
establishing reliability bounds to achieve stable outcomes.

Numerous studies on noise modeling in quantum computing systems have
highlighted the challenges associated with noise estimation [70, 71]. A
natural question arises regarding how can we effectively counteract the
detrimental impacts of non-stationary noise using adaptive algorithms? In the
final chapter of our dissertation, we consolidate the various elements of our
investigation in the context of adaptive probabilistic error cancellation
(PEC) [72]. Our focus is on enhancing the stability of PEC outcomes, using a
Bayesian [73, 74, 75, 76] updating of the quasi-probability distributions, in
the presence of non-stationary noise.

This research focuses on a limited scope. Firstly, the experiments exclusively
uses the superconducting platforms provided by IBM. Other platforms such as
trapped ion, neutral atom, photonic, or quantum dot are not considered.
Secondly, our performance evaluation framework mainly concerns with the output
measured in computational basis and how it is impacted by device noise.
Thirdly, not all the superconducting devices provided by IBM have been
characterized; only a sub-set of the devices, mainly ibm_kolkata, ibm_mumbai,
ibm_washington, ibm_toronto, and ibm_yorktown, are used in this study. Lastly,
we do not focus on the problem of optimal selection of a statistical model for
a given noisy device. Instead, a generic error channel formalism is employed
whenever possible. For verifying the theory, quantum noise channel models like
the Pauli noise [77] channel is chosen, with a specific focus on parameter
estimation. The research does not explore the question of identifying the best
noise model for a given device.

For our research, we have made use of the daily characterization data stored
on IBM’s servers as-is. For data at time-scales of minutes and below, we
collected the data ourselves and offer all the associated collection and
preparation software. For the latter case, the data collection frequency was
limited by network time lags and constraints in the qiskit software, such as
the maximum number of circuits and shots allowed. These limitations have been
gradually improving over time

###  1.3 Organization and notation

The document is structured as follows. Chapter 2 provides background on noise
in quantum computing, quantum channel modeling, and experimental
characterization of quantum decoherence, emphasizing the non-stationary
statistics of noise. Chapter 3 establishes a systematic framework for
assessing noisy quantum computer performance. Chapter 4 focuses on the testing
of reliability. The evaluation encompasses both component-level metrics (such
as individual gates and qubits), and composite-level metrics (such as
circuits). Chapter 5 seeks to determine the bounds on the assessment framework
developed in Chapter 3, using available noise characterization data.
Specifically, it discusses how to bound outcome stability in terms of the
distance between time-varying noise densities [78]. Chapter 6 explores methods
to improve accuracy in the presence of non-stationary noise. Chapter 7 brings
together the various concepts discussed till date in the context of adaptive
probabilistic error cancellation. Chapter 8 provides concluding remarks.

Notations in this dissertation vary in meaning depending on font, although it
should be clear from the context (see Table LABEL:tab:Notations ). This was
needed because the work draws upon concepts from physics, information theory,
computer science, and statistics, each of which has established conventions.
$\ket{\psi}$ represents a pure quantum state, and $\rho$ represents density
matrices. The symbol $\otimes$ signifies a tensor operator, and Tr$(\cdot)$ is
an abbreviation for the trace operator. The system Hamiltonian is denoted as
$\Omega$, with its eigenvalues represented by $\lambda$ (however $\omega$
signifies angular frequency). Quantum observables are typically denoted by
$\hat{O}$, and the uppercase letter $U$ typically stands for a unitary matrix.
A single-qubit rotation by an angle $\theta$ on the Bloch sphere is often
denoted as $R(\theta)$. However, the uppercase Greek letter $\Theta(\cdot)$
denotes the copula [79] function from statistics.
$\mathcal{E}_{\textrm{x}}(\cdot)$ represents a quantum noise channel, while
$\Lambda$ symbolizes the SPAM noise channel, a classical channel operating on
probabilities. The canonical Pauli matrices are denoted as
$\mathds{X},\mathds{Y},\mathds{Z},\mathds{I}$. Note that we do not use the
small Greek letter $\sigma$ for Pauli matrices, reserving it for standard
deviation instead. The Pearson correlation matrix is denoted by the capital
Greek letter $\Sigma$. Note that the small z signifies a standard normal
variable. The identity matrix is represented as $\mathds{I}$, while
$\mathcal{I}(X,Y)$ in calibrated font signifies the mutual information between
random variables X and Y. The letter x typically signifies noise parameter(s),
whereas an uppercase X corresponds to a specific realization of x. If x is not
deterministic, then $f(\textrm{x})$ denotes the probability distribution of x.
This distribution can exhibit temporal fluctuations, denoted as
$f(\textrm{x};t)$, with its cumulative distribution function indicated by
$F(\textrm{x};t)$. Additionally, the curly capital $\mathcal{F}$ stands for
Fisher Information, while $f_{s}$ with the subscript $s$ represents the data
sampling frequency. The symbol $\Pi_{r}$ with a subscript is reserved for the
projector operator onto the eigenstate $\ket{\lambda_{r}}$ (however, when
presented without a subscript, $\Pi$ signifies the normal product operator).
The measurements are conducted in the computational basis (or $\mathds{Z}$
basis), resulting in qubits yielding classical bits. We employ the notation
$b_{i}(t)$ to represent the observed classical bit value on qubit $i$ at time
$t$. The state of an $n$-qubit quantum register is denoted by
$\ket{v}=\ket{v_{n-1}\cdots v_{0}}$, with $n$ generally denoting the quantum
register size. Upon measurement, this state yields an $n$-bit string, with
each $v_{i}\in{0,1}$. We utilize the symbol $W$ to denote the dataset
consisting of collected bit-strings from repeated circuit executions. The
total number of samples collected is typically denoted as $L$, where $l$
denotes the l-th circuit execution (however, note that the curly $\mathcal{L}$
represents the likelihood function). The Hellinger distance between
probability distributions is denoted as $H$, while the curly capital
$\mathcal{H}$ is exclusively reserved for entropy. The symbol $\mathds{H}$
represents the Hadamard gate. A quantum circuit is represented by the capital
$C$, and while, the small c usually signifies the control qubit in a CNOT
gate. Additionally, $c$ serves as a constant in certain information theoretic
results. The symbols $\alpha$ and $\beta$ may assume different meanings
depending on the context, referring either to quantum amplitudes for $\ket{0}$
and $\ket{1}$ or the parameters of the beta distribution, both of which find
common usage in this dissertation. The small greek $\gamma$ is used as a proxy
parameter encapsulating various device noise parameters, while the capital
$\Gamma$ is reserved for the gamma function. The symbol $\eta$ typically
represents a PEC linear combination coefficient. Capital $D$ denotes a noisy
quantum device, while the small letter $d$ typically dimensionality. In the
context of the Bernstein-Vazirani problem, $r$ is used to denote the secret
$n$-bit string, with the latter problem being extensively employed as an
illustrative quantum circuit. The small letter $s$ is used to denote the
stability metric. Absolute time is typically represented by the small letter
$t$, while time intervals are denoted as $\tau$. For instance, $\tau_{c}$
signifies the circuit execution time, and $\tau_{N}$ denotes network delay.
Time duration is represented as $\delta t$. However, the capital letter $T$ is
primarily reserved for parameters related to decoherence characterization,
such as $T_{1}$ and $T_{2}$ (representing qubit relaxation and dephasing time,
respectively). In some instances, $T$ (without a subscript) is used to denote
the target qubit in a CNOT gate, which is generally clear from the context.
Finally, the non-standard abbreviations that have been used in this
dissertation are: NSN (non-stationary noise), PEC (probabilistic error
cancellation), BC (Bhattacharyya coefficient), BV (Bernstein-Vazirani), SPAM
(State preparation and measurement) and WSS (Wide-sense stationary).

Figure 1.1:  Impact of the non-stationary noise, leading to irreproducible
outcomes. The results show histograms (after state preparation and measurement
(SPAM) noise mitigation) upon executing the Bernstein-Vazirani circuit on
ibm_mumbai.  Figure 1.2: Evidence of spatial non-stationarity in the mean
values for quibt decoherence times $T_{1}$ and $T_{2}$ for the 127-qubits of
the ibm_washington device, generated on 14 Jan 2023 10:20 PM UTC.

Figure 1.3: state preparation and measurement (SPAM) fidelity distributions on
ibm_toronto for qubits $0-8$ as measured on 8 April 8 2021, between
8:00-10:00pm (UTC-05:00). Figure 1.4:  state preparation and measurement
(SPAM) fidelity distributions on ibm_toronto for qubits $9-17$ as measured on
8 April 8 2021, between 8:00-10:00pm (UTC-05:00). Figure 1.5:  state
preparation and measurement (SPAM) fidelity distributions on ibm_toronto for
qubits $18-26$ as measured on 8 April 8 2021, between 8:00-10:00pm
(UTC-05:00).

##  Chapter 2 Noise in quantum computing

Practical efforts to realize a quantum computer (e.g. transmons, trapped ions,
silicon quantum dots [18, 80, 81, 1]) introduce various physical processes,
referred to as noise, which deviate from the ideal description of a quantum
computer. Unlike modern classical computers, which boast device components
with extremely low failure rates (e.g., $10^{-17}$ or less), the current
state-of-the-art quantum computers exhibit higher gate-level failure rates
(e.g., $10^{-2}$). In this dissertation, we mainly use experimental data from
transmon[30, 33] based realizations of a quantum computer. Transmon qubits are
a variant of superconducting charge qubits designed to reduce sensitivity to
charge noise[29].

###  2.1 Physical sources of noise

The various noise processes [17, 18, 19] can be classified into three groups:

####  2.1.1 Quantum register

One of the pathways for noise is the implementation of the quantum register,
which encounter phenomena such as: (i) Leakage i.e. unintended energy states
outside the computational subspace, (ii) Undesired coupling to the external
environment (such as spurious charge, magnetic fields, stray photons, lattice
vibrations (phonons), nuclear spins) leading to loss of coherence, (iii)
Spontaneous decay processes that transition a qubit from an excited state to a
lower energy state, (iv) Non-uniformity in qubit’s coupling strength to the
control field, results in errors in quantum gates, (v) Inter-qubit cross-talk
arising from shared control lines or capacitive coupling between neighboring
qubits. In trapped ion systems, cross-talk could arises from motional coupling
between ions, affecting the states of neighboring ions.

####  2.1.2 Control system

The quantum register undergoes four fundamental control operations:
initialization or reset, measurement, single-qubit rotation gates, and 2-qubit
entangling gates. These operations require the application of precisely
calibrated control pulses on the qubits. Imperfections in the control system
used for logic implementation can arise from several sources. Firstly, pulse
distortion occurs when the desired shape and duration of pulses encoding
quantum information are altered due to the finite time resolution and
frequency response limitations, as well as pulse timing errors. Secondly,
control pulses may experience attenuation caused by electromagnetic
interference and material imperfections in the quantum system’s vicinity.
Thirdly, qubits can drift either physically (in the case of trapped ions) or
in parameter space. Fourthly, the noise may be an effect mis-calibration.

####  2.1.3 Thermodynamic isolation system

Transmon qubits, a specific kind of superconducting qubit, require cooling to
approximately 10 milli-kelvin in order to mitigate the presence of thermal
noise. The thermodynamic isolation system [82] helps achieve this using a
system of dilution refrigerators, vacuum chambers, electromagnetic shields,
and vibration suppression mechanisms. The dilution refrigerators employ a
multi-stage cooling process that gradually reaches colder temperatures, using
substances like liquid helium to progressively lower the temperature. The
vacuum chambers effectively eliminate gas molecules and particles that could
potentially couple with the qubits. The electromagnetic shields are
responsible for blocking external radiation and fields from disturbing the
quantum state. Inadequate electromagnetic shielding could allow disruptive
external radiation to interfere with the qubits. The vibration suppression
systems minimize mechanical vibrations and movements that could potentially
jeopardize the quantum states. Noise from imperfect thermodynamic control
systems can be non-Markovian in nature, which are difficult to rectify using
quantum error correction tools.

###  2.2 Cause of non-stationarity

Non-stationary noise refers to noise in a quantum system that exhibits time-
varying statistical properties. The temporal fluctuations of the mean and
variance of the energy relaxation times ($T_{1}$), dephasing times ($T_{2}$),
and qubit frequencies, are well-studied topics[18, 19, 36, 37] that suggest
suggest that noise in NISQ[34] devices can fluctuate unpredictably. For
example, $T_{1}$ times have been found to fluctuate by approximately 50
percent within an hour [18]. Similarly, many advances have also been made for
spatially varying noise in quantum devices [38], its effect on the choice of
circuit geometry [39], as well as the interplay with cross-talk between qubits
for single and two-qubit gates [40, 41].

The causes and mechanisms behind quantum noise non-stationarity are poorly
understood. In transmon registers, potential sources of fluctuations include
TLS (two-level system) defects, quasi-particles, parasitic microwave modes,
phonons, nuclear spins, paramagnetic impurities, spurious resonances, critical
current noise, background charges, gate voltage fluctuations, and the
electromagnetic environment [57]. Among these, TLS defects have been
identified as the primary cause of decoherence [17, 18, 19]. These defects
arise from deviations from crystalline order in the naturally occurring oxide
layers of transmons, resulting in trapped charges, dangling bonds, tunneling
atoms, or collective motion of molecules.

The findings not only highlight the necessity for frequent re-calibration in
qubit setups but also question the reproducibility of device
characterizations, and their use in error mitigation. Consequently, modeling
time-varying quantum noise has become an active area of interest [46], such as
through the inclusion of $T_{1}$ and $T_{2}$ fluctuations in quantum channel
models to investigate the concept of time-varying quantum channels (TVQC).

###  2.3 Decoherence characterization

Decoherence refers to loss of unitarity in state evolution. The traditional
definition of decoherence, which describes the decay of off-diagonal terms in
the density matrix, is now referred to as dephasing and considered one kind of
decoherence [57]. Decoherence studies typically focus on three metrics:
transverse relaxation time $(T_{1})$, longitudinal relaxation time $(T_{2})$,
and dephasing time $(T_{\phi})$.

$T_{1}$, also known as the transverse relaxation time or relaxation time,
measures the attenuation of amplitude in a quantum system. It represents the
probability that an excited state $\ket{1}$ will decay to the ground state
$\ket{0}$ after time $t$, and is modeled by the function:

$\textrm{Pr}(\ket{1}\rightarrow\ket{0})=1-\exp(-t/T_{1})\;.$ (2.1)

The decay-time probability density $f_{T}(t)$ can be described by the
exponential function:

$f_{T}(t)=T_{1}\exp^{-t/T_{1}}\;,$ (2.2)

whose mean is the density parameter $\mathds{E}(T)=T_{1}$.

$T_{2}$ is a measure of how long it takes for a qubit in the superposition
state to decay. Specifically, it measures the decay of the off-diagonal
elements of the density matrix and is modeled by an exponential decay
function. Therefore, it captures the loss of synchronization between the basis
states of an arbitrary quantum ensemble. There are two types of $T_{2}$ time
often quoted in literature [20]:

  * •

Ramsey dephasing time $T_{2}^{*}$: measures the time-scale at which a quantum
register experiences de-phasing effects when left to evolve freely

  * •

Hahn-echo dephasing time $T_{2}^{\textrm{echo}}$: uses intermediate $\pi$
pulses for re-focusing to increase relaxation time.

When simulating noisy circuits, the appropriate $T_{2}$ value to use depends
on whether the physical implementation of the circuit uses Hahn-echo for noise
suppression or not. We will specifically focus on the Hahn-echo with one echo
$T_{2}$ time.

Finally, the pure dephasing time $(T_{\phi})$ is an upper bound on the
decoherence time for a qubit, since thermal fluctuations in the environment
inevitably cause a loss of phase coherence. In practice, the dominant
relaxation time is usually $T_{2}$ (or sometimes $T_{1}$), rather than
$T_{\phi}$ [20].

The three decoherence benchmarks are related by:

$\frac{1}{T_{2}}=\frac{1}{2T_{1}}+\frac{1}{T_{\phi}}\;.$ (2.3)

####  2.3.1 Experimental characterization

We analyzed decoherence (i.e. $T_{1}$ and $T_{2}$ times) in the transmon
processor ibm_kolkata. We had 24-hour access on Tuesday, September 12, to
September 13 (from 12 noon to 12 noon) through OLCF. We chose this time-frame
as it is typical for user program queues for execution on the IBM platform. We
measured fluctuations in $T_{1}$ and $T_{2}$ times for all 27 qubits on the
device during this 24-hour period. We validated our software through numerical
simulations (detailed in the end of this section).

The complete quantum circuit used to gather the decoherence parameters $T_{1}$
and $T_{2}$ is too large to display in its entirety. However, in Fig. 2.1, we
provide a concise representation of a section of the circuit, specifically for
qubit $\ket{0}$. It’s important to clarify that this circuit structure is
replicated for all $27$ qubits in ibm_kolkata, and the entire circuit is
executed multiple times to obtain statistical data.

In the sub-circuit presented in Fig. 2.1, we illustrate only one mid-circuit
reset for the sake of clarity. In reality, the full circuit employs three
conditional resets to ensure a high probability of mid-circuit reset success.
This choice of three resets aligns with qiskit guidelines to optimize the
likelihood of successful resets.

It’s worth noting that mid-circuit measurement allows for the simultaneous
collection of the decoherence parameters $T_{1}$ and $T_{2}$ with a time
interval of just a few hundred microseconds. This simultaneous data collection
facilitates the empirical calculation of temporal correlations between these
parameters.

The basic $T_{1}$ measurement circuit begins by initializing a qubit to the
$\ket{0}$ state and then applying an X gate to transition it to the $\ket{1}$
state. Subsequently, a phase gate is introduced, during which the qubit is
affected by noise. Following this, a measurement is performed in the Z basis.
In the absence of noise, the measurement would yield the $\ket{1}$ state with
complete certainty. However, in the presence of noise, the probability of
obtaining the $\ket{1}$ state is less than 100%, and this probability depends
solely on $T_{1}$. Therefore, by analyzing the observed probability of
measuring $\ket{1}$, we can deduce an estimate for the $T_{1}$ time.

The basic $T_{2}$ circuit starts by setting a qubit to the $\ket{0}$ state. It
then uses the Hadamard gate to create an equal superposition of $\ket{0}$ and
$\ket{1}$. After a brief phase-shift delay, another H gate is applied,
followed by a measurement in the Z basis. In the absence of noise, the circuit
guarantees a 100 percent chance of retrieving the $\ket{0}$ state. However, if
there’s dephasing noise, the probability decreases, and this can be used to
estimate the $T_{2}$ time.

It’s important to note that these calculations also consider SPAM noise, which
will be discussed in detail later.

To calculate $T_{1}$ and $T_{2}$, data is fitted to an exponential decay plot
using four different evolution times: 10 $\mu$s, 50 $\mu$s, 100 $\mu$s, and
160 $\mu$s. For each qubit, the process involves the following sequence:

  1. 1.

basic $T_{1}$ circuit for a 10 $\mu$s, followed by measurement and reset

  2. 2.

basic $T_{2}$ circuit for a 10 $\mu$s, followed by measurement and reset

  3. 3.

basic $T_{1}$ circuit for a 50 $\mu$s, followed by measurement and reset

  4. 4.

basic $T_{2}$ circuit for a 50 $\mu$s, followed by measurement and reset

  5. 5.

basic $T_{1}$ circuit for a 100 $\mu$s, followed by measurement and reset

  6. 6.

basic $T_{2}$ circuit for a 100 $\mu$s, followed by measurement and reset

  7. 7.

basic $T_{1}$ circuit for a 160 $\mu$s, followed by measurement and reset

  8. 8.

basic $T_{2}$ circuit for a 160 $\mu$s, followed by measurement

The statistical estimation of $T_{1}$ and $T_{2}$ are impacted by SPAM noise.
We model the SPAM noise using a binary asymmetric model with the parameters:

$\displaystyle u_{00}+u_{01}$ $\displaystyle=1$ (2.4) $\displaystyle
u_{10}+u_{11}$ $\displaystyle=1$

where, $u_{00}$ denotes the probability of getting 0 given an input of 0,
$u_{01}$ denotes the probability of getting 1 given an input of 0, $u_{10}$
denotes the probability of getting 0 given an input of 1, and $u_{11}$ denotes
the probability of getting 1 given an input of 1. We define $p$ as the
survival probability for the $T_{1}$ circuit, indicating the probability of an
excited state enduring beyond time $t$ in $T_{1}$ measurement. Similarly, $q$
denotes the survival probability for the $T_{2}$ circuit, reflecting the
probability of observing a ground state after time $t$ in $T_{2}$ measurement.

$\displaystyle\tilde{p}_{1}$ $\displaystyle=\frac{\\#1^{\prime}s\text{
observed post-measurement of $T_{1}$ circuit}}{{\\#\text{Circuit
repetitions}}}$ (2.5) $\displaystyle\tilde{q}_{0}$
$\displaystyle=\frac{\\#0^{\prime}s\text{ observed post-measurement of $T_{2}$
circuit}}{{\\#\text{Circuit repetitions}}}$

The survival probabilities in absence of SPAM error are given by:

$\displaystyle p_{1}$ $\displaystyle=e^{-\tau/T_{1}}$ (2.6) $\displaystyle
q_{0}$ $\displaystyle=\frac{1}{2}(1+e^{-\tau/T_{2}})$

While the formula for $p_{1}$ (for $T_{1}$) is straightforward, the derivation
for $q_{0}$ (for Hahn-echo $T_{2}$ with one echo) is a little more involved.
The steps are as follows. First we initialize the qubit in the ground state:

$\rho_{0}=\left(\begin{array}[]{ll}1&0\\\ 0&0\end{array}\right)$ (2.7)

Then, we subject it to a Hadamard gate:

$\rho_{1}=H\rho_{0}H^{\dagger}=\frac{1}{2}\left(\begin{array}[]{ll}1&1\\\
1&1\end{array}\right)$ (2.8)

Then, we evolve the density matrix for time $\tau/2$ using the Hamiltonian
$\Omega=\lambda\ket{1}\bra{1}$. The unitary operator for this phase gate $D$
is:

$D=\exp(-i\Omega\frac{\tau}{2})=\exp(-i\lambda\frac{\tau}{2})\ket{1}\bra{1}+\ket{0}\bra{0}\\\
\equiv\left(\begin{array}[]{ll}1&0\\\
0&e^{-i\lambda\frac{\tau}{2}}\end{array}\right)$ (2.9)

Thus, the state becomes:

$\rho_{2}=D\rho_{1}D^{\dagger}=\frac{1}{2}\left(\begin{array}[]{cc}1&e^{-i\lambda\frac{\tau}{2}}\\\
e^{i\lambda\frac{\tau}{2}}&1\end{array}\right)$ (2.10)

However, during this time-evolution, it is acted upon by an amplitude and
phase damping (APD) channel. The state after taking APD noise into account is:

$\rho_{3}=\mathcal{E}_{\text{APD}}\left(\rho_{2}\right)=\left(\begin{array}[]{cc}1-\frac{1}{2}e^{-\frac{\tau}{2}/T_{1}}&\frac{1}{2}e^{-i\lambda
t}e^{-\frac{\tau}{2}/T_{2}}\\\ \frac{1}{2}e^{i\lambda
t}e^{-\frac{\tau}{2}/T_{2}}&\frac{1}{2}e^{-\frac{\tau}{2}/T_{1}}\end{array}\right)$
(2.11)

After this comes a deliberate bit-flip through a $\mathds{X}$ gate (assumed
noiseless):

$\rho_{4}=X\rho_{3}X=\left(\begin{array}[]{cc}\frac{1}{2}e^{-\frac{\frac{\tau}{2}}{T_{1}}}&\frac{1}{2}e^{i\lambda\frac{\tau}{2}}e^{-\frac{\tau}{2}/T_{2}}\\\
\frac{1}{2}e^{-i\lambda\frac{\tau}{2}}e^{-\frac{\tau}{2}/T_{2}}&1-\frac{1}{2}e^{-\frac{\tau}{2}/T_{1}}\end{array}\right)$
(2.12)

This is followed up with another phase gate $D$ subject to APD noise:

$\rho_{5}=D\rho_{4}D^{\dagger}=\left(\begin{array}[]{cc}\frac{1}{2}e^{-\frac{\tau}{2}/T_{1}}&\frac{1}{2}e^{-\frac{\frac{\tau}{2}}{T_{2}}}\\\
\frac{1}{2}e^{-\frac{\tau}{2}/T_{2}}&1-\frac{1}{2}e^{-\frac{\frac{\tau}{2}}{T_{1}}}\end{array}\right)$
(2.13)
$\rho_{6}=\mathcal{E}_{\text{APD}}(\rho_{5})=\left(\begin{array}[]{cc}1-e^{-\frac{\tau}{2}/T_{1}}\left(1-\frac{1}{2}e^{-\frac{\tau}{2}/T_{1}}\right)&\frac{1}{2}e^{-\frac{2\frac{\tau}{2}}{T_{2}}}\\\
\frac{1}{2}e^{-2\frac{\tau}{2}/T_{2}}&e^{-\frac{\tau}{2}/T_{1}}\left(1-\frac{1}{2}e^{-\frac{\tau}{2}/T_{1}}\right)\end{array}\right)$
(2.14)

Then another Hadamard is applied:

$\rho_{7}=\mathds{H}\rho\mathds{H}^{\dagger}$ (2.15)

Finally, we measure the probability of the qubit being in the ground state.
The probability of observing $\ket{0}$ in the final measurement is:

$\displaystyle q_{0}=\text{Pr}(0)$
$\displaystyle=\frac{1}{2}\left(1+e^{-\tau/T_{2}}\right)$ (2.16)

The survival probabilities in presence of SPAM error are given by:

$\displaystyle\tilde{p}_{1}$
$\displaystyle=p_{1}u_{11}+p_{0}u_{01}=e^{-\tau/T_{1}}\left[1-u_{10}-u_{01}\right]+u_{01}$
(2.17) $\displaystyle\tilde{q}_{0}$
$\displaystyle=q_{0}u_{00}+q_{1}u_{10}=e^{-\tau/T_{2}}\frac{1-u_{10}-u_{01}}{2}+\frac{1+u_{10}-u_{01}}{2}$

####  2.3.2 $T_{1}$ estimation

For the $T_{1}$ circuit, the observed data $y_{l}$ measured in the Z-basis
follows a Bernoulli distribution:

$\displaystyle y_{l}$ $\displaystyle=\begin{cases}1,\;\;\;\;\text{ with
probability }\tilde{p}_{1}\\\ 0,\;\;\;\;\text{ with probability
}1-\tilde{p}_{1}\\\ \end{cases}$ (2.18)

where $\tilde{p}_{1}=e^{-\tau/T_{1}}\left[1-u_{10}-u_{01}\right]+u_{01}$. The
Likelihood function is given by:

$\displaystyle\mathcal{L}=$
$\displaystyle\prod_{l=1}^{L}\text{Pr}\left(Y_{l}=y_{l}\right)=\prod_{l=1}^{L}\tilde{p}_{1}^{y_{l}}(1-\tilde{p}_{1})^{1-y_{l}}$
(2.19)

Setting:

$\frac{\partial\log\mathcal{L}}{\partial\tilde{p}_{1}}=0,$ (2.20)

we get,

$\displaystyle\tilde{p}_{1}^{*}=\sum\frac{y_{l}}{L}$ (2.21)

Thus we have obtained the estimator (denoted by the * sign) for
$\tilde{p}_{1}$. Note that

$\mathds{E}(\tilde{p}_{1}^{*})=\sum\frac{\mathds{E}(y_{l})}{L}=\sum\frac{\tilde{p}_{1}}{L}=\tilde{p}_{1}$
(2.22)

and hence it is an unbiased estimator.

The variance of this estimate is given by:

$\displaystyle\text{Var}(\tilde{p}_{1}^{*})=\sigma^{2}(\tilde{p}_{1}^{*})=\sum\frac{\text{Var}(y_{l})}{L^{2}}=\frac{\tilde{p}_{1}(1-\tilde{p}_{1})}{L}\approx\frac{\tilde{p}_{1}^{*}(1-\tilde{p}_{1}^{*})}{L}$
(2.23)

Using the above, we can get K different equations, one for each phase gate
with evolution time $\tau_{k}$ and whose post-measurement results are denoted
by $\\{y_{l,k}\\}$:

$\displaystyle e^{-\tau_{1}/T_{1}}\left[1-u_{10}-u_{01}\right]+u_{01}=$
$\displaystyle\frac{\sum y_{l,1}}{L}$ (2.24) $\displaystyle
e^{-\tau_{2}/T_{1}}\left[1-u_{10}-u_{01}\right]+u_{01}=$
$\displaystyle\frac{\sum y_{l,2}}{L}$ $\displaystyle\cdots$ $\displaystyle
e^{-\tau_{K}/T_{1}}\left[1-u_{10}-u_{01}\right]+u_{01}=$
$\displaystyle\frac{\sum y_{l,K}}{L}$

Each of these K equations has three unknowns: $T_{1},u_{10}$ and $u_{01}$.
From these K equations, we find the best fit value for $T_{1}$ using the
scipy.optimize.minimize module in python. K has to be at least 3 so that the
problem is not underspecified (overspecified is okay).

####  2.3.3 $T_{2}$ estimation

In an analogous manner, for the $T_{2}$ circuit, the observed data $y_{l}$
measured in the Z-basis follows a Bernoulli distribution:

$\displaystyle y_{l}$ $\displaystyle=\begin{cases}1,\;\;\;\;\text{ with
probability }1-\tilde{q}_{0}\\\ 0,\;\;\;\;\text{ with probability
}\tilde{q}_{0}\\\ \end{cases}$ (2.25)

where
$\tilde{q}_{0}=e^{-\tau/T_{2}}\frac{\left[1-u_{10}-u_{01}\right]}{2}+\frac{\left[1+u_{10}-u_{01}\right]}{2}$.

The Likelihood function is given by:

$\displaystyle\mathcal{L}=$
$\displaystyle\prod_{l=1}^{L}\text{Pr}\left(Y_{l}=y_{l}\right)=\prod_{l=1}^{L}\tilde{q}_{0}^{1-y_{l}}(1-\tilde{q}_{0})^{y_{l}}$
(2.26)

Setting:

$\frac{\partial\log\mathcal{L}}{\partial\tilde{q}_{0}}=0,$ (2.27)

we get,

$\displaystyle\tilde{q}_{0}^{*}=1-\sum\frac{y_{l}}{L}$ (2.28)

Thus we have obtained the estimator (denoted by the * sign) for
$\tilde{q}_{0}$. Note that

$\mathds{E}(\tilde{q}_{0}^{*})=\sum\frac{\mathds{E}(y_{l})}{L}=\sum\frac{\tilde{q}_{0}}{L}=\tilde{q}_{0}$
(2.29)

and hence it is an unbiased estimator.

The variance of this estimate is given by:

$\displaystyle\text{Var}(\tilde{q}_{0}^{*})=\sigma^{2}(\tilde{q}_{0}^{*})=\sum\frac{\text{Var}(y_{l})}{L^{2}}=\frac{\tilde{q}_{0}(1-\tilde{q}_{0})}{L}\approx\frac{\tilde{q}_{0}^{*}(1-\tilde{q}_{0}^{*})}{L}$
(2.30)

Using the above, we can get K different equations, one for each phase gate
with evolution time $\tau_{k}$ and whose post-measurement results are denoted
by $\\{y_{l,k}\\}$:

$\displaystyle
e^{-\tau_{1}/T_{2}}\frac{1-u_{10}-u_{01}}{2}+\frac{1+u_{10}-u_{01}}{2}=$
$\displaystyle\frac{\sum y_{l,1}}{L}$ (2.31) $\displaystyle
e^{-\tau_{2}/T_{2}}\frac{1-u_{10}-u_{01}}{2}+\frac{1+u_{10}-u_{01}}{2}=$
$\displaystyle\frac{\sum y_{l,2}}{L}$ $\displaystyle\cdots$ $\displaystyle
e^{-\tau_{K}/T_{2}}\frac{1-u_{10}-u_{01}}{2}+\frac{1+u_{10}-u_{01}}{2}=$
$\displaystyle\frac{\sum y_{l,K}}{L}$

Each of these K equations has three unknowns: $T_{2},u_{10}$ and $u_{01}$.
From these K equations, we find the best fit value for $T_{2}$ using the
scipy.optimize.minimize module in python. K has to be at least 3 so that the
problem is not underspecified (overspecified is okay).

####  2.3.4 Error bars on decoherence estimates

For the $T_{1}$ circuit:

$\displaystyle p_{1}=$ $\displaystyle e^{-\tau/T_{1}}$ (2.32)
$\displaystyle\Rightarrow\log p_{1}=$ $\displaystyle-\frac{\tau}{T_{1}}$
$\displaystyle\Rightarrow\frac{\delta p_{1}}{p_{1}}=$
$\displaystyle\frac{\tau}{T_{1}^{2}}\delta T_{1}$
$\displaystyle\Rightarrow\sigma^{2}_{T_{1}}=$
$\displaystyle\frac{T_{1}^{4}\left(e^{\tau/T_{1}}-1\right)}{L\tau^{2}}$

Now suppose that we use $K$ different evolution times in the circuit:
$\tau_{1},\tau_{2},\cdots\tau_{K}$. Let the desired standard deviation of the
$T_{1}$ estimate be $\sigma_{\text{desired}}$. Since the data underlying the
estimation obtained at different times are independent:

$\displaystyle\sigma_{\text{desired}}^{2}=\frac{1}{K}\sum\sigma_{i}^{2}$
(2.33)

where $\sigma_{i}^{2}$ is the variance obtained when $T_{1}$ is estimated
using a delay gate with delay time $\tau_{i}$. This gives:

$\displaystyle
L_{\text{min}}=\frac{T_{1}^{4}}{K\sigma^{2}_{\text{desired}}}\sum\frac{e^{\tau_{i}/T_{1}}-1}{\tau_{i}^{2}}$
(2.34)

For the $T_{2}$ circuit:

$\displaystyle q_{0}=$ $\displaystyle\frac{1}{2}(1+e^{-\tau/T_{2}})$ (2.35)
$\displaystyle\Rightarrow\log(2q_{0}-1)=$ $\displaystyle-\frac{\tau}{T_{2}}$
$\displaystyle\Rightarrow\frac{2\delta q_{0}}{2q_{0}-1}=$
$\displaystyle\frac{\tau}{T_{2}^{2}}\delta T_{2}$
$\displaystyle\Rightarrow\sigma^{2}_{T_{2}}=$
$\displaystyle\frac{T_{2}^{4}\left(e^{2\tau/T_{2}}-1\right)}{L\tau^{2}}$

Now suppose that we use $K$ different evolution times in the circuit:
$\tau_{1},\tau_{2},\cdots\tau_{K}$. Let the desired standard deviation of the
$T_{2}$ estimate be $\sigma_{\text{desired}}$. Since the data underlying the
estimation obtained at different times are independent:

$\displaystyle\sigma_{\text{desired}}^{2}=\frac{1}{K}\sum\sigma_{i}^{2}$
(2.36)

where $\sigma_{i}^{2}$ is the variance obtained when $T_{2}$ is estimated
using phase-shift time $\tau_{i}$. This gives:

$\displaystyle
L_{\text{min}}=\frac{T_{2}^{4}}{K\sigma^{2}_{\text{desired}}}\sum\frac{e^{2\tau_{i}/T_{2}}-1}{\tau_{i}^{2}}$
(2.37)

Since we have a circuit that measures $T_{1}$ and $T_{2}$ in one go, we have
to take the max of the $L_{\text{min}}$ for each of the two cases (i.e. max of
$L_{\text{min}}$ for $T_{1}$ and $L_{\text{min}}$ for $T_{2}$).

Using a noisy simulation with known error estimates (detailed in the program
validation section next), we arrived at $L_{\text{min}}$ at each time to be
$20,000$ samples. There were four delay gates (aka four separate decay
experiments), so a total of $80,000$ samples went into the computation of each
$T_{1}$ and $T_{2}$ value for each time-stamp. The standard-deviation we aimed
for is $1\mu s$.

####  2.3.5 Program validation

Program validation is crucial to ensure that, within the assumptions of our
theoretical noise model, our simulations of quantum circuits, utilizing known
noise parameters, yield precise results. Specifically, we aimed for our
statistical analysis to accurately recover the expected values. Additionally,
program validation helps us determine the required number of circuit
repetitions for achieving a specified level of outcome precision. While minor
deviations in the final standard deviation are expected due to various
factors, they should generally align with our target precision.

Since we individually measure the $T_{1}$ and $T_{2}$ times for each of the 27
qubits of ibm_kolkata, validating the program for a single qubit suffices. We
found that a sample size of 10,000 is sufficient to attain the desired
precision of 1 microsecond, but for safety, we opted for 20,000 samples in the
final run on real device. Given the use of four delay gates, each $T_{1}$ or
$T_{2}$ data point estimation relies on a total of $80,000$ samples. For
consistency with the ibm_kolkata processor’s specifications, we set the
readout error at 0.028 $T_{1}$ time at 134 microseconds and $T_{2}$ time at 93
microseconds in our noise simulation.

Our program validation yielded an estimated T1 time of 134.52 microseconds,
which falls within the 1-microsecond precision target. The estimated T2 time
was 94.05 microseconds, slightly exceeding the 1-microsecond precision by a
difference of 1.03 microseconds which we deemed acceptable.

####  2.3.6 Summary of results

In terms of the $T_{1}$ parameter, qubit 15 performed the best, exhibiting a
median $T_{1}$ time of 184 microseconds, while qubit 4 performed the worst
with a median $T_{1}$ time of 74 microseconds. Across all 27 qubits, the
median $T_{1}$ time was 116 microseconds, with a standard deviation of 1.5
microseconds, and the range of $T_{1}$ times spanned from 32 to 297
microseconds (see Fig. 2.2 (a)). As for the $T_{2}$ parameter, qubit 3
demonstrated the highest performance, displaying a median $T_{2}$ time of 76
microseconds, while qubit 19 had the poorest performance with a median $T_{2}$
time of 15 microseconds. Across all 27 qubits, the median $T_{2}$ time
averaged 30 microseconds, with a standard deviation of 1.1 microseconds, and
the range of $T_{2}$ times varied from 3 to 191 microseconds (see Fig. 2.2
(b)). The time-series for all the 27 qubits over the 24 hour period is shown
in Figs. 2.4 (b)- 2.17.

###  2.4 Modeling quantum noise channels

####  2.4.1 Amplitude and Phase damping channel

The amplitude damping channel $\mathcal{E}^{\textrm{AD}}(\cdot)$ and de-
phasing channel $\mathcal{E}^{\textrm{PD}}(\cdot)$ are two fundamental sources
of quantum de-coherence and information loss in transmons [18, 36, 11]. A
realistic model for this noise channel, denoted as APD, involves a combination
of amplitude damping and de-phasing. Amplitude damping can be described by the
Kraus operators $E_{0}^{\textrm{AD}}$ and $E_{1}^{\textrm{AD}}$, while phase
damping can be described by $E_{0}^{\textrm{PD}}$ and $E_{1}^{\textrm{PD}}$,
as follows [11]:

$\mathcal{E}^{\textrm{AD}}(\rho)=\sum\limits_{k=0}^{1}E_{k}^{\textrm{AD}}\rho
E_{k}^{\textrm{AD}\dagger}\;,$ (2.38)
$E_{0}^{\textrm{AD}}=\begin{pmatrix}1&0\\\ 0&\sqrt{1-\gamma}\end{pmatrix}\;,$
(2.39) $E_{1}^{\textrm{AD}}=\begin{pmatrix}0&\sqrt{\gamma}\\\
0&0\end{pmatrix}\;,$ (2.40)
$\mathcal{E}^{\textrm{PD}}(\rho)=\sum\limits_{k=0}^{1}E_{k}^{\textrm{PD}}\rho
E_{k}^{\textrm{PD}\dagger}\;,$ (2.41)
$E_{0}^{\textrm{PD}}=\begin{pmatrix}1&0\\\ 0&\sqrt{1-\lambda}\end{pmatrix}\;,$
(2.42) $E_{1}^{\textrm{PD}}=\begin{pmatrix}0&0\\\
0&\sqrt{\lambda}\end{pmatrix}\;.$ (2.43)

Here, $\gamma=1-\exp(-t/T_{1})$ and $\lambda=1-\exp(-t/T_{2})$, where $t$ is
the time scale of the decoherence process. The relation between $T_{1}$,
$T_{\phi}$, and $T_{2}$ was previously discussed in Eq. (2.3). The Kraus
decomposition of the combined amplitude and phase damping channel
$\mathcal{E}^{\textrm{APD}}(\cdot)$, valid for a single qubit, can be
expressed as $E_{0}^{\textrm{APD}}$, $E_{1}^{\textrm{APD}}$, and
$E_{2}^{\textrm{APD}}$.

$\mathcal{E}^{\textrm{APD}}(\rho)\coloneqq\mathcal{E}^{\textrm{PD}}\circ\mathcal{E}^{\textrm{AD}}=\sum\limits_{k=0}^{2}E_{k}^{\textrm{APD}}\rho
E_{k}^{\textrm{APD}^{\dagger}}\;,$ (2.44)

where,

$E_{0}^{\textrm{APD}}=E_{0}^{\textrm{PD}}E_{0}^{\textrm{AD}}=\begin{pmatrix}1&0\\\
0&\sqrt{[1-\gamma][1-\lambda]}\end{pmatrix}\;,$ (2.45)
$E_{1}^{\textrm{APD}}=E_{0}^{\textrm{PD}}E_{1}^{\textrm{AD}}=\begin{pmatrix}0&\sqrt{\gamma}\\\
0&0\end{pmatrix}\;,$ (2.46)
$E_{2}^{\textrm{APD}}=E_{1}^{\textrm{PD}}E_{0}^{\textrm{AD}}=\begin{pmatrix}0&0\\\
0&\sqrt{[1-\gamma]\lambda}\end{pmatrix}\;.$ (2.47)

Using the fact that:

$\displaystyle
E_{0}^{\textrm{APD}}=\frac{1+\sqrt{1-\lambda-\gamma+\lambda\gamma}}{2}I+\frac{1-\sqrt{1-\lambda-\gamma+\lambda\gamma}}{2}Z\;,$
(2.48) $\displaystyle
E_{1}^{\textrm{APD}}=\frac{\sqrt{\gamma}}{2}X+\frac{\sqrt{\gamma}}{2}iY\;,$
(2.49) $\displaystyle
E_{2}^{\textrm{APD}}=\frac{\sqrt{\lambda-\lambda\gamma}}{2}I-\frac{\sqrt{\lambda-\lambda\gamma}}{2}Z\;.$
(2.50)

the APD channel can be expressed as:

$\begin{split}\mathcal{E}^{\textrm{APD}}(\rho)=&\frac{2-\gamma+2\sqrt{1-\lambda-\gamma+\lambda\gamma}}{4}\rho+\frac{\gamma}{4}X\rho
X-\frac{\gamma}{4}Y\rho Y\\\
&+\frac{2-\gamma-2\sqrt{1-\lambda-\gamma+\lambda\gamma}}{4}Z\rho Z\\\
&+\frac{\gamma}{4}I\rho Z+\frac{\gamma}{4}Z\rho I-\frac{\gamma}{4i}X\rho
Y-\frac{\gamma}{4i}Y\rho X\;,\end{split}$ (2.51)

where $\lambda$ and $\gamma$ are the APD parameters, and $I$, $X$, $Y$, and
$Z$ are the Pauli matrices.

####  2.4.2 Depolarizing channel

Next, the depolarizing channel is a common type of quantum noise channel. It
works by randomly applying one of the Pauli operators
$(\mathds{X},\mathds{Y},\mathds{Z})$ to the quantum state with a certain (but
equal) probability, causing a loss of information about the state. For a
qubit:

$\mathcal{E}_{D}(\rho)=(1-\textrm{x})\rho+\textrm{x}\frac{\mathds{I}}{d},$
(2.52)

where $\rho$ is the input quantum state, x is the probability of noise
occurring, $\mathds{I}$ is the identity operator, and $d$ is the dimension of
the Hilbert space ($d=2$ for a qubit).

####  2.4.3 Pauli noise channel

The Pauli noise channel is a generalization of the depolarizing channel. It
encompasses the effect of the bit-flip ($\mathds{X}$), phase-flip
($\mathds{Z}$), and bit-phase-flip ($\mathds{Y}$) errors with unequal
probabilities.

The impact of Pauli noise on quantum information encoded in an $n$-qubit
register is shown below:

$\mathcal{E}_{\textrm{x}}(\rho)=\sum\limits_{i=0}^{N_{p}-1}\textrm{x}_{i}P_{i}(n)\rho
P_{i}(n)^{\dagger}$ (2.53)

where $N_{p}$ denotes the total number of Pauli coefficients and $P_{i}(n)$
represents $n$-qubit Pauli operators. The coefficients contribute to a
simplex:

$\sum\limits_{i=0}^{N_{p}-1}\textrm{x}_{i}=1,\;\;\;\;\textrm{x}_{i}\geq 0$
(2.54)

Pauli noise channel is widely used in quantum error correction because it is a
simple and natural model for random quantum noise [83, 84, 49]. It is a well-
understood and easily implementable noise model that can simulate a variety of
realistic physical processes that lead to quantum errors, such as dephasing,
amplitude damping, and phase-flip errors. Additionally, the Pauli noise
channel is mathematically tractable and can be efficiently simulated, making
it a useful tool for developing and testing quantum error correction
protocols. It can be used to estimate the average fidelity of a quantum gate
subject to the original APD channel and identify codes that work for the APD
channel [84]. The Pauli noise channel, although not a completely general noise
model, still manages to model many practical situations. It is widely used
because of two reasons: (a) it is efficiently simulatable on a classical
computer (per the Gottesman-Knill theorem) and (b) when used as a proxy for
physically accurate noise models (such as the amplitude and phase damping
noise) which are not efficiently simulatable on a classical computer, it still
manages to preserve important properties like entanglement fidelity [85].

Remarkably, Pauli twirling can map[86, 87, 88] a more complex quantum noise
channel (e.g. APD) to a simple Pauli channel while preserving certain features
such as the average channel fidelity and the entanglement fidelity [85].
Consider a single-qubit amplitude and phase damping channel (APD) [58]. Upon
Pauli twirling [89]:

$\displaystyle\mathcal{E}_{\textrm{twirl}}(\rho)=$
$\displaystyle\frac{1}{4}\sum\limits_{A\in\\{I,X,Y,Z\\}}A^{\dagger}\mathcal{E}_{\text{APD}}\left(A\rho
A^{\dagger}\right)A$ (2.55) $\displaystyle=$
$\displaystyle\sum\limits_{k=0}^{3}c_{k}\sigma_{k}\rho\sigma_{k}$ (2.56)

an APD channel becomes a Pauli noise channel. Here,
$\\{\sigma_{k}\\}_{k=0}^{3}=\\{I,X,Y,Z\\}$ are the Pauli matrices. Thus, the
coefficients of the Pauli noise channel are functions of the coefficients of
the original APD channel, which in turn are functions of the decoherence times
$T_{1}$ and $T_{2}$ [46]:

$\displaystyle c_{1}=c_{2}=$
$\displaystyle\frac{1}{4}\left[1-\exp\left(-t/T_{1}\right)\right]$ (2.57)
$\displaystyle c_{3}=$
$\displaystyle\frac{1}{4}\left[1-\exp\left(-t/T_{2}\right)\right]$ (2.58)
$\displaystyle c_{0}=$ $\displaystyle 1-(c_{1}+c_{2}+c_{3})$ (2.59)

This directly links the estimation of Pauli channels to the decoherence data
collected in Sec. 2.3.

Figure 2.1: Interleaved sub-circuit post transpilation.

(a)

(b) Figure 2.2:  Spatial non-stationarity of decoherence times. Indivudal
error-bars do not capture the variation across qubits. (a) In terms of the
$T_{1}$ parameter, qubit 15 performed the best, exhibiting a median $T_{1}$
time of 184 microseconds, while qubit 4 performed the worst with a median
$T_{1}$ time of 74 microseconds. Across all 27 qubits, the median $T_{1}$ time
was 116 microseconds, with a standard deviation of 1.5 microseconds, and the
range of $T_{1}$ times spanned from 32 to 297 microseconds. (b) For the
$T_{2}$ parameter, qubit 3 demonstrated the highest performance, displaying a
median $T_{2}$ time of 76 microseconds, while qubit 19 had the poorest
performance with a median $T_{2}$ time of 15 microseconds. Across all 27
qubits, the median $T_{2}$ time averaged 30 microseconds, with a standard
deviation of 1.1 microseconds, and the range of $T_{2}$ times varied from 3 to
191 microseconds.

Figure 2.3: This figure presents non-stationary temporal dynamics of
decoherence for a qubit on IBM’s belem device. The top figure shows $T_{1}$
relaxation time series for qubit 0, where two datasets were collected for 5 ms
each on the same day, separated by a vertical line. The blue dataset varies
between 116-126$\mu$s with a mean of 122 $\mu$s and a standard deviation of 2
$\mu$s, while the green dataset varies between 104-114$\mu$s with a mean of
108 $\mu$s and a standard deviation of 2 $\mu$s. The bottom figure displays
Ramsey dephasing time ($T_{2}$) series for qubit 0, where two datasets were
collected for 5 ms each on the same day, separated by a vertical line. The
blue dataset varies between 67-73$\mu$s with a mean of 70 $\mu$s and a
standard deviation of 1 $\mu$s, while the green dataset varies between
60-64$\mu$s with a mean of 62 $\mu$s and a standard deviation of 1 $\mu$s,
collected around different times on the same day. The data shows significant
non-stationarity in decoherence values over a 30-minute interval.

(a)

(b) Figure 2.4:  (a) Schematic of the 27-qubit device ibm_kolkata. (b)
Estimated $T_{1}$ and $T_{2}$ time-series for qubit 0 as collected between
12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023.

(a)

(b) Figure 2.5:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 1 and (b) qubit 2.

(a)

(b) Figure 2.6:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 3 and (b) qubit 4.

(a)

(b) Figure 2.7:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 5 and (b) qubit 6.

(a)

(b) Figure 2.8:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 7 and (b) qubit 8.

(a)

(b) Figure 2.9:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 9 and (b) qubit 10.

(a)

(b) Figure 2.10:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 11 and (b) qubit 12.

(a)

(b) Figure 2.11:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 13 and (b) qubit 14.

(a)

(b) Figure 2.12:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 15 and (b) qubit 16.

(a)

(b) Figure 2.13:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 17 and (b) qubit 18.

(a)

(b) Figure 2.14:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 19 and (b) qubit 20.

(a)

(b) Figure 2.15:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 21 and (b) qubit 22.

(a)

(b) Figure 2.16:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 23 and (b) qubit 24.

(a)

(b) Figure 2.17:  Estimated $T_{1}$ and $T_{2}$ time-series as collected
between 12:00 P.M. ET on Sep 12, 2023 and 12:00 P.M. ET on Sep 13, 2023 for
(a) qubit 25 and (b) qubit 26.

##  Chapter 3 Performance evaluation framework

Performance evaluation of noisy quantum computing in crucial for several
reasons. Firstly, quantum computing is still in its early stages of
development [55, 56, 57], and understanding the sources of errors and
noise[58] is vital to improve the performance of quantum computers. By
conducting rigorous performance evaluation, researchers can identify and
quantify the various sources of noise, such as decoherence, gate errors, and
readout errors. Secondly, this understanding is essential for developing error
mitigation techniques [59], which are necessary for scaling up small-scale
quantum computations. Thirdly, reproducibility of results from quantum
computing is critical. Rigorous performance evaluation ensure that experiments
can be replicated by other researchers, contributing to validation and
verification of quantum algorithms.

However, the task is not simple as the same complexity that gives quantum
technology an advantage over classical computing also hinders its rigorous
checking [54]. Reasons include the inherent randomness in quantum results due
to the Born rule, error accumulation without clear source attribution, the
curse of dimensionality [52], and the inability to step-through program
execution in quantum circuits [53].

Also, the diverse range of terms encountered in quantum computing today can
blur the distinctions between them, making it challenging to appreciate their
nuanced differences. Examples include verification[60, 61, 62] (which pertains
to ensuring correct transpilation), validation[63, 64] (which can have two
connotations: (a) validating correctness of output by comparison to theory or
concurrent classical simulation (akin to accuracy), or validating the quantum
nature of a device), benchmarking[65, 66, 67] (which involves assigning a
performance measure, often a simple scalar number, to a quantum processor/
subsystem/ subroutine, with reproducibility as a key defining characteristic),
accreditation [90, 68] and certification [54].

###  3.1 Distance measures

In this chapter, our focus is on the development of a systematic performance
evaluation framework for the outcomes from noisy quantum computers. To
accomplish this, it is necessary to compare probability distributions. Various
options exist for quantifying the distance between high-dimensional densities,
including Hellinger distance[91], total variation distance[92], Kolmogorov-
Smirnov statistic[93], Mahalanobis distance[94], Wasserstein metric[95], Levy-
Prokhorov metric[96], and non-metric divergence measures like Kullback-Leibler
divergence[97], Jensen-Shannon divergence[98], Renyi’s divergence[97] and
Tsallis divergence[99] for quantum states. Metric measures are advantageous
because they obey the triangle inequality, enabling rigorous comparisons,
while non-metric measures (often called divergence) are useful for obtaining
performance bounds in specific problem settings. However, all distance
measures suffer from the curse of dimensionality[52], which results in
exponentially increasing resource requirements to accurately represent
information as the number of qubits in the quantum system grows [100, 101].

####  3.1.1 Hellinger distance

In our work, we primarily employ the Hellinger distance. The Hellinger
distance between two probability distributions $f_{X}(\textrm{x})$ and
$f_{Y}(\textrm{x})$ for the random variables $X$ and $Y$ is defined by:

$H(f_{X},f_{Y})=\sqrt{1-BC(f_{X},f_{Y})},$ (3.1)

where the Bhattacharyya coefficient BC is

$BC(f_{X},f_{Y})=\int\limits_{\textrm{x}}\sqrt{f_{X}(\textrm{x})f_{Y}(\textrm{x})}d\textrm{x}.$
(3.2)

The Hellinger distance provides a practical and meaningful approach to
measuring the similarity of distributions. Firstly, it operates directly on
observed data, eliminating the need to compute intermediate abstractions like
entropy. Secondly, the Hellinger distance is easy to interpret and apply in
practical scenarios. Thirdly, it can be easily extended to quantum states
through Tsallis divergence [99]. Lastly, the Hellinger distance is
proportional to the Fisher information[102], which quantifies the partial
knowledge a density carries about some unknown.

In particular, suppose x is a realization of the random noise parameter $X$
drawn from the time-varying distribution $f_{X}(\textrm{x};t)$, where $t$
denotes time. Let,

$H_{X}(t_{1},t_{2})=H(f_{X}(t_{1}),f_{X}(t_{2})),$ (3.3)

measure the distance between the densities of $X$ at time $t_{1}$ and $t_{2}$.

For example, the Hellinger distance between two beta distributions (used for
SPAM and CNOT fidelity characterizations) is given by:

$\begin{split}f(\textrm{x};\alpha_{1},\beta_{1})&=\frac{\textrm{x}^{\alpha_{1}-1}(1-\textrm{x})^{\beta_{1}-1}}{\text{Beta}(\alpha_{1},\beta_{1})},\;\;\;\;g(\textrm{x};\alpha_{2},\beta_{2})=\frac{\textrm{x}^{\alpha_{2}-1}(1-\textrm{x})^{\beta_{2}-1}}{\text{Beta}(\alpha_{2},\beta_{2})}\\\
d_{H}&=\sqrt{1-\int\limits_{0}^{1}\sqrt{f(\textrm{x};\alpha_{1},\beta_{1})g(\textrm{x};\alpha_{2},\beta_{2})}d\textrm{x}}=\sqrt{1-\frac{\text{Beta}(\alpha_{1}+\alpha_{2}-1,\beta_{1}+\beta_{2}-1)}{\text{Beta}(\alpha_{1},\beta_{1})\text{Beta}(\alpha_{2},\beta_{2})}}\end{split}$
(3.4)

while that between two gamma distributions (used for duty cycle
characterization) is given by:

$\begin{split}f(\textrm{x};m_{1},\alpha_{1})&=\frac{1}{\Gamma(m_{1})}\textrm{x}^{m_{1}-1}\alpha_{1}^{m_{1}}e^{-\alpha_{1}\textrm{x}},\;\;\;\;g(\textrm{x};m_{2},\alpha_{2})=\frac{1}{\Gamma(m_{2})}\textrm{x}^{m_{2}-1}\alpha_{2}^{m_{2}}e^{-\alpha_{2}\textrm{x}}\\\
d_{H}&=\sqrt{1-\int\limits_{0}^{1}\sqrt{f(\textrm{x};n,\alpha)g(\textrm{x};m,\beta)}d\textrm{x}}=\sqrt{1-\frac{\alpha_{1}^{m_{2}}\alpha_{2}^{m_{1}}}{(\alpha_{1}+\alpha_{2})^{m_{1}+m_{2}-1}}\frac{\Gamma(m_{1}+m_{2}-1)}{\Gamma(m_{1})\Gamma(m_{2})}}\end{split}$
(3.5)

For non-standard distributions $f(x)$ and $g(x)$ where an analytical closed
form solution is not available, we numerically integrate the empirical
distributions using the standard trapezoidal method:

$\begin{split}d_{H}^{2}-1&=\int\limits_{a}^{b}\sqrt{f(x)g(x)}dx=\lim_{n\rightarrow\infty}\sum\limits_{i=0}^{n-1}\left[\frac{1}{2}\left[f(a-ih/2)g(a-ih/2)+f(a+ih/2)g(a+ih/2)\right]\right]^{1/2}h\\\
\end{split}$ (3.6)

where $h=\frac{b-a}{n}$.

Despite its ease of interpretation, the Hellinger distance scales
exponentially in the number of noise parameters. This has the effect that even
small changes in a distribution yield large changes in the distance value.
This is called the curse of dimensionality. To see this, consider $d$
independent and identically distributed noise parameters
$\\{\textrm{x}_{1},\cdots,\textrm{x}_{d}\\}$, whose marginals are given by
$f_{X_{i}}(\textrm{x};t)$. Let $h$ be the Hellinger distance between the
marginals at time $t_{1}$ and $t_{2}$. Thus,

$H_{X_{i}}(t_{1},t_{2})=h\;\;\;\;\forall i$ (3.7)

Since the parameters are independent:

$\begin{split}\log\left(1-H_{X}^{2}\right)&=-d|\log(1-h^{2})|\\\ \Rightarrow
H_{X}&=\sqrt{1-\exp\left[-d|\log(1-h^{2})|\right]}.\end{split}$ (3.8)

Thus the distance approaches $1$ quickly as the number of dimensions
increases.

A more sensitive measure can be defined using $H_{\textrm{avg}}$, defined as
the average over the distances for the $d$ univariate $(X_{k})$ marginal
distributions:

$H_{\textrm{avg}}(t_{1},t_{2})=\frac{1}{d}\sum\limits_{k=1}^{d}H_{X_{k}}(t_{1},t_{2}).$
(3.9)

When the joint distributions are time-invariant, then the marginals must also
be time-invariant, resulting in a small average value for $H_{\textrm{avg}}$.
This test is more sensitive as it mitigates the curse of dimensionality and
offers higher dispersion for improved calibration.

Another sensitive approach is to normalize the distance relative to the
dimensionality $d$ of the distribution:

$H_{\textrm{normalized}}(t_{1},t_{2})=\sqrt{1-BC^{1/d}}.$ (3.10)

We refer to this statistic as the normalized Hellinger distance (note that
although we call it distance, this statistic is not technically a metric as it
does not satisfy the triangle inequality).

####  3.1.2 Moment-Based Distance

We also developed a new disance measure as part of our research which we call
Moment-Based Distance (MBD). The key advantage of MBD is its ability to
incorporate the geometric shape of the underlying noise distribution while
still being a metric. Thus, it takes into account higher order effects like
kurtosis and skewness. Specifically, we define the moment-based metric ($d$)
between two histograms ($f$ and $g$) based on the equality of their moments.

$d(f,g)=\sum\limits_{m=0}^{\infty}S_{m}(f,g)$ (3.11)

where

$S_{m}(f,g)=\frac{1}{(m)!}\int\limits_{a}^{b}\left|\left(\frac{x}{\gamma}\right)^{m}(f(x)-g(x))\right|dx$
(3.12)

and,

$\gamma=\max(|a|,|b|)\\\ $ (3.13)

for a bounded real variable $x$. Here, $a$ and $b$ are the minimum and maximum
values of $x$ and can be derived from theoretical considerations (e.g., when
the random variable is a probability then $\gamma=1$) or from empirical
histogram data.

The moment-based-distance $d(f,g)$ satisfies the following properties:

  1. 1.

$d(f,g)\geq 0$ follows from the definition of $d$.

  2. 2.

$d(f,g)=d(g,f)$ follows from the definition of $d$.

  3. 3.

$d(f,g)=0$ iff $f(x)=g(x)$.
Proof: If $f(x)=g(x)$, then $d=0$ because $S_{m}=0$ for every $m$. Conversely,
if $d=0$, then $S_{m}=0$ for all $m$.If $S_{m}=0$, then for all $x$, the
integrand must satisfy

$\left|\left(\frac{x}{\gamma}\right)^{m}(f-g)\right|=0$

As $(x)^{m}\neq 0$ for all $x$, it must be that $|f(x)-g(x)|=0$ for all $x$
and, hence, $f(x)=g(x)$. $\tiny\blacksquare$

  4. 4.

$d(f,g)\leq d(f,h)+d(h,g)$
Proof: For every $m$,

$\begin{split}S_{m}(f,g)=&\int\limits_{a}^{b}\left|\left(\frac{x}{\gamma}\right)^{m}\frac{1}{m!}(f(x)-g(x))\right|dx\\\
=&\int\limits_{a}^{b}\bigg{|}\left(\frac{x}{\gamma}\right)^{m}\frac{1}{m!}(f(x)-h(x)\\\
&\hskip 56.9055pt+h(x)-g(x))\bigg{|}dx\\\
\leq&\int\limits_{a}^{b}\left|\left(\frac{x}{\gamma}\right)^{m}\frac{f(x)-h(x)}{m!}\right|dx\\\
&+\int\limits_{a}^{b}\left|\left(\frac{x}{\gamma}\right)^{m}\frac{h(x)-g(x)}{m!}\right|dx\\\
\leq&S_{m}(f,h)+S_{m}(h,g)\\\ \end{split}$

and whence the sum satisfies the inequality as well. $\blacksquare$

  5. 5.

The series $d=S_{0}+S_{1}+S_{2}+\cdots$ converges.
Proof: The distance $d$ converges if, after some fixed term, the ratio of each
term to the preceding term is less than some quantity $r$, which is itself
numerically less than unity. If $S_{m+1}<S_{m}$ for all $m\geq 1$, then

$\begin{split}d=&S_{0}+S_{1}+S_{1}\cdot\frac{S_{2}}{S_{1}}+S_{1}\cdot\frac{S_{2}}{S_{1}}\cdot\frac{S_{3}}{S_{2}}+\cdots\\\
<&S_{0}+S_{1}(1+r+r^{2}+r^{3}+\cdots)\\\ =&S_{0}+\frac{S_{1}}{1-r}\textrm{
since }r<1\\\ \end{split}$

To prove that $S_{m+1}<rS_{m}$ for $m\geq 1$, we proceed as follows:

$\begin{split}S_{m+1}&=\int\limits_{a}^{b}\left|\left(\frac{x}{\gamma}\right)^{m+1}\frac{1}{(m+1)!}(f(x)-g(x))\right|dx\\\
&=\int\limits_{a}^{b}\left|\frac{x}{\gamma(m+1)}\right|\left|\left(\frac{x}{\gamma}\right)^{m}\frac{f(x)-g(x)}{m!}\right|dx\\\
&\leq\left|\frac{x}{\gamma}\right|_{max}\frac{1}{m+1}S_{m}\\\
&\leq\frac{1}{m+1}S_{m}\textrm{ since
}\left|\frac{x}{\gamma}\right|_{max}=1\\\ &\leq\frac{1}{2}S_{m}\textrm{ since
}m\geq 1\end{split}$

where $|x|_{max}$ is the maximum of $x$. $\blacksquare$

An important consequence of the latter convergence property is that the
moment-based distance satisfies the practical requirement that lower-order
moments contribute more than higher-order moments to the distance (for $m>1$).
This proves essential to our subsequent use of the moment-based distance
below, as we rely on the approximate distance defined to order $n$ as

$d_{n}=\sum_{m=0}^{n}{S_{m}}$ (3.14)

We next present a series of simulation studies to develop intuition for how
the moment-based distance behaves in the presence of both stable and unstable
distributions. In particular, we will show that moment-based distance is small
but non-zero for distributions that are similar but not identical, while such
deviations grow with dissimilarity. For our studies, we computed the distance
of 10 different distributions with respect to a reference distribution. Table
3.1 summaries the list of tests as well as their moment-based distance from
the reference normal distribution $\mathcal{N}(\mu,\sigma)$.

Table 3.1: Moment-based distance by Distribution Distribution | $d_{4}$ | $d_{20}$ | Error(%)
---|---|---|---
$N(\mu,\sigma)$ | 0.00000 | 0.00000 | NA
$N(\mu+\Delta,\sigma)$ | 2.70868 | 2.70876 | -0.00289
$N(\mu,2\sigma)$ | 0.83252 | 0.83253 | -0.00104
$N(\mu,4\sigma)$ | 1.47301 | 1.47304 | -0.00180
$N(2\mu,\sigma)$ | 2.93489 | 2.93520 | -0.01033
$N(\mu,1.5\sigma)$ | 0.49215 | 0.49216 | -0.00091
$N(1.01\mu,\sigma)$ | 0.11739 | 0.11740 | -0.00079
$SkewedNormal(\mu,2\sigma)$ | 0.80887 | 0.80888 | -0.00140
$Gumbel(\mu,2\sigma)$ | 0.95131 | 0.95134 | -0.00246

For testing purpose, the parameters are $\mu=0.4,\Delta=0.2$ and
$\sigma=0.04$. We note that, as expected, the distribution ‘closest’ to
$\mathcal{N}(\mu,\sigma)$ is $\mathcal{N}(1.01\mu,\sigma)$ and the ‘farthest’
are $\mathcal{N}(2\mu,\sigma)$ and $\mathcal{N}(\mu+2\Delta,2\sigma)$.

We next study how the order of the series $d_{n}$ increases the accuracy of
the distance measured. In our simulation studies of well-defined
distributions, we find that $d_{n}$ converges for $n=4$ when the distributions
are sufficiently dissimilar. As shown in Fig. 3.5, the relative contributions
of each $S_{m}$ to $d_{n}$ decreases with increasing $m$ as expected from the
convergence property. Thus, $m=0$ accounts for about $60\%$ of the total
distance while $m=1$ accounts for $90\%$ and $m=2$ reaches $98\%$. For $m=4$,
$d_{m}$ is nearly $100\%$ of the $d_{\infty}$. Consequently, we will consider
$m=4$ sufficient to accurately characterize the moment-based distance for the
remainder of our analysis. This is certainly an approximation in the sense
that two histograms which start to differ only after the fourth order moment
will be erroneously classified as same. Is $d_{4}$ still a valid distance
metric? Yes. A glance at the proofs will reveal that properties (1) to (4) are
still satisfied when we truncate the d series at a finite m (say m=4).
Moreover, it converges too (i.e. Property (5) is satisfied too) because a
finite number of terms (in this case 5 terms) is by definition convergent when
the individual terms are finite. The latter is true because each $S_{m}$ is
bounded between finite $a$ and $b$ as per Equation (2).

As a point of comparison, we contrast the moment-based distance to total
variation distance (TVD), a state-of-the-art metric which has proven useful in
earlier experimental investigations [103, 104]. We note that the magnitudes of
the moment-based distance and total variation distance are not directly
comparable as they follow very different methodologies but one can compare the
corresponding signal-to-noise ratio (SNR) of the two metrics as the inverse of
the coefficient of variation. For our numerical studies, we generated two time
series, each of length 8192, by sampling two different probability
distributions. The first was a normal distribution with mean 10 and standard
deviation 1, and the second a normal distribution with mean 10 and standard
deviation 4. We calculate the moment-based distance and total variation
distance between these two time series, and then we repeated this numerical
experiment 400 times to generate a distribution of the TVD and MBD distances.
Using the average $\mu$ and standard deviation $\sigma$ of these
distributions, we calculated the respective SNR as

$\textrm{SNR}=\frac{\mu}{\sigma}$ (3.15)

As shown in Fig. 3.6, our results indicate that the moment-based distance has
more statistical power as indicated by a higher SNR. As an aside, a practical
concern is the dependence of precision of the moment-based distance on
sampling size. Although each $S_{m}$ should vanish when two distributions are