Improving Transmon Qudit Measurement on IBM Quantum Hardware

The Hilbert space of a physical qubit typically features more than two energy levels. Using states outside the qubit subspace can provide advantages in quantum computation. To benefit from these advantages, individual states of the $d$-dimensional qudit Hilbert space have to be discriminated during readout. We propose and analyze two measurement strategies that improve the distinguishability of transmon qudit states. Based on a model describing the readout of a transmon qudit coupled to a resonator, we identify the regime in hardware parameter space where each strategy is optimal. We discuss these strategies in the context of a practical implementation of the default measurement of a ququart on IBM Quantum hardware whose states are prepared by employing higher-order $X$ gates that make use of two-photon transitions.


I. INTRODUCTION
Conventional quantum computing is based on qubits which are realized on two-level subspaces of a larger physical Hilbert space.A number of physical realizations of qubits have been proposed and implemented on various platforms.These include superconducting qubits [1], trapped ions [2], cold atoms and Rydberg atoms [3], as well as electron spins in quantum dots [4].On all of these platforms, it is necessary to isolate the twodimensional qubit space from the remaining states of the physical Hilbert space to avoid leakage out of the computation space.However, utilizing qudits, i.e., ddimensional building blocks of quantum computation, can provide advantages [5][6][7][8][9][10][11][12][13][14][15][16][17].Examples include implementing an ancilla qubit within the second and third excited states of a qudit [6].Another example is the socalled shelving [9]: by transferring the population of the first excited state to the second excited state, the error of identifying the ground state decreases.
Superconducting qubits [18,19] are prominent building blocks of noisy intermediate-scale quantum systems.The most promising example is the so-called transmon that can effectively be described as a quantum anharmonic electromagnetic oscillator.In this system, the two lowest-energy levels are identified as the qubit.Taking into account higher-lying transmon levels leads to a natural realization of a superconducting qudit.The smallest extension of the qubit is the qutrit, i.e., a three-level system.Qutrits have been used to implement a Toffoli gate [10] with a significantly lower number of elementary gates compared with a realization based on two-level systems.Another interesting example is the recent experimental demonstration of a qutrit Greenberger-Horne-Zeilinger (GHZ) state [11].
In general, if one is interested in measuring a qudit state, a proper classification of all levels involved is needed.In [20], the qubit state is determined by a fit of the time evolution of the system.In setups which * tobias.kehrer@unibas.chdo not provide time-resolved data, such as the current IBM Quantum [21] devices, other methods of separating qudit states have to be employed [22][23][24].The strategies described in both [22] and [23] involve exciting the qudit-resonator system at readout drive frequencies other than the default frequency.At the default frequency, the distinguishability of the ground state and first excited state is maximized, whereas using the adapted frequencies aims at optimizing distances between different pairs of qudit states.
In this paper, we propose and evaluate improvements of the measurement scheme of transmon qudit states by enhancing their distinguishability.To optimize the readout, we determine the measurement errors from the assignment matrix whose entries denote the probability to classify a measurement outcome to a state |i⟩ even if state |j⟩ was prepared.This assignment matrix is calculated using qudit-state dependent resonator steady-state amplitudes obtained from a model describing the readout of a transmon qudit by driving a coupled resonator.The default measurement schedule of most superconducting quantum hardware consists of a single-tone drive applied to the readout resonator.The frequency of the tone is chosen to maximally separate the ground and first excited states.The strategies we propose are based on modified readout resonator drive frequencies that take into account the separation of all qudit states.These strategies include a single-frequency as well as a multifrequency readout scheme.For a ququart, viz., the four lowest states of a qudit, we compare the proposed strategies in simulation and show that depending on hardware parameters, both strategies can be beneficial.We furthermore compare the model to a measurement of the drive-frequency-dependent resonator states on a current IBM Quantum device.
The paper is organized as follows.In Section II, we present a mean-field model describing the readout sequence of a transmon qudit coupled to a harmonic readout resonator.Based on this model, we calculate the readout drive-frequency-dependent assignment errors between multiple states that in some limits can be expressed analytically.In Section III, we analyze both proposed readout schemes that aim to minimize these errors.We arXiv:2307.13504v2[quant-ph] 19 Jan 2024 compare the data that we generated on current IBM Quantum hardware, see Section IV, to the readout model and strategies discussed in Sections II and III.To improve the state preparation required in this procedure, we propose to add two-photon transitions to the universal gate set of qudit gates [7] and show that this will reduce the execution time of certain qudit circuits and the duration of X-gate calibrations, see Appendices B 2 and C. Finally, we conclude in Section V.

A. Effective Hamiltonian
The fundamental building blocks of a superconducting quantum computer are a quantum anharmonic oscillator, i.e., the transmon qudit, coupled to a harmonic oscillator, i.e., the readout resonator.Following [19], the corresponding effective Hamiltonian obtained by treating the Jaynes-Cummings interaction between the qudit and the resonator as a small perturbation, see Appendix B 1, reads The parameter ω j is the energy (see Appendix A) of the bare qudit state |j⟩, ω r is the energy of the readout resonator, and a ( †) is its annihilation (creation) operator.
Here and in the rest of the paper, we set ℏ = 1.The second-order corrections to the qudit and resonator energies, χ j−1,j and χ j , are defined in Eqs.(B11) and (B12).Additionally, we describe a coherent driving of the resonator at frequency ω d by [1] which enables the readout of qudit states.

B. Readout of States
The readout of a transmon qudit, in short, consists of driving the readout resonator while recording the response signal.We model the time evolution of a general quantum state ρ comprised of a qudit and its readout resonator by the following Lindblad master equation: where κ is the decay rate of the resonator.Using the effective Hamiltonian given by Eq. ( 1) and assuming the qudit to be in state |j⟩, we arrive at the equation of motion of the mean-field amplitude The fact that A depends on the qudit state |j⟩ is used to discriminate different qudit states.If the qudit is in a mixture or superposition of states, this measurement procedure projects the qudit onto one of its Fock states [1].The general form of the complex value returned by an IBM Quantum device is where k(t) encodes the kernel integration instructions, see meas kernel [25] in Qiskit [26], and T is the total duration of the measurement.The choice k(t) = exp(iω d t) corresponds to integrating the measurement signal in the rotating frame of the drive, see Section II B 1, whereas the choice k(t) = exp(iω m t) corresponds to a frame rotating at an arbitrary modulation frequency ω m , see Section II B 2.
In this paper, we mainly consider the offset charge configuration n g = 0.The value of n g influences the transmon qudit energy spectrum, see Appendix A. Note that due to the significant dependence on n g of the third and higher excited states, their corresponding readout resonator states may be smeared out in phase space if charge noise is present.In Fig. 1 we present an overview of the energy dispersion ϵ 3 of the third excited state defined in Eq. (A5) for several IBM Quantum devices.Since ϵ 3 decreases with increasing E J /E C , qudits that lie in the upper-right region are preferred in general.

Rotating Frame of Drive
In this section, we will work in the rotating frame of the drive.Quantities in this frame will be denoted by the superscript (d).Since Eq. ( 5) is defined in the laboratory frame, we choose k(t) = exp(iω d t) to transform the signal into the rotating frame of the drive and obtain Here, A of the readout resonator move on a circle centered at A c = −ie iϕ Ω/2κ with diameter Ω/κ.At resonance ω (j) d,0 = ω r + χ j , the states reach the maximum amplitude 2A c .For qudit readout, it is important that the distance j | between two qudit-statedependent resonator states is large.We can identify two regimes of how the positions of the states in phase space depend on the readout drive frequency.For a large resonator decay rate κ > |χ i − χ j |, all states are close to the position of maximum amplitude within the same frequency range.In this case, d i,j exhibits only one maximum at In contrast, for a small resonator decay rate κ < |χ i −χ j |, the frequency ranges where the state amplitudes A (d) j are close to the maximum amplitude do not match.Here, two drive frequencies where d c denotes the diameter of the circle on which the states move.Thus, at ω (i,j) d,± , both states are located on opposite sides of the circle, which is the maximum separation they can obtain.
If we set the drive frequency to ω (i,j) d,0 (or ω (i,j) d,± ), i.e., maximizing the distance between state |i⟩ and |j⟩, the distance between other pairs of states is in general reduced and hence not optimal for discrimination of these states.Therefore, in Section III, we present two measurement strategies to mitigate this issue.

General Rotating Frame
In a frame of a general rotation frequency ω m , i.e., choosing k(t) = exp(iω m t), the state reached in the longtime limit κT ≫ 1 is time dependent, Ā(m) where sinc(x) = sin(x)/x and the superscript (m) is used to denote quantities in this frame.The difference between Eqs. ( 7) and ( 12) is an additional factor of sinc peaking at ω d = ω m .These resonator state amplitudes and their dependence on the drive frequency ω d are also visualized in Fig. 2a.The states move on the black circle with diameter Ω/κ, whereas the motion of the states A (m) j follows a distorted circle (colored lines).

III. MEASUREMENT STRATEGIES
In the previous section, we presented a model describing the readout on superconducting quantum hardware.The centers of the Wigner functions of the coherent readout resonator states when the qudit is in state |j⟩ are given by A j .Due to intrinsic quantum noise and hardware limitations, the possible readout resonator states for each qudit state overlap.This leads to potential misclassification and thus measurement errors when reading out the qudit states.
In the following, we propose two strategies for improving qudit readout compared to the default measurement scheme that utilizes a single resonator drive frequency that optimizes the classification of |0⟩ and |1⟩.The first strategy consists of finding a single frequency that maximizes the distinguishability between all d qudit Fock states.In the second strategy, we allow for multiple different drive frequencies.

A. Assignment Matrix
To arrive at a measure of the distinguishability of states, we introduce the measurement assignment matrix M [27].The qudit-state-dependent resonator states are defined by their steady-state amplitude A j .We assume their Wigner functions to follow a two-dimensional Gaussian distribution, centered at A j with standard deviation σ j larger than the intrinsic quantum noise.The elements of M are given by  17) and (18).Following Appendices A and B 1, for these plots, we determine EJ /EC by the qubit parameters ω0,1 and α1 from ibm lagos Q4 (July 7, 2023).Moreover, we choose g/2π = 100 MHz, Ω/2π = 100 MHz, κ/2π = 5 MHz, T = 0.35 µs, σj = 0.13Ω/κ, ϕ = 0, and ng = 0. and define the probability to classify a measurement as state |i⟩ even if state |j⟩ was prepared.The region corresponding to each state |i⟩ is defined by the maximum likelihood estimator (MLE) leading to where Θ denotes the Heaviside function.For σ j = σ (valid assumption for this hardware setup, see the discussion about the distribution of σ of Gaussian fits in the second paragraph of Section IV), the MLE is equivalent to the minimum distance estimator (MDE) that implies Using the MDE, a data point z is assigned to the region of state A i if its Euclidean distance to all of the other states A k is larger.In contrast, using the MLE, a data point z is assigned to the region of the state A i that has the largest probability density.For simplicity and since in our measurements all σ j are comparable, we choose the MDE throughout this paper.Ideally, M i,j = δ i,j , meaning perfect measurement.We define two measures ξ j and ξ, where ξ j is the probability of misclassifying the qudit state |j⟩.The dependence of ξ j and ξ on the readout resonator drive frequency is shown in Figs.2b and c.The measurement errors ξ j achieve their minima at different readout resonator drive frequencies.If ω m ̸ = ω d , the locations of the minima cannot be distinguished as well as for ω m = ω d .In the current setup of IBM Quantum hardware, the frequency ω m of the rotating frame cannot be changed.Therefore, the difference between the frequency dependencies of all ξ j is less pronounced.Note that for setups where all qudit states lie on a circle (e.g., always for qutrits) and σ j = σ, M i,j can be expressed in terms of Owen's T function, see Appendix E. This allows for fast numerical calculation of Eq. ( 14).

B. Finite Sampling
In experiments, measuring an unknown state |ψ⟩ = j c j |j⟩ in the Z basis is equivalent to estimating its populations p j ≡ |c j | 2 based on a set of N data points {z j }, also called shots.For each shot, the total state is projected onto one of the d qudit states |j⟩ with probability p j .Therefore, the total probability distribution of measuring one shot at z given ⃗ p is a sum of all d Gaussians defined in Eq. ( 13) weighted by p j .The measurement task can be understood as learning the parameters ⃗ p = (p j ) of this multimodal probability distribution, where A j and σ j are obtained from a separate measurement.The space of possible ⃗ p, can be mapped to a (d − 1)-simplex, using the normalization condition of ⃗ p.
We describe the measurement analysis as follows.Each shot is labeled by a state |j⟩ depending on its phase-space distance (MDE) to the d qudit Fock states.The resulting list of counts ⃗ N can be used to obtain information about ⃗ p.In Appendix D, we show that by using Bayes' formalism, the probability distribution of ⃗ p given ⃗ N is the Dirichlet distribution, The component N j of ⃗ N equals the number of shots classified as state |j⟩.The assignment matrix M reflects the fact that some shots are classified incorrectly.Given ⃗ N , the location of the maximum (also called mode) can be computed analytically, This result is similar to a common procedure known in Qiskit as "measurement error mitigation".Note that applying the inverse of M to ⃗ N can lead to negative components of ⃗ p mode .In Qiskit, this problem is circumvented by approximating ⃗ p mode by the valid ⃗ p ′ that is closest (in two-norm) to ⃗ N , see method least squares in qiskit.utils.mitigation.filters.py[26], Equation ( 23) is the estimate of the state populations ⃗ p after measuring ⃗ N shots.We will use the uncertainty of these estimates, viz., the numerically calculated standard deviations SD[p j ], to decide which of the proposed strategies performs best, i.e., exhibits the smallest standard deviation.

C. Comparison of Strategies
We consider two strategies that make use of either one or multiple drive frequencies.In the default readout scheme of superconducting quantum hardware, measurement pulses with a single drive frequency that maximizes the distinguishability between the qubit states |0⟩ and |1⟩ are applied.
The first strategy we propose replaces the default frequency by the one that optimally separates all qudit states in phase space simultaneously.Since, in general, the state |ψ⟩ that we want to measure is unknown, we suggest to optimize ξ, see Eq. (18), which is the average of the individual measurement errors ξ j .
The second strategy uses N/d shots for each of the d different frequencies at which individual states are most isolated, i.e., ξ j are minimal.We will show that this strategy is advantageous in cases when there is no single frequency at which all states are separated well enough.Hardware parameters and the state to be measured determine which of the two strategies outperforms the other.To compare both strategies, we draw N = 1000 samples from the probability distribution given by Eq. ( 19) for σ j = σ and an equal-superposition state p j = 1/d.The drive frequencies we use are the locations of the minimum of ξ for the single-drive strategy and the minima of ξ j for the multifrequency strategy.Each sample is classified using the MDE, see Eq. ( 16), i.e., by its Euclidean distance to the nearest state |j⟩.The final probability distribution for the p j of the single-frequency strategy is given in Eq. (21).The final probability distribution for the multifrequency strategy is the normalized product of the term in Eq. ( 21) for each measurement frequency ω k , where ⃗ n k is the list of counts of classified shots for the k th measurement frequency.The standard deviation SD[p j ] is computed numerically from this distribution.Figure 3 shows the dependence of the ratio of both averaged standard deviations, of p j on hardware parameters σ j = σ and κ.The blue region corresponds to setups for which the standard deviation SD s of p j using a single-drive frequency scheme is smaller.In contrast, the red region corresponds to hardware configurations where it is beneficial to measure at multiple frequencies, i.e., SD m of the multidrive frequency scheme is smaller.The gray region indicates Here, the elements Mi,j of the assignment matrix equal the relative number of shots, Ni/N , that are classified as |i⟩ even if |j⟩ is prepared.The horizontal gray line denotes the measurement error ξ def obtained using the default measurement pulse.(c) Measurement errors ξj based on Gaussian fits to the data presented in (a) and the assignment matrix M defined in Eq. ( 14).Using the centers and average σ of Gaussian fits for each readout resonator drive frequency, we calculate M numerically.The ξj shown in (b) are larger than those in (c) since they do not only represent assignment errors, but also include additional errors such as qudit decay, leakage, and imperfect state preparation.In both (b) and (c), the minimum of the average assignment error ξ will be smaller than ξ def obtained by the default pulse if the modulation frequency is chosen such that ωm = ω d , cf.Figs.2b and c.
parameter values for which both standard deviations exceed SD s/m ≥ 0.1.Since the expected values p j lie in [0, 1], this threshold corresponds to an uncertainty of at least 10%.
The overall trend is that for small σ, i.e., strongly located Gaussians, the single-frequency strategy performs at a similar, slightly better level than the multifrequency strategy.The multifrequency strategy is preferable for large σ, when the overlap of the Gaussians would be too large using a single drive frequency.Intuitively, this is expected since, for small κ and large σ, only one state is isolated from the others which group together at the origin in phase space, see discussion of regimes κ ≶ |χ j − χ j+1 | in Section II B 1. We also added lines of constant relative uncertainty σκ/Ω.Along these lines, the Gaussian widths σ are fixed in units of the diameter Ω/κ of the circle on which the states move in the rotating frame of the drive.The solid line corresponds to σ = 0.13Ω/κ chosen in Figs.2b and c, whereas the dashed line approximately matches the threshold of SD s/m ≥ 0.1.Following the solid black line, the standard deviation of the singlefrequency strategy appears to exhibit a minimum around κ/2π = 1 − 2 MHz.For fixed σκ/Ω and small resonator decay rates κ, the states move around the circle rather individually, whereas for large κ, the states move as a group, cf.Section II B 1.

IV. MEASUREMENT OF A QUQUART
In this section, we will compare the model described in Section II to data obtained from ibm lagos Q4 on July 7, 2023.We prepare the four lowest Fock states of the transmon qudit and measure them for various readout drive frequencies.In Appendix C, the preparation of the individual qudit Fock states is described.Sequences of X gates defined as simple Gaussian pulses are used to prepare the four lowest Fock states of the transmon qudit.As mentioned earlier, we make use of higher-order X gates, see Appendix B 2. These provide a reduction of the execution time of certain quantum circuits and a reduction of the duration of X-gate calibrations.
In Fig. 4a we show the measurements of the four lowest Fock states for various readout resonator drive frequencies.This plot is the experimental equivalent of Fig. 2a.For each Fock state and for each readout resonator drive frequency, we take N = 2000 shots while keeping the other drive parameters fixed at the default values.For ω 0 = 7.2463 GHz, Fig. 4a shows all shots in the color of the prepared Fock state.This value of ω 0 is −5.5 MHz off the default frequency reported by the IBM Quantum device.Black crosses highlight the centers of the Gaussian fits.For other drive frequencies, we only plot the centers of the Gaussian fits as colored lines (with a white shadow).The straight black lines denote the boundaries of regions (defined via MDE, see Eq. ( 16)) that are as-signed to one Fock state.
We analyze the measurement errors in two ways.First, we define the elements M i,j by the relative number of shots, N i /N , classified as |i⟩ even if |j⟩ is prepared.In this way, M incorporates misclassification errors but also additional errors such as imperfect qudit state preparation.From this matrix, we obtain the errors ξ j , displayed in Fig. 4b.Second, we use the centers of the Gaussian fits for each qudit state and for each value of the resonator drive frequency and a fixed value of σ to compute the assignment matrix defined in Eq. (14).By examining these Gaussian fits, we find a narrow distribution of the σ values: σ = (0.302 ± 0.017) (same arbitrary units as in Fig. 4a).The resulting errors ξ j are shown in Fig. 4c.Here, the ξ j only represent errors that arise from misassignment of shots drawn from the multi-Gaussian distribution, see Eqs. ( 13) and (14).Since real devices feature other sources of error, e.g., qubit decay, leakage, and imperfect state preparation, the values of ξ j presented in Fig. 4b are larger than in Fig. 4c.
Our model, visualized by the theory plots in Figs.2a  and c, shows qualitative agreement with the data presented in Figs.4a and c.In both Figs.4b and c, the horizontal gray line ξ def denotes the average assignment error of the four lowest Fock states using the default readout pulse and should be compared with the solid black line ξ.The corresponding data were taken from Rabi calibration measurements, similar to Fig. 6b, at the drive amplitude that is closest to the fitted optimum.
We find a dependence of the measurement errors ξ j on the readout resonator frequency as expected.The data presented in Figs.4b and c suggest that the default measurement frequency is not ideal to separate all four qudit states.However, the minima appear at only slightly different positions.Note that the difference in positions is only small due to IBM Quantum software/hardware limitations: ω m cannot be set to its ideal value ω m = ω d , see Section II B 2. We expect the impact of varying the readout resonator drive frequency to be much higher if it is possible to analyze all data in the rotating frame of the drive, compare Figs.2b and c.
In this paper, we focused on the analysis of only four qudit states since the readout of higher excited states beyond |3⟩ becomes difficult for several reasons.Higher excited states are more sensitive to charge noise, see Fig. 5. Since χ j depends on the qudit spectrum ω j , charge noise leads to ambiguous steady-state amplitudes.In addition, finding single drive frequency that properly separates all qudit states becomes more difficult with an increasing number of qudit states.For example, for the IBM Quantum device that we utilized in this paper, we estimate χ 1 < χ 4 < χ 2 which indicates that the steady-state amplitude corresponding to |4⟩ lies between |1⟩ and |2⟩.We expect that the more states are involved, the better the performance of a multifrequency strategy in comparison to a single-frequency strategy given a small κ < |χ i − χ j |, cf.Fig. 3.

V. CONCLUSION
We have presented a model that describes phase-space measurement data of qudit states on superconducting quantum hardware.Our model qualitatively matches the data that we generated on a current IBM Quantum device.For qudit-state preparation, we employ higherorder X gates between |j⟩ and |j + 2⟩.This scheme leads to a reduction of the execution time of qudit quantum circuits as well as of the duration of X-gate calibrations.Based on our model, we have compared the performance of two measurement strategies, a single-frequency and a multifrequency scheme, in simulations.For each strategy, we have identified the regime in hardware parameter space where it is optimal.The multifrequency strategy is superior when the qudit-state-dependent resonator states overlap significantly.
To use the full potential of both strategies, it is necessary to adjust the modulation frequency ω m of the device.This is currently not possible on IBM Quantum hardware.Despite these software/hardware restrictions, we still find differences in the frequency locations of the minima of the individual measurement errors ξ j and an improvement over the measurement error ξ def using the default measurement pulse.We expect a better performance of the strategies for setups that operate in the rotating frame of the drive ω m = ω d .
In the future, adaptive measurement schemes that change the drive frequency from shot to shot or between bunches of shots may be possible.This can lead to a further improvement of transmon qudit measurements.E 0 (0) = 0. We define the average transition frequency ω i,j between |i⟩ and |j⟩ of both configurations n g = 0, 1/2 as and the frequency difference ∆ω i,j as The anharmonicity α j and the energy dispersion ϵ j of the transmon qudit are defined by We numerically obtained the fundamental parameter E J /E C of a specific IBM Quantum backend by demanding that the qubit frequency ω 0,1 and anharmonicity α 1 match the values reported by this device.In Fig. 5, we plot the dependence of E n on n g for the five lowest states.The values of the frequency difference vary from ∆ω 0,1 /2π = 25.1 kHz to ∆ω 3,4 /2π = −142 MHz.In Fig. 1, the values ϵ 3 for a number of IBM Quantum devices are displayed.For large ϵ j compared to ω j,j+1 , the Wigner function of the state |j⟩ effectively is smeared out in phase space since both configurations n g = 0, 1/2 exhibit different resonance frequencies.

Appendix B: Schrieffer-Wolff Transformation
In this appendix, we will use Schrieffer-Wolff transformations [29] to obtain perturbative approximations of, first, the Jaynes-Cummings interaction between a qudit and a readout resonator and, second, the two-photon drive of the qudit.
In general, the system of interest can be defined by the Hamiltonian Here, H 0 and H 1 are block diagonal in the subsystems, whereas V is block off diagonal.To find an effective block-diagonal Hamiltonian, i.e., eliminate the block offdiagonal part V , the unitary transformation U = e S is applied to H, Expanding the anti-Hermitian operator S = −S † as S = ∞ n=1 λ n S (n) , H eff can be expressed as The first-order contribution reads H (1) = H 1 + V + S (1) , H 0 .To eliminate the block off-diagonal V in this expression, we impose V = − S (1) , H 0 .Since H 0 is block diagonal, S (1) has to be block off diagonal.As a consequence, S (1) , V is block diagonal.The secondorder contribution reads H (2) = 1 2 S (1) , V + S (1) , H 1 + S (2) , H 0 , and imposing S (1) , H 1 = − S (2) , H 0 guarantees the second order to be block diagonal.The secondorder contribution to the effective block-diagonal Hamiltonian is then given by We choose a superposition of all operators appearing in V as an ansatz for S (1) .

Jaynes-Cummings Interaction
Following the notation of [19], the Hamiltonian describing a transmon qudit coupled to a readout resonator reads where ω j is the energy (see Appendix A) of the bare qudit state |j⟩, ω r is the energy of the readout resonator, and a ( †) is its annihilation (creation) operator.The parameters g j,j+1 denote generalized Jaynes-Cummings coupling strengths between the qudit and the resonator.Using the approximation g j,j+1 = g √ j + 1, the interaction Hamiltonian reduces to g(a † b + ab † ) [18], where b ( †) is the annihilation (creation) operator of the transmon qudit.The qudit and resonator Hamiltonians H q and H r denote two blocks of commuting operators |i⟩⟨j| and a ( †) .The generalized Jaynes-Cummings interaction couples both blocks.
We start with identifying the block-diagonal and block off-diagonal parts, Using a superposition of all operators appearing in V as an ansatz for S (1) , the coefficients C j are obtained as Replacing C j in Eq. (B8) by this expression and using the definition g −1,0 = 0 as well as the sign convention of [1,19] leads to Here, and we have neglected terms proportional to (a 2 |j + 2⟩⟨j| + H.c.).This is justified by the possibility to interpret these terms as perturbations that are eliminated by a second Schrieffer-Wolff transformation.This will lead to terms proportional to |j⟩⟨j|, a † a|j⟩⟨j|, (a † a) 2 |j⟩⟨j|, and also (a 4 |j + 4⟩⟨j| + H.c.).Importantly, for typical values of g j,j+1 , ω j , and ω r , the coefficients of all these terms are a factor of 10 4 smaller than the previous second-order contributions and can therefore safely be neglected.In conclusion, in Eq. (B10), we arrived at corrections to the Hamiltonian that are diagonal in the qudit and resonator states.The shifts of the qudit and resonator energies are ωj = ω j + χ j−1,j , ωr,j = ω r + χ j . (B13) The resonance frequencies of the qudit transitions |i⟩ ↔ |j⟩ can be estimated to ωi,j = ωj − ωi j − i .(B14)

Qudit Drive
In analogy to Appendix B 1, we perform a Schrieffer-Wolff transformation in a system consisting of a qudit and its drive only.We will use this method to predict the Rabi oscillation frequencies of second-order transitions |j⟩ ↔ |j + 2⟩.For previous work on multiphoton transitions, see, e.g., [30,31].Starting with the Hamiltonian [1], an expansion in λ leads to where Note that these expressions only hold for ω d ̸ = ωj+1 − ωj , i.e., drives that are not resonant with transitions between neighboring qudit levels |j⟩ ↔ |j + 1⟩.

Appendix C: Preparation of Qudit States
We prepare the four lowest Fock states of an IBM Quantum transmon qudit, described in Section IV, by applying sequences of calibrated X gates to the qudit ground state |0⟩.For simplicity, we implement these X gates via Gaussian pulses.For each pulse, we first calibrate its drive frequency ω d and second, its drive amplitude Ω q .
The optimal drive frequency is obtained from a Gaussian fit to resonance measurement data shown in Fig. 6a, where we fix the pulse amplitude to an initial estimate.First, the measured N = 2000 shots per qudit drive frequency are averaged.Second, these averages are rotated such that their major principal axis is oriented along the X axis.And third, the means are projected onto the X axis which justifies the axis label "rotated projected data."For the spectroscopy measurements, in addition, we define the origin of the y axis of Fig. 6a to correspond to the initial state of the analyzed transition and the maximum to the final state.Our estimated frequency ωi,j is calculated by Eq. (B14) using g/2π = 65 MHz.
To obtain the initial estimate of the qudit drive amplitude, we define the rotation angle θ of a resonant Rabi Therefore, the rotation angle θ depends on the effective Rabi frequency Λ(j, k) and the pulse duration t.Following Eq. (B15) for λ = 1, the rotation angle θ for Rabi oscillations between |j⟩ and |j + 1⟩ is given by θ = t j,j+1 Ω (j,j+1) j + 1 .
After evaluating the resonance measurement, we calibrate the X-gate drive amplitude via Rabi oscillations, see Fig. 6b.The data are rotated and projected onto the major principal axis as described before for the spectroscopy measurements.As shown in Eq. (B16), the Rabi frequency for transitions |j⟩ ↔ |j + 2⟩ depends nonlinearly on the drive amplitude.Since these transitions are suppressed by the small factor Ω q f j , we choose t j,j+2 = 2t 0,1 such that Ω (0,2) π does not exceed the limits of IBM Quantum software/hardware restrictions.The π amplitude of an X-gate pulse is identified with the location of the first maximum in Fig. 6b, indicated by dashed lines.For transitions between neighboring states, we fit a sine dependence on a linear function of Ω q and, for second-order transitions, we fit a sine dependence on a second-order polynomial of Ω q .Using those fits, any desired rotation angle, e.g., π for an X gate or π/2 for a Hadamard gate, can be mapped back to a corresponding pulse amplitude.
The sequence of calibrating drive frequency and amplitude described above can be iterated several times to improve gate fidelity.For simplicity, we consider only one round of calibrations.To increase fidelity, we chose the initial value for Ω 1,3 based on prior test measurements.
Implementing gates in the |j⟩ ↔ |j + 2⟩ subspace results in two advantages.First, our implementation of an X-gate X j,j+2 between |j⟩ and |j + 2⟩ takes only twice the single-qudit gate duration t 0,1 .In contrast, using single-qudit gates, X j,j+2 consists of three single-qudit operations X j,j+1 X j+1,j+2 X j,j+1 with a total duration of 3t 0,1 .Second, the calibration of the drive frequencies (amplitudes) for |0⟩ ↔ |1⟩ and |0⟩ ↔ |2⟩ are independent of each other and can therefore be combined into a single Qiskit job (set of measurements submitted to an IBM Quantum device).In contrast, the frequency calibration for the transition |1⟩ ↔ |2⟩ depends on the Rabi measurement for the transition |0⟩ ↔ |1⟩.In total, we can perform our calibration in four Qiskit jobs: (i ) drive frequency of X 01 and X 02 , (ii ) drive amplitude of X 01 and X 02 , (iii ) drive frequency of X 12 , X 23 , and X 13 , and (iv ) drive amplitude of X 12 , X 23 , and X 13 .In contrast, the standard sequential calibration of single-qudit X gates would take six jobs: two for each of the three single-qudit X gates between neighboring states.

FIG. 2 .
FIG. 2. (a) Drive-frequency-dependent phase-space positions A (d) j and A (m) j of the coherent state of the resonator given the qudit prepared in |j⟩, see Eqs. (7) and(12).For ωm = ω d , the trajectories of all states A (d) j match, denoted by the black

FIG. 3 .
FIG. 3. Ratio of the standard deviation SDm for the multifrequency strategy and the standard deviation SDs of the single-frequency strategy applied to an equal-superposition state taking N = 1000 shots.The gray region indicates where both standard deviations exceed SD s/m ≥ 0.1.The straight lines denote constant values of σκ/Ω.We take the same qudit parameters as in Fig. 2 and choose g/2π = 100 MHz and Ω/2π = 100 MHz.
FIG. (a)Drive-frequency-dependent ququart measurements on ibm lagos Q4 (July 7, 2023), the experimental equivalent of Fig.2a.The colored curves (with white shadows) correspond to the centers of Gaussian fits to the Fock states |j⟩.Black straight lines indicate the boundaries of regions assigned to individual Fock states.These boundaries are constructed using the minimum distance estimator, see Eq. (16).For the drive frequency ω0 = 7.2463 GHz, we show the measurement results of all N = 2000 shots for each prepared Fock state.Black crosses mark the centers of their Gaussian fits.(b) Measurement errors ξj = 1 − Mj,j based on all measured shots of the data presented in (a).Here, the elements Mi,j of the assignment matrix equal the relative number of shots, Ni/N , that are classified as |i⟩ even if |j⟩ is prepared.The horizontal gray line denotes the measurement error ξ def obtained using the default measurement pulse.(c) Measurement errors ξj based on Gaussian fits to the data presented in (a) and the assignment matrix M defined in Eq. (14).Using the centers and average σ of Gaussian fits for each readout resonator drive frequency, we calculate M numerically.The ξj shown in (b) are larger than those in (c) since they do not only represent assignment errors, but also include additional errors such as qudit decay, leakage, and imperfect state preparation.In both (b) and (c), the minimum of the average assignment error ξ will be smaller than ξ def obtained by the default pulse if the modulation frequency is chosen such that ωm = ω d , cf.Figs.2b and c.
FIG. (a)Qudit resonator spectroscopy of transitions |i⟩ ↔ |j⟩.Colored markers denote measured data and solid curves correspond to Gaussian fits.We plot each resonance spectrum centered around the predicted transition frequency ωi,j, see Eq. (B14), using g/2π = 65 MHz.(b) Rabi oscillations |i⟩ ↔ |j⟩ depending on the drive amplitude for fixed pulse duration.We sweep the readout resonator drive amplitude while keeping all other parameters of the drive fixed.For first-order qudit state transitions, we fit a sinusoidal dependence on a linear function of the drive amplitude in the interval [0, 0.5], see the lines connecting crosses and, respectively, diamonds.For second-order qudit state transitions, we fit a sinusoidal dependence on a quadratic function of the drive amplitude in the interval [0, 1], see the lines connecting dots and, respectively, triangles.The vertical dashed lines indicate the locations of the first maxima obtained from the fits.
FIG. 1. Overview of qudit parameters of the IBM Quantum devices listed in the legend.The qubit resonance frequency and anharmonicity are denoted by ω0,1 and α1, see Appendix A. The energy dispersion ϵ3 of the third excited state given in Eq. (A5) follows from these device specifications that were accessed on May 23, 2023.The labeled straight black lines denote constant values of EJ /EC .