Identifiability and Characterization of Transmon Qutrits Through Bayesian Experimental Design

Robust control of a quantum system is essential to utilize the current noisy quantum hardware to their full potential, such as quantum algorithms. To achieve such a goal, systematic search for an optimal control for any given experiment is essential. Design of optimal control pulses require accurate numerical models, and therefore, accurate characterization of the system parameters. We present an online, Bayesian approach for quantum characterization of qutrit systems which automatically and systematically identifies the optimal experiments that provide maximum information on the system parameters, thereby greatly reducing the number of experiments that need to be performed on the quantum testbed. Unlike most characterization protocols that provide point-estimates of the parameters, the proposed approach is able to estimate their probability distribution. The applicability of the Bayesian experimental design technique was demonstrated on test problems where each experiment was defined by a parameterized control pulse. In addition to this, we also presented an approach for iterative pulse extension which is robust under uncertainties in transition frequencies and coherence times, and shot noise, despite being initialized with wide uninformative priors. Furthermore, we provide a mathematical proof of the theoretical identifiability of the model parameters and present conditions on the quantum state under which the parameters are identifiable. The proof and conditions for identifiability are presented for both closed and open quantum systems using the Schroedinger equation and the Lindblad master equation respectively.


Introduction
The Noisy Intermediate Scale Quantum (NISQ) devices are poised to demonstrate quantum advantage [1].They are, however, very sensitive to a myriad of noise sources such as from interaction with the environment, characterized by coherence times [2,3], and fluctuating quantum systems and interaction with defects [4,5], causing uncertainties in qubit transition frequencies.This instability requires quantum error correction and dynamical decoupling techniques which rely on a larger number of qubits and lengthy circuit depth to encode quantum information and cancel out the noises [6], thereby limiting the quantum algorithms that can be implemented within the device's coherence times.Hence, variational quantum algorithms (VQA) [7,8], which couple computations on the stable classical hardware with the noisy quantum processing units (QPUs), have emerged as a promising strategy to address these drawbacks.Such classical-quantum algorithms typically require accurate modeling of the quantum systems.These numerical models can be parameterized by a few quantities of the QPUs, thus, accurate determination of these parameters is paramount for high-fidelity simulation of the QPU.However, the noisy behavior of current devices poses a great challenge to identifying the testbed characteristics, particularly, providing point-estimates of these parameters.Under such circumstance, posing the characterization problem within a probabilistic Bayesian framework is ideal.Regardless of the approach used for characterization, the identifiability of the model parameters must first be established.Although crucial, such priori analysis is often neglected.Such identifiability analysis not only determine the parameters that can reliably and uniquely estimated, but can also identify the correlated parameters.Furthermore, they establish a set of conditions on the experiments/measurements that make parameters identifiable, thereby informing the design of experiments.One such approach involves characterization (i.e., parameter estimation) through online Bayesian experimental design (BExD).The BExD technique uses classical optimization techniques to identify experimental parameters (or experiments) that provide maximum information about the parameters being identified and iteratively updates the prior distribution of the parameters.Estimating complete probability distributions provide insight on the correlation between parameters and allow for design of control pulses that are robust to slight variances in parameters.Furthermore, such distributions can better inform development of noise models and has been an active area of research in the NISQ era [9,10,11].
The BExD approach has been used to identify parameters in electro-chemical models of Lithium-Ion cells [12], pharmacokinetic model [13], Thomson scattering model for nuclear fusion [14] and electrical impedance tomography [15].The probabilistic nature of quantum mechanics naturally lends itself to Bayesian characterization.The BExD technique has been used to estimate the resonant frequency of a single two-level closed system under time-independent Hamiltonian by minimizing the Bayes' risk [16].Granade et al. [17] demonstrated the BExD approach to infer the transition frequency and T 2 dephasing term on a canonical, two-level system, and built the open-source framework QInfer [18].McMichael et al. [19] utilized BExD to determine optimal controls in Ramsey interferometry to estimate precession frequency and T ‹ 2 decoherence time, and demonstrated a five-factor speed up over random sampling methods.Wang et al. [20] demonstrated the BExD framework for quantum sensing of the environment around a Nitrogen-Vacancy (NV) center in diamond and reported a 90% speed up over their heuristic method in estimating the nearby nuclear spins and oscillating magnetic field.Hincks et al. [21] estimated five-parameters in their Hamiltonian and Lindblad model for NV in diamonds with Ramsey and Rabi style experiments where they showed the BExD approach results in a posterior variance that is two orders of magnitude smaller than those obtained using the heuristic methods.Wang et al. [22] performed characterization of an NV center in diamonds by estimating the Rabi frequency and updating their Hamiltonian model to include missing interactions (chirping).Both Wang et al. [22] and Ferrie et al. [16] demonstrated exponential convergence of the BExD (when minimizing posterior variance) up to a saturation point, beyond which the ansatz is no longer able to model the missing interactions.Gester et al. [23] recently employed the BExD approach to calibrate two-qubit entangling operators for a closed trapped-ion systems using 1200 ˘500 experiments in less than one minute.Although the BExD approach significantly reduced the number of experiments required for accurate estimation of the model parameters, the significant speed up is mainly due to the use of analytical solution to their Hamiltonian model.Stace et al. [24] demonstrated the BExD technique on single qubit and synthetic two qubit trapped-ion systems, where they estimated the detuning frequency, Rabi frequency and the coupling strength using optimized bang-bang pulses.They reported difficulty in estimating model parameters where the algorithm fails to converge which is attributed to the identifiability of the parameter due to the states prepared under the pulse parameterization.Although paramount, most studies for characterization of the quantum system often do not consider the identifiability of the parameters under a given experimental setup.Wang et al. [25] demonstrated that the identifiability of closed quantum systems with a single qubit in density matrix formalism (Liouville-von Neumann equation).Zhang et al. [26] investigated the identifiability of both closed and open, dephasing qubit systems and demonstrated the identifiability can be greatly increased by performing appropriate experiments that can be parameterized by control pulses driving the quantum system.
However, the identifiability of multi-level closed and open quantum systems has yet to be explored.Hence, this work demonstrates the identifiability of parameters in the Schroedinger and Lindblad models of multi-level closed and open (with energy decay and dephasing) systems, respectively, and the sufficient conditions on quantum state for the parameters to be identifiable under a set of measurements.Furthermore, it demonstrates the invariance of the conditions to unitary transformations, hence their applicability to a different set of measurements.The BExD approach is then used to identify optimal experiments.Identifying these optimal experimental parameters amounts to maximizing functions of the experiment controls, parameters, and prior that quantifies how informative the experiment would be if performed.These utility functions will be discussed in the following sections.Although demonstrated on a single, qubit problems, its performance for characterizing quantum systems with higher energy levels has yet to be thoroughly investigated.This work aims to address this scarcity and apply the BExD technique for characterizing multi-level quantum systems.We present an approach for iterative pulse extension to accurately estimate parameters that govern processes operating on different time scales.

Quantum Mechanical Model
Consider an open quantum system whose dynamics are modeled by the Lindblad master equation [27] Bρ Bt " ´i rH, ρs `Lpρq where ρ P C 3ˆ3 is the density matrix, r¨, ¨s is the commutator, Hptq P C 3ˆ3 is the Hamiltonian governing the unitary dynamics, Lpρq is the Lindbladian governing the dissipative dynamics.The Hamiltonian in the lab frame consists of the time-independent system Hamiltonian H s and a time-dependent control Hamiltonian H c ptq Hptq :" ωa : a ´χ 2 a : a : aa looooooooomooooooooon Hs `f ptqpa `a: q loooooomoooooon where a and a : are the lowering and raising operators, ω is the 0-1 transition frequency, χ is the anharmonicity, and f ptq " 2RepΩptqe iw d t q is the control function with drive frequency ω d and amplitude Ωptq :" pptq `iqptq where pptq, qptq P R.Under the rotating wave approximation (RWA), we obtain the Hamiltonian in the rotating frame Hptq " pω ´ωd q a : a ´χ 2 a : a : aa loooooooooooooomoooooooooooooon Hs `Ωptqa `Ωptqa : loooooooomoooooooon Hcptq " ∆a : a ´χ 2 a : a : aa `pptq `a `a: ˘`iqptq `a ´a: where the rotation frequency in the RWA is taken to be the drive frequency ω d , and ∆ :" ω ´ωd is the detuning between the qubit 0-1 transition and the drive frequencies.This model was experimentally validated for a transmon qutrit at Quantum Device and Integration Testbed (QuDIT) at LLNL [28].The decoherence is introduced through the Lindblad jump operators for decay (L 1 :" 1 ?T1 a) and dephasing (L 2 :" 1 ?T2 a : a) where T 1 and T 2 are the decay and dephasing times, respectively.The Lindbladian is then written as Lpρq :" where t¨, ¨u is the anti-commutator.The control pulse is modeled as a set of piece-wise constant pulses, each with pptq, qptq, and duration ∆t.The pulse resolution can be increased by increasing the number of segments.In general, real experiments require 1 ns resolution to play pulses on hardware.This generalization can be seen as a bang-bang control where each application uses a different amplitude and duration for the control.Equation 1 is vectorized with ρ " vecpρq and evolved until final time T using matrix exponentiation by means of the Padé approximant where N T is the number of piecewise-constant segments and the Liouvillian superoperator Lptq is defined as and I is rank-3 identity matrix.Given a density matrix, the probability of measuring the system in a particular state y is given as Prpy|ρq " Tr pM y ρq (7) where M y is the measurement operator corresponding to the observable y.In this work, we assume projective measurements in the z-basis and take M y " |e y y xe y | , y " 0, 1, 2 where |e i y is the i th canonical unit vector; then each element of Prpy|ρq " diagpρq denotes the probability of measuring that state and will be used to compute the likelihood in the subsequent sections.
3 Identifiability of Quantum Mechanical Models Before characterizing a system described by a mathematical model, it is imperative to determine the global and unique (termed global structural) identifiability of the parameters in such a model.Some techniques for determining the identifiability of a model include those based on Laplace transforms [29], similarity transformation [30], generating series [31], differential algebra [32,33] and the Taylor series approach [34].To determine the identifiability of the models of open and closed quantum systems, we will employ the Taylor series approach [34,35,36].The approach relies on the observations being a unique analytic function of time, hence their derivatives must also be unique.The uniqueness of this representation will guarantee the structural identifiability of the model.Consider an observable/measurable quantity ypt; θq that is dependent on model parameters θ.The Taylor series approach expands this quantity about time t 0 A sufficient condition for global structural identifiability is given by the uniqueness of the expansion coefficients where m is the smallest positive integer such that the resulting system is solvable for the parameters.We emphasize that the order of expansion here is used only to establish identifiability conditions.This minimum order of expansion is not known a priori but as will be shown, first order and second order expansions are sufficient for closed and open quantum systems, respectively.We should mention that the identifiability conditions are derived for models in the lab-frame and, therefore, do not make any approximations such as neglecting the rapidly oscillating terms in the RWA.

Model Identifiability: Schroedinger Equation
Consider a closed quantum system modeled by the Schroedinger equation 9 |ψy " ´iH |ψy, driven by a known pulse f ptq with a pulse amplitude Ωptq and the probability of measuring state |iy by y i pt; θq :" Prp|iy , t; θq " xψpt; θq| M i |ψpt; θqy.The pulse amplitude Ωptq is taken to be constant for all times.An alternative parameterization for Ωptq can also be used, in which case, the resulting Taylor series expansion will yield additional pulse dependent terms that simply scale the coefficients a i,k in Eq. ( 8) and, therefore, do not alter the identifiability conditions.We expand the measurement about t 0 " 0 with F :" f pt 0 q, which yields the first three expansion coefficients We note the zeroth order terms (a i,0 pθq, i " 0, 1, 2) do not yield any information on the parameters of the model.Defining ψpt 0 q " tψ 0 , ψ 1 , ψ 2 u P C 3 , the first order terms can be written as which uniquely define F when Impψ i ψ j q ‰ 0 for i ‰ j.The second order terms are given by Collecting all the known and identified terms in a right-hand-side S P R 2 , we can write the identifiability conditions in matrix form Notice that the rows of A c are linearly dependent.Since the drive F is known (or can be identified as indicated by the first order terms), the sufficient conditions to identify ω and χ are Repψ 0 ψ 1 q ‰ 0 and Repψ 1 ψ 2 q ‰ 0, respectively.Furthermore, the sufficient conditions for identifiability of χ requires that the conditions for ω also be satisfied.Hence, the parameters are globally and uniquely identifiable and the model is globally structrually identifiable, when the control pulse is able to prepare a state that satisfies these conditions.

Model Identifiability: Lindblad Master Equation
We now determine the identifiability of the Lindblad master equation.Consider an open quantum system defined by Eq.
for the remaining unknown parameters where the remaining terms involving the known or identifiable parameters (F and τ 1 ) are grouped in S 2 P R 3 .These second order terms show that ω, χ, τ 2 are globally identifiable if A 2 is non-singular.Hence, the control function should be able to drive the system to such a state.Similarly, the identifiability of the density matrix formalism for closed system (i.e., Liouville-von Neumann equation) can be determined by fixing τ 1 " τ 2 " 0. In this case, the first order conditions only show identifiability of F, like in the case of state vector formalism (i.e., Schroedinger equation §3.1), and a second order expansion is needed to determine identifiability of the remaining parameters ω, χ.Note that these identifiability conditions in both density matrix and state vector formalism are equivalent.Letting ρpt 0 q " |ψpt 0 qy xψpt 0 q|, we obtain the condition for Impρ ij q " Impψ i ψ j q ‰ 0 for i ‰ j for parameter F , Repρ 01 q " Repψ 0 ψ 1 q ‰ 0 for parameter ω and Repρ 21 q " Repψ 2 ψ 1 q ‰ 0 for parameter χ.We should mention that these conditions also hold for other measurement basis Ă M by applying an appropriate unitary transformation U such that Ă M i " U M i U : since y i ptq " Tr `U M i U : ρptq ˘" Tr `Mi U : ρptqU ˘" Tr pM i r ρptqq where the conditions on the quantum state now apply in the transformed basis.These conditions show that the model parameters are identifiable by measuring the population in a single basis without needing to estimate the full density matrix via state or process tomography, which can be expensive.

Bayesian Experimental Design
Having demonstrated the identifiability of the model parameters is dependent on the quantum state, and therefore, the control applied (or experiment performed), we present an approach for automatic, online design of experiments that allow the parameters to be identified.Consider a set of model parameters θ to be identified, set of measurements y from an experiment parameterized by a set of controls (i.e.controllable parameters) ξ.We begin with the Bayes' theorem where Prpθq, Prpθ|y; ξq, and Prpy|θ; ξq are the prior, posterior distribution and likelihood, respectively.The normalizing constant is the marginalized probability distribution over parameter space (also referred to as the evidence) and is defined as The Bayesian experimental design (BExD) framework aims to estimate the probability distributions of the model parameters (i.e. the posterior distribution) by identifying the optimal controls of an experiment that provides maximum information on the parameters.The measure of effectiveness of the experiment (parameterized by ξ) if the outcome y is observed is called the utility function U py|ξq.Hence, an optimal experiment is one that maximizes the expectation of this utility function over all possible observations U pξq " E Prpy|ξq rU py|ξqs " Typical choices for the utility function include: 1) negative variance of the posterior, 2) Kullback-Leibler (KL) divergence and 3) Wasserstein distance; other utility functions are also possible [37,38].Minimizing the variance seeks to minimize the spread about the mean, thereby reducing the uncertainty in the model parameters.Meanwhile, maximizing the KL-divergence and Wasserstein distance functions seeks controls that maximize the distance between the prior and the posterior distributions, thereby yielding the largest update over the prior distribution.
For all such utility functions, the integrals in Eq. ( 22) are often analytically intractable and the curse of dimensionality renders standard numerical integration inefficient.Hence, the integral is often approximated using Monte Carlo sampling by drawing N samples θ i " Prpθq where w i are the weights and ř i w i " 1.These MC samples are also used to construct an approximate discrete representation of the PDF where δp¨q is the Dirac delta function.The search for such optimal control ξ ‹ requires the solution to an optimization problem Eq. ( 25) that can be high-dimensional and typically non-convex and multi-modal [39].
For these reasons, we solve the optimization problem using a non-gradient based differential evolution algorithm.The applicability of the evolutionary method for optimization of utility function in BExD was demonstrated by Hamada et al. [40].Schoneberger et al. [39] demonstrated the robustness of non-gradient based optimization algorithm for optimal experimental design.An experiment is then performed using the optimal control ξ ‹ , and the measurements y are used to construct the likelihood and update the posterior distribution at the n th epoch/iteration as To reduce computational cost associated with evaluating the likelihood in regions of low posterior density, the Monte Carlo samples are resampled to redistribute the particles to regions of high posterior density.We make use of the resampling algorithm proposed by Liu and West [41].The expectation of the inferred posterior p θ :" E Prpθ|yq rθs yields the unbiased estimator for the model parameters.Figure 1   The Kullback-Leibler (KL) divergence is a measure of distance between two probability distributions and is defined as Note that the KL divergence is not symmetric in the two distributions.At the n th epoch, Prpθ|y; ξq " Prpθ|y n ; ξq is the posterior (i.e.Eq. ( 26) for some control vector ξ) and Prpθq " Prpθ|y n´1 ; ξ ‹ n´1 q is now the prior (i.e.posterior from previous epoch).Taking the expectation of the KL divergence over all possible outcomes y yields an expression for the mutual information Hence, maximizing the mutual information amounts to identifying control parameters that provide the maximum information gain over the prior distribution.

Likelihood Function
For clarity of presentation, we denote by r θ the true distribution of the parameters (i.e. the parameters of the quantum testbed).For a total of N independent shots, let y i " Prpy| r θ; ξq be a possible outcome of a shot (i.e., the state measured for a shot) after applying control ξ, n i the number of times y i is observed with ř i n i " N .Then, the likelihood can be modeled by a multinomial distribution [17,38] 7) and the numerical model for a given θ and ξ.It should be mentioned that the multinomial coefficient can be omitted since the weights in Eq. ( 24) are renormalized at each epoch.For a large number of draws N , this omission may lead to numerical underflow as the resulting expression would produce near zero probabilities, particularly when Prpy i |θ; ξq is small.Alternatively, under certain assumptions, the multinomial distribution asymptotically tends to a Gaussian distribution [42] Prpy|θ; ξq " N ´Prpy| r θ; ξq, Σ ¯(30 where Σ is the covariance that accounts for errors between the numerical model and the 'true' model (i.e.experiment).This model of the likelihood has been used extensively in Bayesian experimental design [37,43,38] and in Markov Chain Monte Carlo (MCMC) methods.The exact form of Σ is not known a priori and must be estimated.Let Ypθ; ξq " tPrpy i |θ; ξq, i " 1, . . ., N u and r Ypξq " ) be vectors of probabilities computed using the numerical model and experimentally measured, respectively, after applying control ξ.Then, we estimate Σ by means of Monte Carlo sampling over the prior θ with N θ samples and random control ξ with N ξ samples.The sample covariance matrix is then computed as ´Ypθ i ; ξ j q ´r Ypξ j q ¯´Ypθ i ; ξ j q ´r Ypξ j q ¯T (31)

Case Studies for Quantum Characterization
We now present the problem formulation for verification of our characterization technique using synthetic test problems.
We consider an open system modeled by the Lindblad master equation with true value of the parameter set θ " tω, χ, T 1 , T 2 u to be that of the standard QPU in Quantum Design and Integration Testbed (QuDIT) at LLNL [44].We mainly consider uncertainties in these four parameters, which are the most pronounced noise sources in the experiments, especially for a quantum gate in multi-level system.The QPU has a single transmon in 3D aluminum cavity whose parameters are given in Tab. 1.We consider the cases where the exact model parameters are assumed to be point- estimates or following some distribution.When treating the parameters as point estimates, we take r θ " ) .Each experiment is then performed using this fixed parameter set.In the case where the parameters are assumed to be random variables, we consider their exact distribution to be r θ " N pµ θ , Σ θ q :" N pr ω, 10 ´6q ˆN pr χ, where each experiment is performed using a random sample from the distribution which introduces stochasticity into the model problem.The prior distribution (i.e.distribution at the initialization) is taken to be Prpθq " U ω p3.5, 4.5q ˆUχ p0.1, 0.2q ˆUT1 p30, 60q ˆUT2 p20, 40q Since the 'true' distribution is not known in practice, the use of an uninformative prior, like the uniform distribution, avoids bias in the inference.In addition to the uncertainty introduced by imposing a distribution over the model parameters, we also introduce shot noise (i.e.sampling uncertainty) by sampling a finite number of outcomes N (in Eq. ( 29)).When omitting shot noise and using exact Prpy| r θ; ξq (i.e. when N Ñ 8), we employ the Gaussian model for the likelihood (Eq.( 30)).We refer to those experiments performed using point-estimates of parameters and exact measurements (i.e.N " 8) as exact experiments and those experiments with uncertainty (either parameter or sampling) as stochastic experiments.We use a total of 2000 Monte Carlo particles and limit the total number of epochs to 500.The results with different number of particles is given in Appendix A. It was previously mentioned that the topology of the utility function can be highly multimodal [39], hence we employ the non-gradient based differential evolution algorithm, since it is easily parallelizable and less susceptable to convergence to local extrema.Each experiment will be defined by a different control pulse that is parameterized as a sequence of piecewise constant segments in the rotating frame, with a drive frequency ω d and each segment defined by pptq, qptq and ∆t as shown in Fig. 2.Although other pulse parameterization, such as spline parameterization [45], can also be used, the piecewise-constant pulses simplies matrix exponentiation since it can be performed for each segment simultaneously and in parallel.For a finely resolved pulse, this piecewise-constant parameterization will lead to an overwhelming number of parameters that must be optimized, and therefore, a parameterization is more applicable.Hence, for an N p segment pulse, the total number of parameters defining the experiment are N T " 3N p `1.The bounds on the pulse parameters are pptq, qptq P r´12, 12s MHz, ∆t P r1, 30s µs and ω d P r3.5, 4.5s.We assume each experiment begins with the qutrit in the ground state |0y.For verification purposes, we compare the converged distributions to the exact, imposed distribution and analyze the error in the model predictions for a reference quantum gate using Swap-02 gate.
We define the error in the mean as δpp µq " E rr µs ´p µ where p µ is the unbiased esimator for µ P θ.Furthermore, we quantify the errors in the model prediction at each time t as Ep Ă |iy, |iyq " E r θ " Prp|iy | r θ; ξ 02 q ı ´Eθ rPrp|iy |θ; ξ 02 qs where ξ 02 is the control pulse for a Swap-02 gate computed using quantum optimal control framework Quandary [45] with the 'true' parameters taken to be point estimates r θ.The control pulse and the evolution of the population for the Swap-02 gate are shown in Fig. 3.The short duration of the control pulse (100ns) does not clearly display the dissipative dynamics, such as qubit energy decay over time.For this reason, we repeat the gate 100 times, for a total duration of 10µs to allow the dissipative dynamics to be more prevalant and amplify the errors.

Results: Characterizing Qutrit Systems
We consider Bayesian experimental design under a fixed control pulse and an approach for iterative control pulse extension.In the fixed case, the number of control pulse segments remains constant, whereas in the iterative pulse extension case, the number of pulse segments increase with epochs.

Fixed Pulse Parameterization
We first investigate the effectiveness of our approach under a fixed pulse parameterization and the problem formulation presented in Sec. 5. Figure 4 shows the mean error and the standard deviation of the estimated distributions using different number of control pulse segments under different sources of uncertainty.Immediately we see that the mean error and the standard deviation for ω and χ distributions are orders of magnitude smaller than those for the distributions of T 1 and T 2 .Furthermore, variances of all four parameter distributions decrease near-monotonically.In both the exact and stochastic experiments, the transition frequency and anharmonicity show a faster rate of convergence than the decoherence parameters, with several orders of magnitude reduction in only the first few epochs.We see a better overall convergence in T 1 and T 2 when using a more robust pulse parameterization with a larger number of segments, but with a slight decrease in the accuracy of ω and χ estimates.Control pulses with a greater number of segments lead to longer pulse durations and better expose the dissipative dynamics that are present on larger timescales.This suggests that an adaptive pulse parameterization -where short pulse sequences are used to accurately estimate ω and χ, then longer sequences are used to improve the T 1 and T 2 estimates -may be beneficial.Figure 5 shows the exact (imposed) and the converged distribution for the four parameters using a different number of pulse segments N p and shots N s .We see a good agreement in the distributions for ω and χ whereas the distributions of T 1 and T 2 feature larger variances.The mean of T 1 has better agreement with that of the exact distribution, whereas the distribution of T 2 displays both a larger error in the mean and a wider spread.As expected, the uncertainty introduced by sampling a finite number of shots introduces additional errors; this often manifests as a larger shift in mean or larger variances.Furthermore, we see that shorter pulse sequences yield better estimates for ω and χ whereas longer pulse sequences result in better estimates for T 1 and T 2 .The statistical moments of the converged distribution reported in Tab. 2 show that the frequency parameters are less susceptable to shot noise than the decoherence times.When model parameters are not unique, quantifying the accuracy using only errors in the parameters is not sufficient.
In such cases, it is imperative to quantify the accuracy of model response/predictions. Figure 6 shows the probability density of the errors in the population due to different sources of uncertainty, pulse segments, and number of shots.Furthermore, the mean and standard deviation of the errors is also presented in Tab. 3. In the case of exact experiments, we see that longer pulse sequences yield model parameters that result in more accurate predictions (i.e. higher densities for smaller errors).This is due to the more accurate estimation of T 1 and T 2 and therefore, more accurate representation of the dissipative dynamics that are present in longer duration pulses (i.e., larger gate repitition).When performing estimation using stochastic experiments, we see that the error distribution for all population and pulse parameterization are similar, and are an order of magnitude larger than that of exact predictions.Furthermore, despite the lack of accuracy in the decoherence parameters compared to the frequency parameters, the mathematical model with estimated parameters yields sufficiently accurate predictions with mean errors on the order of 10 ´4 and variances on the order of 10 ´3.  Figure 8 shows the distribution of the parameters and prediction errors whose mean and standard deviation are presents in Tab. 5. Firstly, we see that in the absence of shot noise, the parameter mean is better estimated for all four parameters, whereas the presence of shot noise results in a shift in mean value and a larger spread.However, the distribution (and the statistical moment) of the error shows that iterative pulse extensions lead to an order of magnitude reduction in error mean and variance (in the absence of shot noise) when compared to the fixed pulse parameters; the mean errors are also an order of magnitude smaller with comparable variances on the order of 10 ´3.This suggests that the BExD approach is able to reliably estimate the model parameters, even though the shot noise significantly affects the accuracy of the estimated parameters Furthermore, despite severe undersampling (using only 500 shots), the errors in the predictions from the calibrated model are well within typical bounds « 99.5% [46] for several applications of NISQ testbeds.

Conclusion
The conventional techniques for quantum characterization, whether probabilistic or deterministic, rely on measurements from a set of experiments that are designed a-priori.These measurements are then used to perform an offline calibration of the mathematical model to reconstruct the measurements from the experimental testbed.In most cases, the theoretical identifiability of the model parameters is neglected.Even when the identifiability conditions are known, manual design of experiments satisfying such conditions is not trivial.Hence, this work presented an automatic, online Bayesian framework for characterizing open qutrit systems described by the Lindblad master equation.The approach automatically identifies optimal robust experiments that provide maximum information on the model parameters to be estimated.Furthermore, we proved the theoretical identifiability of the Schroedinger and Lindblad equations given a measured observable, and derived the sufficient conditions for the quantum state to establish the global, unique identifiability of the model parameters.
We demonstrated our approach to identify the frequencies and decoherence parameters on a numerical test problem subjected to both parameter uncertainty and shot noise.The characterization study demonstrated that different parameters exhibit different rates of convergence and are dependent on both the pulse parameterization and magnitude of shot noise.It was seen that accurate estimation of the frequencies required shortened pulse sequences whereas estimation of decoherence parameters required longer pulse sequences.We introduced an iterative pulse extension approach to improve the convergence of all estimates.The results show that the Bayesian experiment design framework with iterative pulse extensions is able to accurately estimate all four parameters in only 500 epochs (i.e.500 different experiments).The approach is better able to estimate the distribution of the system frequencies than the decoherence parameters.Desipite this, the BExD approach yielded robust estimates of the parameters such that the mean prediction error and variances were on the order of 10 ´4.The results showed that the online BExD approach is able to design robust, optimal experiments to accurately characterize the quantum system described by the Lindblad master equation.
Ongoing and future work includes experimentally demonstrating the technique on multi-qubit systems.

A Robustness of Bayesian Experimental Design To Parameter Settings
Figure 9 shows the distribution of parameters and prediction errors from an ensemble (with different random initialization of 2000 particles) of six runs using the same control parameterization as those in Section 6.2.Here, the experiments are considered to be stochastic with parameter uncertainty but without shot noise.The distributions are similar to those in Fig. 8 and shows the estimation is robust to different initialization strategies.Figure 10 shows the distribution of parameters and prediction errors when using different number of Monte Carlo particles.Here, the experiments are again considered to be stochastic with parameter uncertainty but without shot noise.The distributions show that the frequency parameters are estimated more accurately and the decoherence times, with only a few number of samples.The under-resolved parameterization of the distribution (Eq.( 24)) and estimation of expectations by Monte Carlo sampling (Eq.( 23)), leads to this decrease in estimation accuracy, which in turn leads to an increase in prediction errors.However, this shows that the BExD approach does not require an overwhelming number of particle but a sufficient number to cover the parameter space.

Figure 1 :
Figure 1: The workflow for our Bayesian experimental design framework

Figure 2 :
Figure 2: The control pulse parameterization used in this work with Ωj " pjptq `iqjptq

Figure 3 :
Figure 3: Time evolution of: (a) controls in rotating frame and (b) population on each state while one Swap02 gate is being played.

Figure 4 :
Figure 4: Convergence of mean error δpp µq and standard deviation p σ of the model parameters when performing exact (top) experiments and experiments with parameter uncertainty (bottom) using pulses with: (a,c) three segments and (b,d) six segments.The solid and dashed lines represent the convergence in mean error and standard deviation, respectively.

= 6 Figure 5 :
Figure 5: Converged parameter distribution when performing stochastic experiments using a different number of pulse segments (Np) and different number of shots Ns.Note the distributions are translated by the mean of the exact distributions.

Figure 6 :
Figure 6: Distribution of the error in the population when using a different number of pulse segments and performing: a) exact experiments, b) experiments with parameter uncertainty, and c) experiments with parameter uncertainty and shot noise.

Figure 8 :
Figure 8: (Top) Converged parameter distribution when performing stochastic experiments and (Bottom) distribution of the error in the population when using different number of shots.Note the parameter distributions are translated by the mean of the exact distributions.

Figure 9 :
Figure 9: (Top) Converged parameter distribution over an emsemble of six runs with when performing stochastic experiments with parameter uncertainty and (Bottom) distribution of the error in the population.Note the parameter distributions are translated by the mean of the exact distributions.
presents the complete workflow for Bayesian experimental design.
defined by Prpy|θ; ξ, N q "ˆN n 1 , n 2 , . . ., n M ˙N ź from experiments, of observing y i and Prpy i |θ; ξq is the probability of observing y i computed using Eq. (

Table 1 :
Parameters of the QPU at LLNL's Quantum Design and Integration Testbed.

Table 2 :
Mean and standard deviation (in the brackets) of the parameters estimated using a fixed pulse parameterization

Table 4 :
Mean and standard deviation (in the brackets) of the parameters using iterative pulse extensions.

Table 5 :
Mean and standard deviation (in the brackets) of the error using iterative pulse extensions.