Automated analysis of single-tone spectroscopic data for cQED systems

TL;DR

This paper introduces an automated, robust, and fast tool for analyzing single-tone spectroscopy data in cQED systems, enabling precise calibration of quantum devices through maximum likelihood estimation.

Contribution

The authors developed an open-source Python tool that automates the analysis of STS data for cQED systems, improving calibration efficiency and accuracy.

Findings

The tool accurately extracts physical parameters from noisy STS data.

It significantly reduces calibration time for transmon qubits.

The approach is validated on real cQED devices.

Abstract

Physical systems for quantum computation require calibration of the control parameters based on their physical characteristics by performing a chain of experiments that gather most precise information about the given device. It follows that there is a need for automated data acquisition and interpretation. In this work, we have developed a tool that allows for automatic analysis of single-tone spectroscopy (STS) results for a single cell consisting of qubit and resonator in the circuit quantum electrodynamics (cQED) architecture. Using analytic approaches and maximum likelihood estimation, our algorithm is capable of finding all relevant physical characteristics of the cell by using only the measured STS data. The described approach is fast and robust to noise, and its open-source Python implementation can readily be used to calibrate transmon qubits coupled to notch-port resonators.

Figures

Tables

Table 1. Table 1: Parameters of the model and their physical meaning.

Parameter	Physical meaning
$f_c$, GHz	bare cavity frequency
$g$, GHz	cavity-qubit coupling constant
$\Pi$, A	period in current of the transmon
$I_{ss}$, A	transmon upper sweet-spot current
$f_{ge}^{max}$, GHz	transmon frequency at sweet-spot
$d$	transmon asymmetry parameter

Table 2. Table 2: Grid specifications for the brute force algorithm for STS detection.

Parameter	Value range	Steps
$f_c$	${⟨f_r⟩}_I\pm 1$ MHz	3
$g$	20 - 40 MHz	5
$f_{ge}^{max}$	4 - 12 GHz	80
$d$	0 - 0.9	9

Table 3. Table 3: Analysis of optimal parameters found for the three cases in Fig. 5 after the brute and Nelder-Mead (N-M) optimizations. min ⁡ σ 𝜎 \min\sigma columns stand for the lower bounds of the standard deviations of the MLE estimators found using the Fisher information matrix at the optimum.

Parameter	(a)				(b)				(c)
	Brute	N-M	$\min \sigma$	TTS	Brute	N-M	$\min \sigma$	TTS	Brute	N-M	$\min \sigma$	TTS
$f_c$, GHz	6.5003	6.5007	4$\times {10}^{-6}$	n/a	6.962	6.9631	1.6 $\times {10}^{-4}$	n/a	6.47	6.465	2$\times {10}^{-4}$	n/a
$g$, MHz	24	35.8	0.17	n/a	36	64.9	15	n/a	36	86.1	1	n/a
$f_{ge}^{max}$, GHz	8	8.97	0.013	9.04	10.3	11.3	1.35	9.08	6.3	5.89	0.01	5.9
$d$	0.5	0.09	0.013	0.13	0.5	0.49	0.07	0.6	0.1	0.25	0.05	0.3
Loss, kHz	251	20			338	51			2038	149

Table 4. Table 4: Performance of the algorithm on a 2.5 GHz 2-core i5-5337U CPU for the case (a). Most of the time is taken by the f r subscript 𝑓 𝑟 f_{r} extraction and the brute search steps.

	Extract $f_r$	$\Pi$, $I_{ss}$	Brute	N-M	Total
Time, s	2.72	0.3	2.79	1.37	7.34
Ratio, %	37	4	38	18	100

Equations

$$f_{\pm }(I)\equiv f_r(I)=\frac{f_c+f_{ge}(I)}{2}\pm \sqrt{g^2+(f_{ge}(I)-f_c{)}^2/4},$$

$$f_{ge}(I)=f_{ge}^{max}{\left[{\cos }^2\left(\frac{\pi (I-I_{ss})}{\Pi }\right)+d^2{\sin }^2\left(\frac{\pi (I-I_{ss})}{\Pi }\right)\right]}^{\frac{1}{4}},$$

$$S_{21}^{notch}(f_p)=a{e}^{i\alpha }{e}^{2\pi i{f}_p\tau }\left[1-\frac{Q_l/Q_e^{\prime }}{1+2i{Q}_l(f_p/f_r-1)}\right],$$

$$I_{ss}=\varphi +\Pi (1+D)/2.$$

$$I_{ss}=\varphi +\Pi D/2.$$

$$\mathcal{L}$$

$$\mathcal{M}$$

$$f_{\pm }=f_{\pm }(I_i,\Pi ,I_{ss},f_c,g,f_{ge}^{max},d)$$

$$r/\sigma =r/\sqrt{{\sigma }_0^2+{\sigma }_1^2},$$

$${\hat{H}}_{tr}=4E_C{\hat{n}}^2-E_J({\Phi }_e)\cos \hat{\phi },$$

$$E_J({\Phi }_e)$$

$$k(d,{\Phi }_e)$$

$$E_m=m\sqrt{8E_J({\Phi }_e)E_C}-\frac{E_C}{12}(6m^2+6m+3).$$

$$f_{ge}({\Phi }_e)\approx f_{ge}^{max}\cdot \sqrt{k(d,{\Phi }_e)}$$

$$\left[\sqrt{k(d,{\Phi }_e)}-1\right]\cdot E_C$$

$$\hat{H}/h=\frac{f_{ge}}{2}{\hat{\sigma }}_z+f_c{\hat{a}}^{†}\hat{a}+g({\hat{\sigma }}^{-}{\hat{a}}^{†}+{\hat{\sigma }}^{+}\hat{a}),$$

$$E_{g,0}/h$$

$$E_{\pm ,n}/h$$

$$g({\Phi }_e)\propto |⟨g|\hat{n}|e⟩|\approx \frac{1}{\sqrt{2}}{\left(\frac{E_J({\Phi }_e)}{8E_C}\right)}^{1/4}.$$

$$\mathcal{L}(\mathbf{p}|\mathbf{w})=\prod _{i=1}^N\frac{1}{\sqrt{2\pi }\sigma }\exp [-(f_{r,i}-{\mu }_i{)}^2/2{\sigma }^2]$$

$$-\ln \mathcal{L}(\mathbf{p}|\mathbf{w})$$

$$\equiv {\chi }^2/2+N\ln \sqrt{2\pi }\sigma$$

$${\text{Var}}_{\mathbf{w}}[{\mathbf{w}}^{\ast }]\ge \text{diag}[\hat{I}(\mathbf{w}{)}^{-1}],\forall \mathbf{w}.$$

$$\hat{I}(\mathbf{w})$$

$$\displaystyle=E_{\mathbf{w}}\left[-\mathbb{H}_{\mathbf{w}}\ln\mathcal{L}\right],\$$

$${\mathbb{H}}_{\mathbf{w}}=\left(\begin{array}{ccc}{\partial }^2/\partial w_1^2& {\partial }^2/\partial w_1\partial w_2& \dots \\ {\partial }^2/\partial w_1\partial w_2& {\partial }^2/\partial w_2^2& \dots \\ ⋮& ⋮& \ddots \end{array}\right)$$

$$\hat{I}(\mathbf{w})=\sum _i\frac{{\mathbb{H}}_{\mathbf{w}}{\mu }_i^2(I_i,\mathbf{w})-2{\mu }_i(I_i,\mathbf{w}){\mathbb{H}}_{\mathbf{w}}{\mu }_i(I_i,\mathbf{w})}{2{\sigma }^2}.$$

Full text

††preprint: AIP/123-QED

I Introduction

Computation on a quantum computer involves operating large numbers of physical quantum bits (qubits). One of the most significant challenges here is that individual man-made qubits are not identical and cannot all be operated same wayKelly et al. (2018); Chen (2018). Typically, control parameters for a device depend on its physical parameters and are determined by a sequence of calibration steps. Since the devices cannot be produced identical, the optimal parameters vary among them and should be found in situ. Manual adjustment of these parameters is not a scalable solution, and operating just a handful of qubits already requires a computer system to handle the calibration on its own. Moreover, a lot of research is still being done to optimize chip designs and fabrication methods. This research includes gathering statistics on many different samples that should be measured in a comparable way and consistently evaluated. An automated measurement system does not only speed up this process, but also excludes human errors.

In this work, we are proposing and implementing several methods of data processing and computer vision that aid automatic calibration of circuit quantum electrodynamics (cQED)Blais et al. (2007) cells used to read out the quantum states of superconducting qubits. Computer vision is understood here as an enterprise that employs statistical methods to disentangle data using models constructed with the aid of geometry, physics, and learning theory.Forsyth and Ponce (2011) As it will be shown below, the outputs of circuit calibration on certain steps can to be represented as images, and it is hardly possible to automatically extract the required information from them without resorting to methods of computer vision. The approach that we propose is accurate, fast and robust to noise, it is applicable to superconducting qubits with a spectrum periodic in magnetic (or electric) field. The generic single cell here is described with closed-form model equations, and is compatible with the paradigm of one readout resonator per one qubit.Versluis et al. (2017); Kelly et al. (2015)

Within the considered paradigm, the flow of an experiment that characterizes a single qubit-resonator cell of a quantum processor chip is depicted schematically in Fig. 1(a). It is similar to the procedures for measuring cQED samples of qubits with tunable frequency described in Refs.[Jerger, 2013; Chen, 2018]. For non-tunable qubits, a different algorithm of automatic calibration was described in the patent application Ref. [Bloom et al., 2018].

We start from briefly discussing each step shown in Fig. 1(a). Usually, one starts by locating the dip in the transmission corresponding to the resonance frequency of the readout resonator. Having found the dip, the experimenter can set the range to scan the probe frequency $f_{p}$. Scanning this range during the next step produces single-tone spectroscopy (STS) data. STS tracks the frequency shift of the resonance dip while changing the control parameter (e.g., the magnetic field applied to the sample that is proportional to some control current $I$). By processing the results of the STS, one can get a coarse estimate on the qubit frequency dependence on current $I$ denoted as $f^{(0)}_{ge}(I)$ (here $g$ and $e$ stand for ground and excited state of the qubit, respectively). Besides that, during this step one also gets accurate estimates of the resonator frequency $f_{r}(I)$ and the resonator-qubit coupling strength $g$. Based on these parameters, one can later set the scan ranges for the excitation frequency and the current ($f_{exc}$ and $I$) to perform the two-tone spectroscopy (TTS). Finally, by processing the TTS results, one can obtain the precise dependence of the qubit transition frequency $ge$ on current $f_{ge}(I)$ with accuracy sufficient for pulsed experiments in order to perform accurate single-qubit gates. Overall, the process can be structured so that the results of one measurement step define the parameters of the next, successively acquiring more accurate information about the physical system.

As one can see, the results that are obtained in this scheme may be divided in two groups: some contain single curves (data are 1D manifolds), and others contain colour scale images (heatmaps, 2D manifolds). The work presented in this paper is confined to the automatic analysis of the latter. More precisely, we are interested in processing STS data. In Fig. 1(b) we illustrate three types of correlated measured data and the theoretical model. Depending on the range of the qubit transition frequencies relative to the resonator frequency, we need to interpret STS data as either a single curve or two curves defined by the cQED model. Here we are aiming at building a method that is able to distinguish between these three cases and to estimate quantitatively the physical parameters of the system based on the best-fit model.

The model to fit the data is based on the theoretical description of the cQED systems and physical qubits. In Appendices A and B we derive all the necessary equations to form the model curves that are expected to appear in single-tone spectroscopy. We use a transmonKoch et al. (2007) with an asymmetric SQUID (superconducting quantum interference device) as a qubit. Since our methods do not depend on the particular shape of the magnetic (or electric) field dependence of the qubit energy levels, without any loss of generality the presented approach can also be used for other types of qubits.

For the singe-tone spectroscopy, the theoretical curves are described by the function (see Fig. 1(b)):

[TABLE]

where $f_{c}$ stands for the field independent resonator frequency; $f_{r}(I)=f_{c}$ when $g=0$ (it is also defined as the bare cavity frequency, or the frequency of the cavity with no qubit coupled to it). Next, the qubit transition $f_{ge}(I)$ is expressed as

[TABLE]

where $\Pi$ is the period of the spectrum in current $I$, $d$ is the SQUID asymmetry parameter and $f_{ge}^{max}$ is the qubit frequency at the “sweet spot” corresponding to current $I=I_{ss}$, where the applied magnetic field exactly compensates the non-zero local stray field that is usually present on the chip in a real experiment.

We note that for the avoided crossings pattern both curves from (1) are necessary. Alternatively, for the cases when the qubit frequency is entirely below (above) the bare cavity frequency only $f_{+}(I)\ \left(f_{-}(I)\right)$ is used.

Overall, we thus have six parameters for our model: $f_{c}$, $g$, $\Pi$, $I_{ss}$, $f_{ge}^{max}$ and $d$. Since they are the parameters in the Hamiltonian for the compound cavity-qubit system described in Appendix A, hereafter we call them the Hamiltonian parameters.

The structure of the rest of the paper is as follows. First, we describe our approach to extract the Hamiltonian parameters from the STS data. Then we address issues related to the accuracy, performance, and reliability of the algorithm. Finally, we discuss the limitations of the algorithm, further applications and future work.

II Methods

Before describing the methods, we first highlight important peculiarities of the data and the essential experimental details. A sample in the cQED architecture interacts with incident microwave radiation. Therefore, the standard mathematical concept to describe its properties is the scattering matrix $S_{ij}$ which shows how harmonic signals scatter between $i$th and $j$th ports of the sample. The absolute value of the complex scattering parameter $S_{21}(f_{p})$ shows the amplitude of the signal transmission between the ports 1 and 2 at some probe frequency $f_{p}$. Alternatively, $\angle S_{11}(f_{p})$ would show the change of the signal phase after the reflection from the port 1. By varying the probe frequency and other control parameters, one can record changes in $S_{ij}$ and from that draw conclusions about the properties of the studied device.

As an example, in Fig. 2(a) we show a heatmap for the amplitude $|S_{21}|$ obtained in STS for a notch-type (side-coupled) coplanar waveguide resonator connected to a tunable transmon qubit. $|S_{21}|$ is shown in colour, and two other axes show the probe frequencies $f_{p}$ and the coil current $I$ values for which it was recorded. Looking at the data slices in Fig. 2(b), one can say that the plot shows the dependence of the resonator dip frequency on the applied coil current.

While in Fig. 2(a) we have shown only the amplitude $|S_{21}|$, the phase of the signal $\angle S_{21}$ is also important. The amplitude has a local minimum at resonance, while the phase has a local maximum of its derivative there. For certain kinds of readout resonators the phase provides better contrast than the amplitude and should be preferred. This means that, in general, the interpretation of the data requires more advanced analysis than simply identifying the location of the minimal transmission amplitude.

Now, knowing the internal structure of the data, we outline our algorithm of extracting the Hamiltonian parameters from them. The main idea is to fit the model (1) to the data in Fig. 2(a) in the sense of maximum likelihood estimationBishop (2006) (MLE, see also Appendix C). To reduce the complexity of the estimation in terms of the number of parameters, we first perform the MLE of the resonator frequency by fitting only along the $f_{p}$ axis using an external fitting library called circlefitProbst et al. (2015) (see Section II.1 for details). This removes the extra dimension from the data which now can be compared directly with the model leaving only 6 fit parameters. With an assumption111We note that the normality of noise in estimated $f_{r}$ is asymptoticJennrich (1969) and would only be attained if the number of data points during the resonator fit was infinite. Therefore, this assumption is in fact an approximationAnastasiou, Reinert et al. (2017). It is not essential for our method to be statistically consistent,Jennrich (1969) but allows to explicitly find the likelihood function and calculate the Cramer-Rao lower bounds for the fit parameters variances. that the distribution of noise in the estimated $f_{r}$ is normal, the MLE is equivalent to a nonlinear least-squares problem with global optimization. It is a complex task, firstly, due to the periodic dependence of the frequencies on current $I$ and the unknown position of the qubit sweet spot $I_{ss}$ and, secondly, due to the strong nonlinear dependence of Eqs. (1) and (2) on other parameters. Fortunately, it is possible to get a good initial guesses for $\Pi$ and $I_{ss}$ and then do a brute force search for the solution in other variables due to the well-defined bounds for each of them. Finally, we can fine-tune the result in all six parameters using a local optimizer.

The detailed description of the method which uses Fig. 2(a) as an example is presented below. As it has already been mentioned, the chosen type of the resonator-qubit arrangement which yields the avoided crossings pattern is not unique; however, the other two cases are treated the same way.

II.1 Extracting resonance frequency

The resonator curve MLE is based on the circlefit library which does the least-squares fitting of various types of microwave resonators. For each current $I$, it fits the complex transmission $S_{21}(f_{p})$ which, for a single ideal resonator, is a circle on the complex plane. However, when an real experimental device is considered, the response model becomes more complex. For example, for a notch-type resonator used in this work the following model is employed:Probst et al. (2015)

[TABLE]

where $a$ denotes the overall attenuation (amplification) in the transmission lines connecting the resonator to the measurement apparatus, $\alpha$ is the frequency-independent phase shift in the transmission lines, $\tau$ is the propagation delay far from the resonance; $Q_{l}$, $Q_{e}^{\prime}$ are the usual loaded and complexKhalil et al. (2012) external quality factors, respectively; $f_{r}$ is the resonance frequency. Fitting procedure of the $circlefit$ library is to sequentially extract all of the model parameters by a succession of partial approximations (see Probst et al. (2015) for details). For brevity reasons, we do not present here the full complex fits for our data; however, as an illustration in Fig. 2(b) one can see the absolute values of the optimized models (solid black lines) that match well with the absolute values of the experimental data taken from slices 1 and 3 of Fig. 2(a). Finally, from the optimal model parameters we obtain the resonance frequency $f_{r}(I)$ depending on current. A possible caveat is that for some $I$ values the resonance dip may disappear so the resonator fit would fail (this happens, e.g., at the avoided level crossings; see Fig. 2(b), slice 2). Therefore, such slices are excluded in advance via a threshold condition.

The resulting plot of the extracted resonator frequency versus current $I$ is shown in Fig. 2(c) (blue dots). There are some current values located between the branches where the data points are missing, as expected, due to the absence of the resonance. Additionally, we plot here the mean value of the detected frequencies shown as a dashed black line. This parameter is important since it will, firstly, serve as an initial guess for the cavity frequency of the model (1) with $f_{c}\approx\langle f_{r}\rangle_{I}$ and, secondly, will be used in the period and phase extraction algorithm which tracks the changes of the sign of the value $\Delta f_{r}=f_{r}-\langle f_{r}\rangle_{I}$ (marked as grey dashed lines in Fig. 2(c)). Locating these sign changes allows one to find the qubit sweet spot without fitting the full model.

II.2 Extracting field period and sweet spot locations

As we already have stated above, one of the obstacles for the fitting is the periodicity of the data in one of the fitting parameters. Along with the equally unknown phase of the signal, this leads to the presence of many local minima in the loss function which impede the progress of iterative optimization algorithms. In other words, the unknown parameters $\Pi$ and $I_{ss}$ are preventing us from finding the global minimum. So it would be very convenient to determine the period and the sweet spot location without fitting the full model (1) to alleviate the stated problem. In the absence of noise, this would be a trivial task. Indeed, one would just need to find the sign changes in $\Delta f_{r}$ directly by traversing the array of its values, and then compute the period as the distance between two similar sign changes, and the sweet spot location as the point between the adjacent ones. However, in the conditions of experiment, the noise to some extent is always present, therefore below we suggest a solution that is insensitive to local sporadic perturbations of the data.

An obvious method to find analytically the period and the phase (here, the phase is equivalent to $I_{ss}$) of a signal is the Fourier transform. Fast Fourier transform (FFT) should be applied to the data, and then the peak of the largest magnitude in the spectrum would give the likely frequency and the phase of the signal. FFT works well for the continuous curves, but it becomes inconvenient for the avoided crossings pattern due to the complexity of its spectrum caused by the discontinuity. Moreover, FFT may be not sufficiently accurate when the data contains only a couple of periods. A more robust tool to find the period in a given dataset $y$, especially when it contains just a few periods, is the autocorrelation function $\mathcal{R}_{yy}(l)=\sum_{n}y_{n}y_{n-l}$. The location of the largest of its local maxima (except for $l=0$) equals exactly the sought periodParthasarathy, Mehta, and Srinivasan (2006). However, this is only true when the mean value of the function is zero; otherwise, $\mathcal{R}_{yy}(l)$ would be linear with a steep slope that may smear out all the extrema. Therefore, instead of directly calculating the autocorrelation of $f_{r}(I)$ we consider the standardized function $\Delta f_{r}(I)$ introduced above which has zero mean and the same period. For this function, the autocorrelation function $\mathcal{R}_{\Delta f_{r}\Delta f_{r}}(\Delta I)$ depending on the current offset $\Delta I$ is shown in Fig. 3(a). As one can see, $\Delta I$ spans 200 $\mu$A just as the data itself. This means that to calculate $\mathcal{R}_{\Delta f_{r}\Delta f_{r}}(\Delta I)$ the data is being zero-padded at all $\Delta I$ except for $\Delta I=0$, and this is why we get diminishing correlation peaks at $\Pi$, $2\Pi$, etc. The orange dot in the plot Fig. 3(a) shows the highest local extremum of $\mathcal{R}_{\Delta f_{r}\Delta f_{r}}(\Delta I)$ at 88 $\mu$A. It is a very prominent peak and can be easily distinguished among all others. There is as well a small peak at 176 $\mu$A which corresponds to $\Delta I=2\Pi$. Note that the autocorrelation function is not ideally smooth and has some abrupt bends on the sides of the peaks. This happens because of some missing points in the resonator frequency data (corresponding to the plateaus of Fig. 2(b)) that were added back as zeros since the step between adjacent data points should be identical for the autocorrelation method to be accurate. This zero padding may be observed directly in Fig. 3(b).

The problem with the autocorrelation function is that it does not tell us anything about the phase of the signal, i.e. the exact value of the sweet spot current. Therefore, after having found the period, we need another method to precisely determine $I_{ss}$. This method consists of finding the global maximum of the zero-lag correlation function $\mathcal{R}_{\Delta f_{r}S}(0)$ between $\Delta f_{r}(I)$ and a square wave $S(I,\Pi,\phi,D)$ having the same period but unknown phase $\phi$ and duty cycle $D$ which is the relation between the horizontal spans of the “high” and “low” levels. The values of the “high” and “low” levels must be opposite, i.e. 1 and -1, while their magnitude does not matter. An illustration of a square wave function satisfying the optimal condition is presented in Fig. 3(b) in orange. The phase $\phi$ (in $\mu$A) denotes the $I$-coordinate of the first point after the rising edge, and is marked with a triangle. Generally, the idea behind this approach is to robustly detect sign changes that were shown back in Fig. 2(c). Since $\mathcal{R}_{\Delta f_{r}S}(0)=\sum_{I}\Delta f_{r}(I)\cdot S(I)$, it gains value if $\Delta f_{r}$ and $S$ are of the same sign for a certain $I$. Therefore, the combination of $\phi,D$ that ensures the maximal number of points where $f_{r}(I)\cdot S(I)>0$ (the one that we are looking for) delivers the global maximum to $\mathcal{R}_{\Delta f_{r}S}(0)$. This approach will still work correctly even if the mean value $\langle f_{r}\rangle_{I}$ does not lie exactly between the branches in the avoided crossings pattern and intersects one of them since the period $\Pi$ is already found and fixed.

The optimal values for $\phi$ and $D$ are found using the brute force algorithm. It calculates the loss function $\mathcal{L}=-\mathcal{R}_{\Delta f_{r}S}(0)$ on a $50\times 50$ grid of $(\phi,D)$ and takes the minimal value of all. This method is stable and universal due to the evident boundaries on $\phi\in[-\Pi/2,\Pi/2)$ and $D\in[0,1]$. The loss function topography for the avoided crossing patterns is nicely structured and for our example is shown in Fig. 3(c). One peculiarity is that instead of a single minimum it has an area of the same minimal value. Again, this effect comes from the missing zero-padded $f_{r}$ points at some $I$ values. However, any value from this valley suits well enough for our purposes, and the algorithm finds no difficulty in locating it.

Having found values for $\Pi,\ \phi$ and $D$ we now can finally calculate the current of the transmon sweet spot in the case of the avoided crossings pattern:

[TABLE]

Alternatively, for the continuous patterns (which are expected for $f_{ge}(I)<f_{c}$ or $f_{ge}(I)>f_{c}$) we need to apply a different formula:

[TABLE]

Since at this point the algorithm does not yet know which type of the pattern it observes, it may apply some heuristic criterion to guess which equation to apply. When the noise is not too large, such heuristic can be to calculate the maximal absolute differential of the frequency data: $\max_{i>0}|f_{r,i}-f_{{r,i}-1}|$, and compare it to the peak-to-peak amplitude: $\max_{i,j}|f_{r,i}-f_{r,j}|$. For the avoided crossings these values are close and for the smooth dependencies they are not. However, in case of a strong noise this indicator may fail, and the program will have to check both current values to be $I_{ss}$ by fitting the full model two times, and then choose between the two possibilities based on the maximal likelihood.

II.3 Full model fitting

Having performed the aforementioned preliminary steps, it is now possible to fit the full model to the extracted spectral points. To do this, we employ brute force optimization combined with the Nelder-Mead simplex downhill algorithmNelder and Mead (1965) to polish the brute force result. For both methods, we use a common loss function that is based on the maximum likelihood method (again, we presume Gaussian distribution of the points around the model). For the known probe frequency span of the data $\Delta f_{p}$ (Fig. 2(a), whole y-axis) and the set of $N$ extracted points $\{p_{i}\}=\{(I_{i},f_{r,i})\}$, we calculate the loss function as

[TABLE]

where, as in the previous section,

[TABLE]

The condition (5) means that we choose only the model frequencies that lie within the window $\Delta f_{p}$ around the model $f_{c}$ parameter. This ensures, firstly, that in the optimum case we will not take any excess points outside the experimental frequency scan, and, secondly, that we have a singular model value for each current in case of avoided crossings.

To substantiate the choice of the optimization algorithms, we present in Fig. 4 three visualizations of the defined loss function. The plots show how the function behaves if a certain pair of 6 model parameters is varied while others are kept optimal. From the plots it is obvious that the loss function is ill-defined and has a lot of local minima. Moreover, it may not be smooth everywhere because of the condition (5). The $f^{max}_{ge}$ and $d$ parameters present significant difficulty in terms of false minima and low curvature as can be seen from Fig. 4(a). In contrast, $f_{c}$ presents the least difficulty, as can be seen from Fig. 4(b). The last plot Fig. 4(c) shows a very complex structure of the loss function and narrow optimal valleys and serves as an illustration of why the period and phase extraction algorithm is important.

The brute force algorithm acts on the grid specified in Table 2. The ranges in the grid are based on the usual design parameters of the qubit samples in our database, and the number of steps is chosen so that the algorithm reliably finds the optimal valley. After the coarse brute force optimization is done and the optimal valley is located, we apply the Nelder-Mead search on all 6 parameters.

III Results

The above described algorithm was implemented in Python using routines of the SciPyJones et al. (01) library. We present here examples of the detection for all possible qubit-resonator dispositions. The resulting fits are shown in Fig. 5(a)-(c) along with the original data and the error plots characterizing the quality of fit quantitatively; we will further reference them as (a), (b) and (c) respectively.

III.1 Fit quality analysis

In Table 3 we analyse the fitting results for the three cases shown in Fig. 5(a)-(c). Comparing the brute estimation and the final result after the Nelder-Mead search, one can see that the chosen algorithm order works as expected: the exhaustive search finds good initial conditions for the Nelder-Mead, and there is a significant improvement in the loss value of the polished result compared to the brute estimation which is due to the more accurate determination of the cavity frequency $f_{c}$ which leads as well to major shifts in optimal $g$, $f_{ge}^{max}$ and $d$.

However, in turns out that the optimal qubit parameters and the coupling strength may have some errors. For example, in case (b) for $f_{ge}^{max}$ we see a more than 2 GHz difference between the optimal and the more accurate value found for the sample with TTS. For (a) and (c), the differences in $f_{ge}^{max}$ are much smaller (about 10-70 MHz) but instead we see higher errors in $d$ (around 10%). Inaccuracy in the qubit parameters may be qualitatively explained by the low sensitivity of the resonator frequency to the qubit frequency when they are far away (see (1)) and by the strong correlations between the parameters.

To quantify the inaccuracy, we have calculated analytically the Fisher information matrices $\hat{I}$ at the optimal points for the three cases using the SymPyMeurer et al. (2017) package (see Appendix C for details). In Table 3, columns $\min\sigma$ we show the square roots of the Cramér-Rao lower bounds (CRLB) for the variances of our MLE estimators. Let us consider the estimates for $g$ and $f_{ge}^{max}$ parameters. It can be seen, in case (b) $\min\sigma$ are higher than $\min\sigma$ for them in the cases (a) and (c). This is in accordance with the high observed discrepancy between the TTS and estimated $f_{ge}^{max}$. Additionally, using the principal component analysis of the Hessians of the log-likelihood, we find that the most unconstrained component eigenvector contains the parameters $f_{c},\ g,\ f_{ge}^{max},\ d$, and the second one (being an order of magnitude smaller) consists of $\Pi$ and $I_{ss}$. To reduce the correlations, one may choose another parametrization of the model, i.e. replace $d$ and $f_{ge}^{max}$ by $E_{JJ}^{max}$ and $E_{JJ}^{min}$. However, the accuracy of the fit using the presented parametrization is sufficient for practical applications.

III.2 Testing noise robustness

To stress test the algorithm, we artificially deteriorate the data from Fig. 5 by adding complex Gaussian noise of various variances to the complex data at each ($f_{p}$, $I$) before fitting. However, since the original pictures already contain an experimental noise, we first have to define a common scale for all three cases (a), (b) and (c) to be able to compare them. The signal-to-noise ratio (SNR) turns out to be a sensible universal measure of noise, and is defined for our data in a similar manner to the SNR in the resonator fitting testsProbst et al. (2015). Since the data for a resonator lie on a circle in the complex plane, it is convenient to take its radius $r$ as the signal amplitude. Then, if the total noise is of the form $(\xi_{1}+i\xi_{2})/\sqrt{2}$, where $\xi_{1},\ \xi_{2}$ are distributed normally with zero mean and standard deviation (SD) $\sigma$, the SNR is defined as

[TABLE]

where $\sigma_{0}$ is the SD of the original noise and $\sigma_{1}$ is the SD of the added noise. Obviously, the resulting SNR in this scale can’t be larger then the original SNR=$r/\sigma_{0}$.

Using the described common scale for the measurements (a), (b) and (c), we test the algorithm stability by launching it 50 times for different added noise realizations with various standard deviations $\sigma_{1}$. The added noise power is steadily increased, and the resulting SNR is being calculated according to Eq. (6). We then record all fitting results, and calculate their mean value, 25th and 75th percentiles for each $\sigma_{1}$. Since the original noise is fixed, the variance of the fitting results will go to zero at SNR=$r/\sigma_{0}$ where the added noise is negligible and the fit converges to the same values as in Table 3.

In Fig. 6 we plot the results of the test outlined above. As one can see from the graph, the algorithm becomes stable above SNR=3, and the lowest SNRs before divergence occurs are around 2. The original data SNR for (a) is 19, and thus in the tests $\sigma_{1}>5\sigma_{0}$. Conversely, the original data SNR for (b) and (c) are much lower, and thus $\sigma_{1}\approx\sigma_{0}$ in the tests. This means that for (b) and (c) the variance of the fitting results shown with the clouds in Fig. 6 is lower than it would be if all the noise was generated.

In overall, the detection algorithm is stable even at very low SNRs where the patterns are barely visible by the human eye. The stability is limited mostly by the resonator frequency detection procedure.

III.3 Performance

We have tested the performance of the algorithm on a 5-year old 2-core Intel i5-3337U CPU @ 2.5 GHz. The results are presented in Table 4 for the case (a). As can be seen, the procedure takes around 7 seconds to complete. On a contemporary Intel Core i7-7700 CPU, it takes about two seconds.

The heaviest computational tasks are the frequency extraction part with the $circlefit$ library (2.7 seconds, most time spent in iterative calibration fits) and the brute search and Nelder-Mead optimization part where the time costs mostly come from the square root calculation necessary for Eqs. (9) and (13). However, this is still fast enough, firstly, because recording a resonator spectrum such as one in the example takes around 20 seconds, and, secondly, it is much faster than doing it manually.

Due to the complexity of calculating the value of the loss function, the frequency extraction scales at least as $\mathcal{O}(NM)$ if $N$ and $M$ are the number of points for current $I$ and probe frequency $f_{p}$, respectively, and the following optimization scales linearly with $N$.

IV Summary and discussion

In summary, we have developed and implemented a method for automatic processing of single-tone spectroscopic data that uses maximum likelihood estimation to extract the physical parameters of cQED systems from measurement results. We have shown that it is possible to estimate all relevant circuit parameters by tracing only the behaviour of the resonance dip (or peak) for any kind of qubit-resonator disposition. Our method is very robust with respect to noise down to SNR=3 and is fast enough to be used in practice. The accuracy of the algorithm is not sufficient for starting qubit gate calibrations right after the single-tone spectroscopy step, but it facilitates much faster two-tone spectroscopy with significantly reduced scan ranges. Additionally, our automated evaluation method allows to distinguish between functional and defective samples.

The software implementation of the algorithm in Python can be found on GitHub222https://github.com/vdrhtc/Measurement-automation/blob/master/lib2/fulaut/AnticrossingOracle.py. The code as well accepts data as a Python dictionary with corresponding keys. Currently, only transmon qubits are supported; however, implementing this analysis for another types of qubits seems rather straightforward.

Considering the downsides of the suggested approach, we note that the implementation of the algorithm is not perfectly optimized. The brute force search is not fully parallelized and might be replaced with a more sophisticated global optimization algorithm such as simulated annealing or basin hopping. Another problem that may slow down the execution is that for large qubit-cavity detuning and high noise two very different sets of parameters may minimize the loss function equally well. An illustration for this statement can be seen in Fig. 7(a),(b) where the algorithm was launched twice using different grid ranges for $f_{q}^{max}$. In one case (orange) the search was done below, and in the other (green) above the mean resonance frequency $\langle f_{c}\rangle_{I}$. The resulting fits are almost equally accurate, and without additional information about the system it is not possible to reliably decide between the two. Therefore, some hints should be passed to the algorithm to avoid ambiguity; for instance, a prior probability density may be imposed on the parameter values and the MLE replaced with maximum a posteriori method. Otherwise, it would be necessary to check both possibilities with other methods, i.e. via two-tone spectroscopy. Since for the transmon qubits the case in Fig. 7 is the only type of ambiguity that may occur, we introduce an optional flag for the program that specifies where to look for the qubit transition frequency during the brute force search, above or below the frequency of the resonator.

In the future we are planning to use the proposed algorithm along with automated analysis of two-tone spectroscopic measurements to build a full system that will automatically calibrate the qubit samples from scratch.

V Acknowledgements

We gratefully acknowledge valuable discussions with I. Besedin, V. Ryazanov, A. Dmitriev and P. Fedorov. This work was partially supported by the Ministry of Education and Science of Russian Federation in the framework of Increase Competitiveness Program of the NUST MISIS (contract no. K2-2017-081). Experiments and analytical tools were developed with the financial support from the Russian Science Foundation (contract no. 16-12-00095). Numerical simulations were performed with the support of Russian Foundation for Basic Research, grant No. 19-12-80006.

Appendix A Transmon Hamiltonian

The simplest version of this qubit consists of a Josephson junction shunted with a large capacitor. Flux tunability of the frequency is attained by replacing a single Josephson junction with a SQUID as in Fig. 8(a) and applying external magnetic flux $\Phi_{e}$ to its loop. This scheme can be equivalently represented by a capacitively shunted junction of tunable Josephson energy, Fig. 8(b). The Hamiltonian for such system is as follows:

[TABLE]

where $E_{C}=e^{2}/2C_{\Sigma}$, $C_{\Sigma}=C_{s}+C_{1}+C_{2}$, is the charging energy, $\hat{n}$ and $\hat{\varphi}$ are the operators for the Cooper pair number and the phase of the qubit island. For the equivalent Josephson energy $E_{J}(\Phi_{e})$ one obtains

[TABLE]

where $E_{J\Sigma}=E_{J1}+E_{J2}$ ($E_{J1},\ E_{J2}$ are the single junction energies), $d=\frac{E_{J1}-E_{J2}}{E_{J1}+E_{J2}}$ is the asymmetry of the SQUID. As one can notice, the dependence is periodic in $\Phi_{e}$.

It is also possible to derive analytical expressions for the energy levels and transition frequencies for this type of qubits. The energy of the $m$th level is (in the limit of $E_{J}\gg E_{C}$) Koch et al. (2007)

[TABLE]

From this equation, the qubit transition frequency ($\ket{0}\rightarrow\ket{1}$, or $\ket{g}\rightarrow\ket{e}$, or the $ge$ transition) may be approximated as

[TABLE]

where $f_{ge}^{max}=\sqrt{8E_{J}(0)E_{C}}-E_{C}$ is the maximal possible frequency. The error of this approximation is

[TABLE]

which is non-negligible when $d$ is small and $\Phi_{e}\approx\Phi_{0}/2$. This is satisfactory since the expression (9) itself is no longer valid in these conditions because $E_{J}\approx E_{C}$ there. On the contrary, this approximation simplifies the model since now it depends just on two parameters ($f_{ge}^{max},d$) instead of three. Therefore, we choose (10) as the theoretical curve in fitting.

One final note is that in the real-life applications is not possible to know directly the flux $\Phi_{e}$ that is threaded through the SQUID. The experimenter usually knows only the current $I$ (or voltage) which he applies to some coil that is connected inductively to the SQUID. Then $\Phi_{e}=MI+\Phi_{r}$, where $M$ stands for the mutual inductance of the coil and the SQUID, and $\Phi_{r}$ is some residual flux inherent to the sample. So, in the main text we only use the current $I$, the period in current $\Pi=\Phi_{0}/M$ and the sweet spot current $I_{ss}=-\Phi_{r}/M$ and not the corresponding magnetic fluxes.

Appendix B Circuit QED

The readout of the superconducting qubits is now predominantly done using an ancilla system which is usually implemented as a superconducting microwave resonator which acts as an electromagnetic cavity in the standard cavity quantum electrodynamics (QED). Truncating the qubit to two levels ($ge$ transition), one may obtain the following Hamiltonian for the compound cavity-qubit system (in RWA):

[TABLE]

where $f_{ge}$ is the qubit frequency, $f_{c}$ is the cavity frequency and $g$ is the coupling strength.

As long as the RWA is done, this Hamiltonian may be diagonalized analyticallyBlais et al. (2004):

[TABLE]

By substituting the dependence of the qubit frequency $f_{ge}\equiv f_{ge}(\Phi_{e})$ into these equations, we can get straightforwardly the full system spectrum in dependence on the magnetic flux. We note here that, in general, $g$ may as well not be constant since the relevant matrix elements of the qubit may depend on the control parameters just as its frequency. For the transmonKoch et al. (2007)

[TABLE]

There are two useful observations from the expressions (13) and (B). The first one is that for the avoided crossings the exact value of $g$ is only relevant when $(f_{ge}-f_{c})^{2}/4\lesssim g^{2}$. Therefore, if its dependence on the control parameters is slow enough (like for the transmon case), a constant $g$ fixed to its value $g(\Phi_{e}^{*}),\ \Phi_{e}^{*}:f_{ge}(\Phi_{e}^{*})=f_{c}$ is a good approximation in the case (a) of Fig. 9. For the cases (b) and (c), in contrast, the constant $g$ model is, in general, not valid. Fortunately, for the transmon $g^{2}(\Phi_{e})=g_{max}^{2}\sqrt{k(d,\Phi_{e})}\propto f_{ge}(\Phi_{e})$ and thus equations (13) may be rewritten in terms of $f_{c}^{\prime}=f_{c}-g_{max}^{2}/f_{ge}^{max}$ and $g^{\prime}=g_{max}\cdot\sqrt{f_{c}/f_{ge}^{max}}$ (note that now $g^{\prime}$ is flux-independent). These approximations are accurate for Fig. 9 (b), (c) when $g\ll f_{c}$. Our algorithm thus finds $f_{c}^{\prime}$ and $g^{\prime}$ while $f_{ge}^{max}$ and $d$ are not affected. For other types of qubits, in the linear coupling regime $g$ should be proportional to the matrix element of the qubit-resonator coupling operator. Therefore, the model for the algorithm needs to include $g(\Phi_{e})=g^{\prime}\cdot\bra{g}\hat{V}\ket{e}(\Phi_{e})$ explicitly where $g^{\prime}$ is the fitting parameter and $\hat{V}$ is the dimensionless qubit-resonator coupling operator term in the Hamiltonian.

In Fig. 9 we have used the equations (10), (12) and (13) to model a tunable transmon interacting with a cavity for various $\Phi_{e}$ and various $f_{ge}^{max},\ d$. In the lower row of the figure, one can see that it is possible to extract the dependence of the modified cavity frequency $f_{c}^{\prime}$ on $\Phi_{e}$ (in the main text, we call it $f_{r}$ as the experimentally observed resonator frequency); for example, the avoided crossings pattern can be directly observed in Fig. 9(a), and the other two possible behaviours for the qubit entirely above or below the resonator in Fig. 9(b),(c). To shorten the notation, in the following we will define the corresponding branch frequencies of (13) as $f_{\pm}=(E_{\pm,0}-E_{g,0})/2\pi$.

Appendix C Maximum likelihood, Fisher information

For a vector $\mathbf{w}$ of parameters, we can calculate the likelihood $\mathcal{L}(\mathbf{p}|\mathbf{w})$ to observe $N$ sample data points $\mathbf{p}=\{I_{i},f_{r,i}\}_{i=1}^{N}$ if they are presumed to be normally distributed around the model ($\mathbb{N}_{\mu_{i},\sigma}$, $\mu_{i}=\mathcal{M}(I_{i},\mathbf{w}),\ \sigma=\text{const}$) and uncorrelated. The sought value is the product of the probability densities for observing the individual data points:

[TABLE]

The maximum likelihood method consists of finding the optimal parameter vector $\mathbf{w}^{*}$ such that $\mathcal{L}(\mathbf{p}|\mathbf{w}^{*})$ is the highest among all allowed $\mathbf{w}$. The vector $\mathbf{w}^{*}$ is called an estimator for the true parameter vector $\mathbf{w}^{0}$ and is itself a random value. If we take the negative logarithm of the likelihood, we will get the negative log-likelihood function:

[TABLE]

Minimizing this function is equivalent to maximizing the likelihood since the logarithm is monotonic. As the second term in the right part does not depend on $\mathbf{w}$, the task is reduced to minimizing the sum of squares, or the $\chi^{2}$; this argument motivates the choice of the form for the loss function (4).

Next, we consider the variance of the MLE estimator $\mathbf{w}^{*}$ that we get from the fit. Neglecting its bias (in general, for nonlinear MLE it is always present but can only be calculated approximatelyCox and Snell (1968) and scales as $1/N$), we can use the multivariate Cramér-Rao inequality to estimate the lower bound of its variance. In the general case, the inequality is as followsKeener (2011); Schervish (2012):

[TABLE]

Here $\text{Var}_{\mathbf{w}}$ means the variance of the estimator when the true model distribution parameters are $\mathbf{w}$, and $\hat{I}(\mathbf{w})$ is the Fisher information matrix. We show below that this matrix depends on $\mathbf{w}$, so to get the relevant bound we need to calculate it at $\mathbf{w}=\mathbf{w}^{0}$. Since $\mathbf{w}^{0}$ is not known for the real data, one will have to use $\mathbf{w}^{*}$ as an approximation for $\mathbf{w}^{0}$. Using the definition of $\hat{I}(\mathbf{w})$, the properties of $\ln\mathcal{L}$ and assuming the possibility of double differentiation, we can writeKeener (2011); Schervish (2012):

[TABLE]

where $E_{\mathbf{w}}$ means the expectation under $\mathcal{L}(\mathbf{p}|\mathbf{w})$, i.e. the average over possible realizations of $\mathbf{p}$ distributed according to (14), and

[TABLE]

is the Hessian operator. Next, using (14) and the fact that $E_{\mathbf{w}}[f_{r,i}]=\mu_{i}(I_{i},\mathbf{w})$, we can derive the analytic expression for $\hat{I}(\mathbf{w}$):

[TABLE]

Notably, this formula is valid for any type of nonlinear MLE with Gaussian noise. In overall, to obtain the Fisher information matrix, one needs, firstly, to calculate the Hessian of the model function $\mathcal{M}(I,\mathbf{w})$ at $\mathbf{w}=\mathbf{w}^{*}$ (or $\mathbf{w}^{0}$ if the data is synthetic). Secondly, the measurement error variance $\sigma^{2}$ is required. A common method to estimate it is to estimate the variance based on the sample $\{f_{r,i}-\mu_{i}\}$. However, for the nonlinear models it is not possibleYe (1998); Andrae, Schulze-Hartung, and Melchior (2010) to derive a valid analytic expression for the corresponding estimator $(\sigma^{2})^{*}$ due to the unknown number of the effective degrees of freedom. One way to solve this problem is to find $\sigma^{2}$ independently by measuring the variance of $f_{r}$ at some fixed $I$. Another way is to use the expression valid for the linear regression $(\sigma^{2})^{*}=\chi^{2}/(N-M)$ where $M=\dim\mathbf{w}$ and then test its validity using the Monte-Carlo approach. Our numerical tests with 5000 runs show that even when having $M=6$ and $N=18$, $(\sigma^{2})^{*}$ agrees well with the known $\sigma^{2}$ of the synthetic data: $|E[(\sigma^{2})^{*}]-\sigma^{2}|<0.01\sigma^{2}$ and therefore can be used in error calculation.