Resonant Coupling Parameter Estimation with Superconducting Qubits

Today's quantum computers are comprised of tens of qubits interacting with each other and the environment in increasingly complex networks. In order to achieve the best possible performance when operating such systems, it is necessary to have accurate knowledge of all parameters in the quantum computer Hamiltonian. In this article, we demonstrate theoretically and experimentally a method to efficiently learn the parameters of resonant interactions for quantum computers consisting of frequency-tunable superconducting qubits. Such interactions include, for example, those to other qubits, resonators, two-level state defects, or other unwanted modes. Our method is based on a significantly improved swap spectroscopy calibration and consists of an offline data collection algorithm, followed by an online Bayesian learning algorithm. The purpose of the offline algorithm is to detect and roughly estimate resonant interactions from a state of zero knowledge. It produces a square-root reduction in the number of measurements. The online algorithm subsequently refines the estimate of the parameters to comparable accuracy as traditional swap spectroscopy calibration, but in constant time. We perform an experiment implementing our technique with a superconducting qubit. By combining both algorithms, we observe a reduction of the calibration time by one order of magnitude. We believe the method investigated will improve present medium-scale superconducting quantum computers and will also scale up to larger systems. Finally, the two algorithms presented here can be readily adopted by communities working on different physical implementations of quantum computing architectures.

Today's quantum computers are comprised of tens of qubits interacting with each other and the environment in increasingly complex networks. In order to achieve the best possible performance when operating such systems, it is necessary to have accurate knowledge of all parameters in the quantum computer Hamiltonian. In this article, we demonstrate theoretically and experimentally a method to efficiently learn the parameters of resonant interactions for quantum computers consisting of frequency-tunable superconducting qubits. Such interactions include, for example, those to other qubits, resonators, two-level state defects, or other unwanted modes. Our method is based on a significantly improved swap spectroscopy calibration and consists of an offline data collection algorithm, followed by an online Bayesian learning algorithm. The purpose of the offline algorithm is to detect and roughly estimate resonant interactions from a state of zero knowledge. It produces a square-root reduction in the number of measurements. The online algorithm subsequently refines the estimate of the parameters to comparable accuracy as traditional swap spectroscopy calibration, but in constant time. We perform an experiment implementing our technique with a superconducting qubit. By combining both algorithms, we observe a reduction of the calibration time by one order of magnitude. We believe the method investigated will improve present medium-scale superconducting quantum computers and will also scale up to larger systems. Finally, the two algorithms presented here can be readily adopted by communities working on different physical implementations of quantum computing architectures.
Small-and medium-scale quantum computers built out of superconducting circuits [14][15][16][17] and comprising up to a few tens of physical qubits have been in operation for the past few years [18][19][20][21][22][23]. The two main reasons for the success of superconducting technologies in the quantum computing arena are the pre-existing facilities for scaling up fabrication, due to methods being similar to silicon (Si) technology, as well as the favorable coherenceto-gate-time ratio. For state-of-the-art superconducting qubits with very short gate times of less than 100 ns (e.g., see Ref. 24), such a ratio reaches values close to 1000. This means that thousands of one-and two-qubit gates can be performed within the qubit lifetime, with fidelities in excess of 99.4 % [22]. Very recently, a major milestone in quantum computing has been reached -quantum supremacy -where a 53-qubit superconducting * Corresponding author: matteo.mariantoni@uwaterloo.ca quantum computer has been utilized to sample the output of a pseudo-random quantum circuit [25]. As a result of these advances, quantum computing has transitioned into the so-called noisy intermediate scale quantum (NISQ) era [26].
Despite these successes, the high-fidelity operation of medium-and large-scale quantum computers is accompanied by the daunting task of calibrating numerous physical qubits. In particular, calibrating tunable qubits requires the estimation of resonant interaction parameters such as the resonance frequency and coupling coefficient between pairs of interacting qubits [27,28], a qubit and a resonator [18], two-level state (TLS) defects [29][30][31], or substrate and box modes [32]. These calibrations are necessary to implement two-qubit gates and avoid loss of quantum information due to spurious interactions leading to coherent or incoherent errors. TLS defects, especially, are a pervasive source of errors in superconducting architectures which must be remediated [33].
In this article, we study theoretically and demonstrate experimentally a data-efficient and automated method for identifying and estimating the parameters of resonant interactions based on swap spectroscopy [18,34]. We realize swap spectroscopy by performing energy relaxation time T 1 measurements of a frequency-tunable Xmon transmon qubit [24] at different qubit frequencies.
The identification and estimation method is divided into two parts: an offline data collection algorithm [35] and an online Bayesian learning algorithm [36,37]. Both algorithms are based on the dynamics of interacting quantum systems. The former is used from a state of zero knowl-edge to roughly identify resonance parameters; the latter focuses on improving the estimate of those parameters. In this context, the term "online" means that measurements taken during the execution of the algorithm inform the subsequent ones. For the "offline" method, the execution of the entire algorithm is predetermined.
By means of our parameter-estimation method, we can to shorten the calibration time of an Xmon transmon qubit significantly. The offline data collection algorithm makes it possible to reduce the number of measurements by a square-root factor when compared to a traditional swap spectroscopy calibration. In our experiment, this algorithm takes ≈ 30 min to detect resonances in a 1 GHz bandwidth; one order of magnitude less time that it would take previously. The online Bayesian learning algorithm runs in ≈ 25 s per resonance and improves the estimation accuracy to be equivalent to that of a high-resolution traditional swap spectroscopy.
Our method is not confined to the realm of superconducting quantum computing. In fact, it can easily be adopted by practitioners working on different physical implementations of quantum computing architectures such as trapped ions and semiconductor qubits [1].
This article is organized as follows: in Sec. II, we explain qubit calibration in frequency-tunable architectures; in Sec. III we summarize the working principle of traditional swap spectroscopy, explaining why it is inefficient; in Sec. IV, we introduce the offline octave sampling algorithm (Subsec. IV A) and show our experimental implementation of the method to find the interaction parameters between an Xmon transmon qubit and three resonance modes (RMs) (Subsec. IV B); in Sec. V, we explain the online Bayesian learning algorithm (Subsec. V A) and demonstrate its performance for refining the estimate of the resonance parameters (Subsec. V B); in Sec. VI we discuss additional concerns with the algorithms and the relevance of our methods for quantum computing; finally, in Sec. VII, we provide an outlook and conclusions.

II. QUBIT CALIBRATION IN FREQUENCY-TUNABLE ARCHITECTURES
A fundamental requirement to the operation of a quantum computer is the proper calibration of the physical qubits in the system. This calibration includes many specific operations. One of the most basic tasks, for example, is to run a Rabi experiment on each qubit. This allows the determination of some experimental parameters needed to set up, e.g., a π-pulse and perform a measurement. Once this first task has been realized, further experiments can refine the knowledge of the pulse amplitude, rotation axis, measurement parameters, etc. Finally, a full calibration requires knowing the precise parameters of the total system Hamiltonian, allowing for the systematic optimization of the fidelity of one-and two-qubit gates as well as measurement.
In a frequency-tunable superconducting qubit architecture such as the Google architecture [25] or the one used in this work, an additional degree of freedom must be considered during calibration: the qubit frequency. Xmon transmon qubits are one example of tunable qubits [24]. In this design, an on-chip capacitive island made from aluminum (Al) is coupled in parallel to a DC superconducting quantum interference device (SQUID) comprised of two Josephson tunnel junctions in parallel, forming a superconducting loop [17]. Simply put, an Xmon transmon qubit is a quantum anharmonic oscillator, characterized by a non-equally spaced ladder of quantum states. The frequency (i.e., energy) difference f q between the energy ground state |g and first excited state |e differs from that between |e and the second excited state |f by the so-called qubit anharmonicity α [38]. The qubit transition frequency f q is controlled in situ by applying a local external magnetic flux that threads the DC SQUID, tuning the Josephson energy E J and therefore the level separation. This degree of freedom leads to a few distinct advantages to the operation of a quantum computer.
For instance, frequency tunability allows for tunable qubit-qubit interactions because the effective coupling strength between two qubits depends on the frequency difference between them. This enables the implementation of several types of two-qubit gates such as the controlled-phase (CPHASE) gate, which takes advantage of state |f as an auxiliary state [18,27,28,39], as well as the √ iSWAP and iSWAP gates [18]. In addition, setting the frequency of spatially neighboring qubits away from each other helps avoid control crosstalk and frequency crowding issues, with the latter being endemic in fixedfrequency systems [40].
Another advantage inherent to frequency tunable architectures has to do with energy relaxation. On-chip superconducting qubits interact with a distribution of TLS defects, which are present in the various amorphous dielectric materials surrounding the qubit metallic structures (e.g., Si and Al oxides). While the microscopic origin of TLS defects is still under debate, their effect on the qubit can be broadly categorized either as causing frequency noise or excitation noise [29]. In particular, TLS defects interacting semi-resonantly with a qubit can coherently or incoherently exchange energy with it resulting in a lossy channel. The ability to set the frequency of a qubit away from that of TLS defects is therefore desirable, and realizable only with tunability.
Calibrating qubits to implement two-qubit gates or to avoid TLS defects is a parameter estimation problem. We need to determine the Hamiltonian parameters that define the resonant interactions between a qubit and another system. In all the aforementioned cases, two parameters must be found: the resonance frequency and coupling strength of the interaction.
Historically, this kind of calibration has been realized using swap spectroscopy. Unfortunately, traditional swap spectroscopy is inefficient in the amount of data it requires, and therefore slow. This is inconvenient for mul-tiple reasons. First, as the number of qubits in a system grows, so does the number of calibrations that must be performed. This is particularly relevant to qubit-qubit coupling calibration, which cannot be performed in parallel on all qubits. Second, TLS defects in the environment are known to fluctuate over time [31,33,41,42]. Similarly, f q itself can shift in time. The identification of resonant interactions must therefore be repeated at regular intervals. We thus require a robust, accurate, and time-efficient method to identify the parameters associated with resonant interactions.
In order to test our parameter-estimation method and compare it against traditional swap spectroscopy, we use an Xmon transmon qubit and make it interact with a synthesized resonance mode. The mode, which emulates the interaction with another qubit, resonator, or TLS defect, is created by applying a coherent drive with a microwave source to the qubit under test. The synthesized resonance mode is a convenient and flexible tool to test our method since we can arbitrarily change its resonance frequency by tuning the source frequency as well as its coupling strength by changing the emitted source power.

III. TRADITIONAL SWAP SPECTROSCOPY
Swap spectroscopy is an experimental method that allows exploring the environment of a qubit at various frequencies by using the qubit itself as a probe. Traditionally, swap spectroscopy has been used to select the operating frequency of qubits, making it possible to avoid TLS defects or regions of low T 1 . Additionally, it has been used to explore resonant interactions, such as those with other qubits [28] or resonators [34]. Performing swap spectroscopy requires a minimally calibrated qubit and, thus, is suitable as a tune-up experiment.
In a swap spectroscopy experiment the qubit is initialized at the so-called idle frequency. A π-pulse is then applied to the qubit, energizing it from |g to |e . At the end of the π-pulse, a flux pulse is applied to the DC SQUID in order to tune the qubit to a different frequency, the probe frequency f p (A). This procedure requires knowledge of the correspondence between the qubit frequency and pulse amplitude A, which can be calibrated via reg- The initial π-pulse (red) excites the qubit, which is initialized at the so-called idle frequency. The flux pulse (blue) changes the qubit transition frequency fq to the probe frequency fp(A) for a duration t. The qubit is then measured (green) after being set back to its idle frequency. ular pulse spectroscopy (see App. B). After a time t, the flux pulse is turned off and the qubit is measured back at the idle frequency. This pulse sequence is illustrated in Fig. 1. Note that using a flux pulse to set f p (A) presents advantages over quasi-statically changing the idle qubit frequency by means of a DC current to the SQUID. Namely, it avoids recalibrating the π pulse and measurement pulse at each qubit frequency.
In a traditional experiment, t and f p (A) are swept linearly over a desired range and the qubit is measured at each point, recording how many shots correspond to an excited or ground state, n e or n g , respectively. As a result of measuring the qubit in the Z basis, swap spectroscopy is mostly insensitive to qubit dephasing. Figure 2 (a) shows the result of a typical swap spectroscopy experiment, with data taken between 4.146 GHz and 5.170 GHz and for times up to 500 ns. Resonant couplings appear as oscillations, or chevron patterns, of the measured average population P e = n e /(n e + n g ) in time. For example, on the far right of the spectrum it is possible to observe very fast oscillations, corresponding to a strong coupling of g ≈ 40 MHz between the qubit and the measurement resonator (see Appendix A for details on the sample layout and experimental setup). To the left of the resonator we observe a slower oscillation corresponding to a weaker interaction between the qubit and synthesized resonance mode. Finally, at an even lower frequency, around 4.35 GHz, we observe a "streaky" structure. In this region, the qubit excitation is lost faster than elsewhere, and we cannot observe any oscillation.
The features observed in Fig. 2 (a) demonstrate a selection of possible resonant interactions: strong interactions, where g 1/T 1 resulting in multiple coherent oscillation cycles and weak interactions appearing as regions of lower T 1 . Neither is ideal for the operation of a qubit. In the case shown in Fig. 2, the best choice for the qubit idle frequency is around 4.6 GHz; far away from any unwanted couplings.
The data in Fig. 2 (a) gives us a rough idea about the parameters of any possible resonance modes coupled to the qubit within the measured spectrum. It is hard to tell, however, that there are in fact two resonance modes at 4.8 GHz, or what the frequency of the oscillation for the resonator is. A more detailed scan, such as the one in Fig. 2 (b), might be necessary to estimate the parameters with sufficient accuracy. Traditional swap spectroscopy, with data taken in a linear grid, is an effective method to detect and estimate resonance modes. We will show in Sec. IV, however, that the traditional method is inefficient, and that there exists a much better way to perform this task.

IV. OFFLINE OCTAVE SAMPLING
The offline octave sampling algorithm has a similar objective as swap spectroscopy, i.e., to determine if there are any systems interacting resonantly with the qubit and provide an estimate for their coupling parameters. However, we want to achieve this purpose in a more efficient fashion by acquiring fewer data and therefore saving time. Note that the basic experiment employed, that is the pulse sequence illustrated in Fig. 1, is the same as for swap spectroscopy. The difference lies in how the data is sampled. Whereas traditional swap spectroscopy samples the frequency-time space in a regular grid, octave sampling takes advantage of TLS dynamics to acquire as few data as possible.

A. Theoretical Method
In order to explain the data collection strategy, we analyze the time dynamics of the systems at play. Since we are searching for resonant interactions with a qubit, we work in a single-excitation manifold (|g ↔ |e ). Thus, even if a system is characterized by multiple energy levels (e.g., a resonator), we can still treat it as a two-level system, regardless of the underlying physics. This is our working assumption throughout the rest of the article.
Note that we can probe the environment of a qubit within a different single-excitation manifold. For example, if we have an anharmonic oscillator as a probe instead of a true qubit, we can populate the second excited state and look for systems coupled to the |e ↔ |f transition. This allows for the calibration of a CPHASE two-qubit gate [18,22,28]. In either case, because we consider the exchange of a single excitation, the effective Hamiltonian remains unchanged.
After a rotating wave approximation, the Hamiltonian of a qubit at the probe frequency f p interacting with a resonance mode at a frequency f RM readŝ where g is the coupling strength of the qubit-resonance mode interaction,σ z,1 (2) are Pauli matrices for the qubit (1) and resonance mode (2), andσ + 1(2) andσ − 1 (2) are rising and lowering operators for the qubit and resonance mode. We solve for the time evolution of the qubit when it is initialized in state |e and with the resonance mode starting in |g . The theoretical probability of finding the qubit in the excited state after a time t is then given byP where (2) is plotted in Fig. 3 (a) as contours. Close to resonance, the excitation swaps between the qubit and the resonance mode with frequency Ω increasing at larger δf , resulting in the familiar chevron pattern. Both the width of the pattern, which we quantify by the full-width half maximum of the amplitude, 4g, and the frequency of the oscillation depend on g, the coupling strength. Crucially, the width depends linearly on g, while the period of the oscillation (and therefore the time location of the first minimum) depend on 1/g. With these observations in mind, we choose to divide the frequency-time space into bins within which we take a constant number n s of swap spectroscopy measurements. Instead of naively sampling the spectrum in a uniform grid, we adapt the measurement based on the coupling strength g we are trying to detect. g determines the time t at which we measure and the bin size. On the one hand, a resonance mode with large coupling strength g has a large width and a short period. For small t then, we choose bins to be wide and short. On the other hand, a more weakly coupled resonance mode appears later in time and has a narrower frequency width and a longer period. In this case, the bins are longer and narrower [see Fig. 3 (a)]. Because we choose the size of these bins to vary by factors of two, we refer to this method as octave sampling.
The goal of octave sampling is to detect a resonance mode by finding the first minimum of an oscillation, where measurements give P e ∼ 0 because the excitation has swapped into the resonance mode. In order to make this division systematic, we introduce the concept of a coupling octave, characterized by the coupling strength g o , where o is the octave number which ranges from 0 to o f [35]. The final octave number o f is determined by the frequency or time resolution that is desired.
For each octave, the full frequency spectrum to be analyzed ranges between a minimum and maximum frequency f min and f max . This range is divided into 2 o bins of equal size, with frequency width ∆f = 2g o and time length ∆t = 1/(4g o ). The location in time of the bins' lower edge is t = 1/(4g o ). One such bin, with g o = g, is highlighted in red in Fig. 3 (a). Note that the highlighted bin is not centered on the oscillation minimum. This is because a low-P e measurement in that area corresponds to a range of possible coupling strengths, namely, those for which g o /2 ≤ g ≤ g o . The resonance mode plotted in Fig. 3 (a) is at the upper end of this range, and, hence, at the lower edge of the bin.
The execution of the algorithm is determined by ∆ = f max − f min . The zeroth octave consists of a single bin which spans the whole spectrum. This width corresponds to a coupling octave g 0 = ∆/2, and therefore to a time, t = 1/4g 0 . For the next octave, we divide the width of the bins by two such that the scan has twice as many bins as for the previous step. The length of the bins in time is correspondingly doubled. An example of this division is shown by the orange grid in Fig. 3 (a).
Explicitly, starting from the zeroth octave o = 0: 1. Divide the frequency range in 2 o bins, each with 3. Take n s swap spectroscopy samples within each of the 2 o bins, sampling uniformly in frequency and inverse uniformly in time, where the notation X ∼ U(a, b) signifies that X is drawn randomly from a continuous uniform distribution U between a and b; 4. Increment o and start over for the next octave.
The total number of bins to be measured depends both on the size of the spectrum ∆ and the final octave number o f . To set o f , we can choose either a maximum time For the same resolution, traditional swap spectroscopy divides the frequency axis in ∆/∆f o f points, and the time axis in t o f /∆t min points, where ∆t min is the time resolution. While octave sampling reaches a time resolution of 1/∆, it would be unfair to compare the traditional method directly to that number. Instead, we assume that 1/∆t min is on the order of 100s of MHz, in order to be able to detect strong couplings, such as those to other qubits or resonators. The total number of points is then

B. Experimental Results
The result of an experimental implementation of the octave sampling algorithm is shown in Fig. 3 (b), which is the efficient version of Fig. 2. The frequency ranges from f min = 4.146 GHz to f max = 5.170 GHz. We choose a final octave number o f = 8. Each bin is color coded according to the average excitation probability P e measured over n s = 5 samples. We are able to discern the same features as in Fig. 2, while acquiring much fewer data. For an equivalent frequency and time resolution, octave sampling requires an order of magnitude less points than traditional swap spectroscopy, making it much more efficient. For a 1024 MHz spectrum divided in 4 MHz, and the ability to detect oscillations at up to 200 MHz in a 250 ns time interval, traditional swap spectroscopy requires 1024/4 × 250/2.5 = 25600 points. The equivalent octave sampling uses 511 bins; if we take 5 samples per bin, as in Fig. 3 (b), we need a total of 2555 measurements.
Given the octave sampling results of Fig. 3 (b), we now intend to determine if there is one or more resonance modes interacting with the qubit. If there are no resonance modes at all, the qubit does not undergo any swap, and we should always measure it to be in |e with P e = 1, notwithstanding T 1 decay. Hence, a measurement of P e < 1 indicates energy loss due to a resonance mode interacting with the qubit.
In practice, however, other spurious experimental effects can lower the measured P e below the theoretical value of one, even in the absence of a resonance mode. Those include, for instance, energy relaxation, qubit state preparation, measurement visibility, bin averaging, and statistical fluctuations. We must therefore set a threshold probability P t < 1 to account for these effects.
Qubit energy relaxation unavoidably makes the qubit decay to |g . In order to avoid classifying this effect as a resonance mode, the threshold should be chosen to be below the T 1 decayed qubit excitation at the measurement time t o f . Explicitly, we must have P t < e −to f /T 1 . The qubit state preparation and measurement visibility affect the maximum qubit excitation that can be measured. The threshold must therefore be chosen below this value. We also need to account for the possibility that some of the n s samples measured within each bin do not fall exactly at the oscillation minimum. Thus, even is there is a TLS visible in a particular bin, it is likely that P e > 0. Consequently, we cannot select an overly low threshold. Finally, since the measurement is a binomial process, we must account for the variance in P e .
In order to cope with these issues, we determine a threshold by inspecting the average excitation of each bin in Fig. 3 (b). We find the bin with the highest average excitation value, in our experiment max(P e ) = 0.976. This quantity is a good estimate for the visibility. We then subtract an empirically selected buffer value of 0.3 to account for the other effects mentioned above and obtain the threshold P t = 0.676.
After setting P t , the average excitation P e for each bin must be compared to it. If for a certain bin P e < P t , we have detected a resonance mode.
Before estimating the coupling parameters, however, we need to verify whether there are multiple resonance modes within a frequency-time scan or not. This task is achieved by considering all bins meeting the condition P e < P t and grouping neighboring bins. If multiple such groups exist, then there are multiple resonance modes interacting with the qubit. In this case, we must split the octave data into multiple frequency-time subscans such that there is only one resonance mode per scan. This procedure is explained in Appendix C. Once the data is split by resonance mode, the parameter estimates are given by searching for the below-threshold bin at the lowest octave.
We performed an experiment where a qubit interacts with three resonance modes, two of which are synthesized with a microwave source and one is an actual resonator. As expected, our estimation method detects three resonance modes. Their estimated frequency and coupling parameters are reported in Table IV. We purposely choose two of the resonance modes to be close in frequency to illustrate an important feature of our estimation method: that two distinct resonance modes are detected as independent modes only if their frequency separation is at least twice as large as the largest coupling strength and twice as large as the minimum frequency width of the bins. The resonator is characterized by a strong coupling coefficient (see Appendix A for the coupling capacitance). For all three modes, the octave algorithm estimates the frequency properly, though the precision is limited by the bin size. The estimated coupling strength is given by a range. Obtaining more precise and accurate results necessitates the online estimation algorithm to be explained in the next section.  Fig. 3

V. ONLINE BAYESIAN LEARNING ALGORITHM
The offline octave sampling algorithm is data efficient and can be performed from a state of zero knowledge of the qubit's spectrum. However, it does not provide a very accurate estimate of the coupling parameters. To improve accuracy, we can use the rough estimate given by the offline method to execute an online Bayesian learning algorithm and refine the parameters in a very short time.
Given an initial probability distribution over the coupling parameters with a resonance mode, the online algorithm successively selects swap spectroscopy measurement settings to increase knowledge. After the result of a measurement is recorded, the distribution is updated according to Bayes' theorem and a new measurement setting is generated. This procedure is repeated iteratively until the distribution converges as desired [43].

A. Theoretical Method
The online estimation algorithm is an experimental implementation of the one proposed in Ref. 36. It employs a particle filter method to efficiently represent the prior and posterior distributions and compute Bayes' theorem at each iteration.
A particle distribution is a discretized representation of a probability distribution. The denser the distribution in a particular region of the parameter space, the higher the probability of those parameters being valid. Here, each particle represents a value for the coupling parameters (f RM , g) of a resonance mode. At the beginning of an iteration, we compute the means, µ(f RM ) = f RM and µ(g) = g , and standard deviations σ(f RM ) = f 2 RM − f RM 2 and σ(g) = g 2 − g 2 of the prior distribution. The next step is to perform a single swap spectroscopy measurement to determine the excited population P e at a particular probe frequency f p and time t. These measurement settings are heuristically selected to increase information gain [36]. In practice, t should scale inversely with σ(g), while f p should be within a factor of µ(g) on either side of µ(f RM ).
We chose the following measurement settings: where a = π/2, c = 5, r 1 is picked from U(−1/2, 1/2) and r 2 is picked from U(0, 1). These parameters are empirical constants determined in Ref. 36, though they were slightly adjusted for this experiment in order to have a larger distribution for f p and t when M > M 0 . Note that, unlike the method proposed in Ref. 36, we choose to limit t to a maximum value well under T 1 . This is done for the same reason as was explained in Subsec. IV A.
For this purpose, we use the hyperbolic tangent function as it has a linear behavior for small arguments, such that tanh(a/σ(g)t max )t max ≈ a/σ(g) when σ(g) is large. After a certain number of iterations M 0 = 25, we modify the heuristic slightly to accelerate convergence. Initially, we choose probe frequencies coarsely according to the value of g. Then, as our knowledge improves, σ(f RM ) decreases and can be used to select frequencies in a narrower range around µ(f RM ). The factor c here is used to avoid choosing measurement frequencies too narrowly. The time t is always weighted by ∼ 1/σ(g), but we bias the selection to larger values after M 0 iterations. The last step in the iteration is to apply Bayes' theorem to update our knowledge of the coupling parameters. We want to obtain the posterior distribution based on the measurement result P e (f p , t). This is achieved in two sub-steps: first, we compute the likelihood of obtaining the measurement value given each particle's (f RM , g) parameters; second, we resample the distribution according to these likelihoods.
We compute the likelihood from the measurement result P e = n e /n, which is the proportion of excited state outcomes we recorded n e for n individual measurement shots. Since we know the theoretical fraction that we should expect to measure,P e :=P e (f p , t, f RM , g), which is given by Eq. 2 in the decoherence-free case [44], we know that the result is a binomial random variable n e ∼ B(n,P e ). Accordingly, the likelihood of obtaining a particular measurement result given measurement settings (f p , t) and a particle (f RM , g) is the probability mass function L(n e |f p , t, f RM , g) = n n e P e ne (1 −P e ) (n−ne) , where n ne is the binomial coefficient. In effect, we are computing the probability that the measurement result obtained corresponds to a resonance mode with coupling parameters (f RM , g). The next step is to resample the distribution to keep only those parameters that are most probable. Although this task can We apply Bayes' theorem to determine the posterior distribution. This task is achieved by resampling the particles according to their likelihoods. The distribution is split in two "clouds". After resampling, the posterior distribution can be used as the next iteration's prior.
be achieved in a variety of ways, the general idea is to pick particles from the prior at random, weighted by the likelihood. To avoid duplicate particles in the posterior distribution, we add normally distributed random noise proportional to the covariance of the prior. The exact procedure used in this work is described in Appendix D.
The iteration process is visualized in Fig. 4, allowing us to understand more intuitively how the particle filter technique works. If the measurement is useful, i.e., the resulting likelihood favors a subset of the prior, the posterior distribution is shrunk or filtered, improving knowledge of the parameters. Otherwise, if the likelihood does not discriminate the particles, the distribution is not modified significantly. After resampling, the next iteration starts with the last iteration's posterior as prior.
Compared to octave sampling, the task of the online Bayesian learning algorithm is simpler since we must already have knowledge of the presence of an interaction. The particle filter can therefore fit to the parameters that are the most likely given the measurements. At the end of the final iteration, the parameters and their uncertainties are given by the mean and standard deviation of the particle distribution. If the algorithm converges, the final particle "cloud" is small, resulting in an estimate that is accurate: both true and precise. If the algorithm does not converge, meaning that the final particle cloud is not tightly concentrated in a single region, running the experiment again might be necessary. Regenerating the initial particle distribution via the octave sampling method could also improve detection performance.
Note that the algorithm does not take into account any potential errors on the value of the measurement settings. Therefore, directly taking the standard deviation of the distribution to be the uncertainty on the estimated parameters is valid only in the case where there is no uncer-tainty on the measurement settings. This will generally not be the case in an experimental setting.

B. Experimental Results
We run the online Bayesian learning algorithm on the three distributions generated from the octave data for each detected resonance mode. The generation of those distributions is discussed in Appendix C. For each mode, we perform 35 iterations of the algorithm, at which point the distribution has converged. The runtime of the algorithm for a single resonance mode is approximately 23 s. We then compute the mean of the distribution after the final iteration. The results are shown in Table II. As expected, the parameters of the synthesized modes and of the measurement resonator are correctly identified. Standard errors are not shown because the standard deviation of the particle distribution after the final iteration converges to a small value that is not representative of the true error. The error must instead be estimated according to other potential limitations. For our experiment, this should be done according to the characteristics of the flux pulse used to set f p and t. More details can be found in Appendix B.
To test the performance of the estimation algorithm, we run it 1000 times on the second resonance mode (RM2) with slightly different initial distributions. We plot the convergence of the parameters in Fig. 5. As shown by the histograms, more than 99% of the runs converge successfully to properly estimate the frequency and coupling strength, with just a few failures. The average of the parameters after the 35 th iteration is f RM = 4.8301(4) GHz and g = 1.45(9) MHz.
We note that one possible cause of failure is the overes-  timation of g by an integer multiple; under these conditions, in fact, crests in the oscillations partially overlap.
Other possible causes include variance in the parameters of the qubit itself (our probe). If the frequency of the qubit changes abruptly, for example due to appearance of a semi-resonantly coupled TLS defect, the fluxto-frequency calibration becomes invalid; this results in inaccurate measurement settings [33,41,42].
The mean of each of the 1000 initial particle distributions is distributed uniformly at random within a 10 MHz and 1.5 MHz range for f RM and g, respectively. Each individual particle distribution is uniform, with a width of 15 MHz in f RM and 2.5 MHz in g.

VI. DISCUSSION
Both algorithms presented above depend on a few parameters that are crucial to their function. For the octave sampling algorithm, the choice of the frequency range to be measured is naturally determined by the properties of the device: superconducting qubits have a limited frequency range within which they operate optimally. For the device in this work, the upper end of the measurement range f max simply corresponds to the maximum attainable frequency. The lower limit f min is chosen to be as low as is possible while still being able to control and measure the qubit. Since we use a resonator for readout, the farther away the qubit is in frequency from the resonator, the lower the fidelity of the measurement. Hence, a qubit cannot operate far away from its readout resonator. Other constraints, e.g., pulse control bandwidth, might dictate even tighter limits.
A second important parameter for the octave sampling method is the number of measurements made per bin n s . In principle, a single high-quality (many shots) measurement of P e at the center of the bin should be enough. This would also make for a fair comparison with the traditional swap spectroscopy technique. However, because of the efficiency of the octave method, we can afford to take a few more measurements per bin. This is what we have chosen to do, by randomly distributing n s = 5 measurements per bin. Another way to take more than one sample could be to sample within bins regularly. In either case, this redundancy serves both to marginally increase the detection sensitivity (by increasing the resolution), or to protect against possible statistical fluctuations in the measurement.
The last parameter we discuss is the number of octaves to be measured, which corresponds to the maximum time t o f , or equivalently, the final frequency resolution ∆f o f , as explained in Sec. IV A. Generally, this parameter should be determined by the requirements of the experiment for which the calibration is made. If a long gate sequence is needed, e.g., for randomized benchmarking, detecting weakly coupled RMs is important. This would not necessarily be the case for shorter experiments, like process tomography. A frequency resolution ∆f o f ∼ 1 texp , where t exp is the length of the experimental gate sequence, is therefore generally a good choice.
For the particle filter algorithm, the choice of a, c, and M 0 is discussed in Ref. 36. Other parameters of interest include t max , and the number of particles to be used. The time t max is used in Eq. 4 to restrict the maximum measurement time t. This is necessary because the qubit eventually decays to the ground state. To obtain reliable results, t max should be set well below T 1 . Note that another way to limit the maximum measurement time would be to replace tanh(a/σ(g)t max )t max with a/σ(g) in Eq. 4 (as in the original proposal), and simply stop the algorithm once a sufficiently small σ(g) is reached.
The number of particles to be used is constrained mainly by the performance of the computer running the resampling procedure and, potentially, numerical accuracy issues [45]. As a rule of thumb, at least 10000 particles should be used; in this work, we have used 40000.
In Sec. IV A, we explained that our method can be used not only for a simple qubit swap spectroscopy experiment, but also with a double-excitation protocol, to look for resonances with the |e ↔ |f transition. The algorithms discussed in this work are, in fact, more general than that. As long as the qubit dynamics look like a chevron pattern, we can use the same methods to efficiently detect the location and estimate the parameters of the interaction. This is the case, for example, with a whole class of parametric two-qubit gates, where instead of varying the probe frequency f p of the qubit, we vary the frequency of a radio-frequency flux drive applied to the DC SQUID of a qubit or tunable coupler [46][47][48].
Finally, we briefly discuss the problem of choosing qubit operating frequencies. Indeed, once the calibration showcased in this work has been accomplished and all resonant couplings have been identified, the goal is then to make use of the information to optimize performance of a quantum processor. This process is architecture dependent. For an array of directly-coupled superconducting qubits, we want to avoid crosstalk between neighbor qubits, and minimize interactions with TLS defects. We therefore need to choose the idle fre-quencies of all qubits at the same time, taking into account both wanted and unwanted couplings. Additional concerns apply for choosing the operating frequencies of two-qubit gates: we must consider account the frequency path that the qubits will take during the gate. Having a qubit cross through a resonance with a TLS, for example, is not desirable.
While this work does not explain the process needed to perform this performance optimization (for that see, e.g. Ref. 49), we emphasize that the runtime improvements of the offline and online algorithms when compared to traditional swap spectroscopy enables several advantages. First, the calibration may be run more often. Second, the calibration is affordable enough to be run on a larger spectrum, and so potentially more resonance modes, giving the frequency optimization process more information to work with.

VII. CONCLUSION
In conclusion, we explain two methods for the Hamiltonian parameter estimation of resonant couplings in the context of tunable superconducting qubits. The octave sampling technique can be run without knowledge about the environment of the qubit, and allows efficient detection of coupled resonance modes. The parameter estimation algorithm can be performed or omitted depending on whether a more accurate knowledge of the coupling parameters is desired. Using these algorithms reduces the number of measurements needed by a square-root factor when compared to the traditional method. This translates to a reduction in runtime by one order of magnitude in typical conditions.
We experimentally demonstrate both techniques on an Xmon transmon superconducting qubit and evaluate their performance. We are able to detect the resonance with the qubit's measurement resonator, as well as with synthesized resonance modes. Overall, we determine that the methods are efficient, reliable, and readily automated. We expect this type of calibration to be critical to the operation of large-scale quantum computers, superconducting and not. Future work includes integrating the information we acquire by our methods into a comprehensive optimization process for selecting the operating frequency of each qubit in a quantum processor.

ACKNOWLEDGMENTS
This research was undertaken thanks in part to funding from the Canada First Research Excellence Fund (CFREF) as well as the Discovery and Research Tools and Instruments Grant Programs of the Natural Sciences and Engineering Research Council of Canada (NSERC). The device used was fabricated at the Quantum-Nano Fabrication and Characterization Facility at the University of Waterloo, Canada. We would like to acknowledge the Canadian Microelectronics Corporation (CMC) Microsystems for the provision of products and services that facilitated this research, including CAD tools and design methodology. The authors thank Christopher Warren for useful discussions. Performing swap spectroscopy requires the ability to set the qubit's frequency to a desired value f p for a particular duration. In our experiment, this is done with a flux pulse applied to the DC SQUID with an AWG. As shown in Fig. 8, the flux pulse reaches the sample after going through multiple stages of filtering, attenuation, and connections. This means that the waveform will be modified compared to what is generated by the AWG.
In addition, while we control the amplitude of the pulse, we are ultimately interested in the resulting frequency of the qubit. We therefore need a way to convert between the amplitude of the flux pulse A and the qubit frequency f p . This can be done, for example, with pulse spectroscopy, where we send π-pulses to the qubit at different frequencies while it is detuned by a flux pulse of a particular amplitude. The qubit frequency for that amplitude can then be fit. This is repeated for many amplitudes in order to get a map between A and f p .
The above considerations mean that the measurement settings (f p , t) that we select may contain multiple kinds of potential errors. This must be taken into account when estimating the coupling parameters with either the offline or online algorithms.
For example, if f p is higher than the true probe frequency of the qubit due to some systematic error, the result f RM reported by the online algorithm will be higher than the true value as well. Similarly, an error on the value of t will lead to a wrong estimation of g. In practice, this kind of systematic error is not a major problem as long as the error is consistent between experiments.

Appendix C: Details on Octave Analysis
We explain here the method employed to (1) split the raw octave data into spectra containing a single RM, if any, and (2) generate a particle distribution from spectra containing a single RM.
As explained in the main text, we want to make sure that there is only a single resonance in the data when we estimate a frequency and coupling strength. If there are no resonance, or if there is more than one, the particle distribution that we generate is nonsensical. We therefore need to check for bins below the threshold, and split the data accordingly.
To identify resonances, we start from the deepest octave (the one corresponding to the longest time), and look for bins with average excitation below the threshold. If there are multiple consecutive bins below the threshold, we keep only the one with the lowest excitation.
This gives us a list of independent peaks found in the data. We then go to the next deepest octave and perform the same thresholding. At this point it is possible that, if any resonance was found in the previous octave, a bin corresponding to the same resonance is found as well. If that is the case, we keep only the lower octave (shorter time) bin information.
We continue this procedure until the zeroth octave, at which point we are left with the smallest octave (the lowest time) bin for every resonance detected. To split the spectrum, we simply cut the data at the midpoint of each detected resonance. If this cut happens to fall within a bin, we slice the bin in two, making sure to keep track of each part's proportion. For the data in Fig. 3 (b), we detected three resonances, tabulated in Table I.
Following this procedure, we can restrict ourselves to a frequency-time scan containing a single resonance mode. We can specify a discrete distribution representing our knowledge of the frequency and coupling parameters between the qubit and the resonance mode. We thus choose all the bins below threshold and weigh them following a heuristic procedure to find how much below threshold they are. This allows us to draw a particle distribution associated with the bins as explained in the following.
Each particle represents a frequency-couplingstrength pair. We choose a bin according to its weight and generate a particle within it. In this experiment, the weight of each bin is chosen to be proportional to max(P e ) − P e . Thus, bins for which the average excitation P e is low have a high weight. The frequency of the generated particle is chosen uniformly at random within the bin bounds and its coupling strength is selected from U (∆ o , 2∆ o ), where ∆ o is the bin width. We then repeat this procedure to generate a full distribution.
The three resulting distributions can be seen in Fig. 7, with their statistics tabulated in Table IV. We use a distribution comprising 40000 particles, though this number can be adjusted depending on the capabilities of the computer controlling the experiment. These particle distributions can be used as the starting point for the online algorithm discussed in Sec. V. TABLE IV. Frequency and coupling parameters for three resonance modes computed from the data in Fig. 3  3 (2) 2(1) 4 (7) 4. FIG. 7. Initial particle distribution generated from the octave data in Fig. 3 (b) after splitting the spectrum in three parts.

Appendix D: Bayesian Resampling Procedure
The exact resampling procedure is described in detail in Chap. 2 of Ref. 35. We reproduce a shortened version here for completeness.
As inputs, we require the particle locations, as an array of frequency-coupling-strength pairs, in addition to their likelihoods, as computed in the main article. Note that the likelihoods must be normalized to sum to unity, at which point they become weights. This normalization ensures that the weights array is a valid discrete probability distribution that can be drawn from.
The algorithm then draws particles from the input distribution according to the weights, and, from the position of the particles drawn, generate new ones by adding "noise". This prevents having duplicates in the output distribution, even if a particle from the input is drawn multiple times. The resampling algorithm is shown in pseudocode in Alg. 1.
Algorithm 1: Particle Resampling Input: Array of particle positions { x k } Input: Array of particle weights {w k } Output: New particle positions { y k } Function resample({ x k },{w k }) a = 0.98 µ = mean({ x k }) Σ = (1 − a 2 )cov({ x k }) for i ∈ 1 : n do l = rand(Discrete({w k })) µ l = a x l + (1 − a) µ y i = rand(Normal( µ l , Σ)) end return { y k } end Note that the particle at index i in the { x k } array has its corresponding weight at the same index in the {w k } array. In addition, the amount of noise added to the position of the particle drawn is controlled by the spread of the input distribution, quantified by taking the covariance. For more details, see Ref. 35. The sample is mounted in an aluminum package at the mixing chamber stage of a dilution refrigerator, at a temperature of approximately 10 mK. We use a dedicated microwave line to drive the qubit. Two coaxial lines are joined at the mixing chamber for the flux biasing of the DC SQUID. We use a homodyne readout setup, with the in-phase quadrature (IQ) mixers configured for image rejection.