<figure><img src="image_1.png" /></figure>

HAL
open science

Dynamics of Transmon Ionization

Ross Shillito, Alexandru Petrescu, Joachim Cohen, Jackson Beall, Markus

Hauru, Martin Ganahl, Adam G.M. Lewis, Guifre Vidal, Alexandre Blais

► To cite this version:

Ross Shillito, Alexandru Petrescu, Joachim Cohen, Jackson Beall, Markus Hauru, et al. Dynamics of Transmon Ionization. Physical Review Applied, 2022, 18 (3), pp.034031. 10.1103/PhysRevApplied.18.034031. hal-03630643

HAL Id: hal-03630643

https://hal.science/hal-03630643v1

Submitted on 27 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dynamics of Transmon Ionization

Ross Shillito,¹,²,* Alexandru Petrescu,¹,³ Joachim Cohen,¹ Jackson Beall,² Markus Hauru,²

Martin Ganahl,² Adam G.M. Lewis,² Guifre Vidal,²,⁴,⁵ and Alexandre Blais¹,⁵

¹Institut Quantique and Département de Physique, Université de Sherbrooke, Sherbrooke, Quebec, J1K 2R1, Canada

²Sandbox@Alphabet, Mountain View, California 94043, USA

³Laboratoire de Physique de l'Ecole Normale Supérieure, Mines-Paristech, CNRS, ENS-PSL, Inria, Sorbonne Université, PSL Research University, Paris 75005, France

⁴Google Quantum AI, Mountain View, California 94043, USA

⁵Canadian Institute for Advanced Research, Toronto, Ontario M5G 1M1, Canada

<figure><img src="image_3.png" /></figure>

(Received 21 April 2022; revised 20 July 2022; accepted 11 August 2022; published 13 September 2022)

Qubit measurement and control in circuit quantum electrodynamics (QED) rely on microwave drives, with higher drive amplitudes ideally leading to faster processes. However, degradation in qubit coherence time and readout fidelity has been observed even under moderate drive amplitudes corresponding to a few photons populating the measurement resonator. Here, we numerically explore the dynamics of a driven transmon-resonator system under strong and nearly resonant measurement drives and find clear signatures of transmon ionization where the qubit escapes out of its cosine potential. Using a semiclassical model, we interpret this ionization as resulting from resonances occurring at specific resonator-photon populations. We find that the photon populations at which these spurious transitions occur are strongly parameter dependent and that they can occur at low resonator-photon population, something that may explain the experimentally observed degradation in measurement fidelity.

DOI: 10.1103/PhysRevApplied.18.034031

I. INTRODUCTION

Dispersive readout in circuit quantum electrodynamics (QED) is realized by driving a measurement resonator coupled to the qubit [1]. In principle, increasing the drive amplitude, and thereby the resonator-photon population, increases the measurement rate, something that is expected to lead to fast, high-fidelity, and quantum nondemolition (QND) readout [2]. However, experimentally the fidelity and QND character of the measurement of transmon qubits [3] is often observed to decrease beyond a photon-number threshold, which can be as small as a few photons [4,5]. Perturbative models have been made in attempts to explain these observations but have limited applicability beyond low photon numbers [6–8].

In this paper, we go beyond perturbative treatments by numerically investigating the full dynamics of a strongly driven transmon-resonator system. At distinct resonator

populations, we find clear signatures of transmon ioniza-
tion where transmon states above the Josephson-junction
potential are occupied [9,10]. Because these states are not
strongly influenced by the Josephson potential, they are
well described by charge states. Consequently, for states
above the transmon well, the transmon-resonator coupling
appears longitudinal and the system dynamics are con-
sequently modified. Importantly, the term *ionization* is
used loosely here and refers only to the escape of the
transmon population to states lying above the Josephson
cosine-potential well.

Accurate simulation of the dynamics of transmon ion-
ization requires description of the density matrix of the
system on a truncated Hilbert space of very large dimen-
sion, something that is made possible here by the use
of large-scale computational accelerators known as tensor
processing units (TPUs). While other studies have been
limited to steady-state calculations with strongly detuned
drives [9], the computational power of TPUs allows us
to simulate the full time dependence with drives that are
resonant with the resonator, as is relevant for qubit mea-
surement. Accounting for tens of transmon levels and hun-
dreds of resonator states, we moreover see signatures of the
high-power readout [11–13]. We interpret these numerical
results using a semiclassical model capturing the nonlinear

<footnote>*Ross.Shillito@USherbrooke.ca</footnote>

<footnote>Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license. Fur-
ther distribution of this work must maintain attribution to the
author(s) and the published article's title, journal citation, and
DOI.</footnote>
response of the driven system. Using this model, we identify parameter regimes where ionization is expected to occur at sufficiently small photon number to affect dispersive readout, observations that are in qualitative agreement with experiments [4].

This paper is organized as follows. In Sec. III, we introduce the model, provide details on its TPU implementation, and present the numerical results on the dynamics of transmon ionization. Next, in Sec. III, we formulate a semiclassical theory that allows us to interpret the dynamics of the coupled transmon-resonator system as governed by transitions between qubit-state-dependent effective resonators (Sec. III A) obeying nonlinear equations of motion (Sec. III B). In Sec. III C, we choose system parameters that illustrate how ionization can reduce readout fidelity even at the low drive powers that are typical of dispersive readout in circuit QED. We summarize our findings in Sec. IV.

## II. MASTER-EQUATION SIMULATIONS WITH TENSOR PROCESSING UNITS (TPUs)

### A. Model

We consider a transmon capacitively coupled to a resonator [see Fig. 1(a)]. In the presence of a drive of amplitude $\mathcal{E}$ and frequency $\omega_d$ on the resonator, the Hamiltonian of the system takes the usual form ($\hbar = 1$) [2]

$$ \hat{H}(t) = 4E_c \hat{n}_t^2 - E_J \cos \hat{\phi}_t + \omega_r \hat{a}^\dagger \hat{a} - ig\hat{n}_t(\hat{a} - \hat{a}^\dagger) - i\mathcal{E}(\hat{a} - \hat{a}^\dagger) \sin(\omega_d t). \quad (1) $$

The first two terms correspond to the free transmon Hamiltonian with charging energy $E_c$, Josephson energy $E_J$, charge operator $\hat{n}_t$, and phase operator $\hat{\phi}_t$. We denote the eigenenergies and eigenstates of the free transmon Hamiltonian $E_i$ and $|i\rangle$, respectively. The first two of those eigenstates, labeled $\{|0\rangle, |1\rangle\}$, span the computational basis of the qubit. Of the higher excited states, an approximate number $2E_J/\omega_p$ are bound states lying within the cosine-potential well illustrated in Fig. 1(b). Here, $\omega_p = \sqrt{8E_cE_J}$

<figure><img src="image_4.png" /><figcaption>FIG. 1. (a) A transmon (green) capacitively coupled to a driven resonator (blue). (b) The cosine potential of the transmon with its first nine eigenstates. For our choice of parameters, the last three states are ionized states lying above the top of the potential.</figcaption></figure>

is the plasma frequency, which is approximately the 0–1 transition frequency of the transmon [3]. Moreover, we label as $|n\rangle$ the eigenstates of the free-resonator Hamiltonian of frequency $\omega_r$, corresponding to the bosonic annihilation operator $\hat{a}$. The transmon-resonator coupling of amplitude $g$ in the second line of Eq. (1) includes fast-rotating terms that are beyond the rotating-wave approximation (RWA) and are important to capture the contribution of high-energy states [14]. In the absence of the drive, the dressed energies and states of the coupled system are denoted $E_{\bar{i},n}$ and $|i, n\rangle$.

Including cavity loss at a rate $\kappa$, the driven transmon-resonator system is described by the usual Lindblad master equation [2],

$$ \dot{\hat{\rho}} = -i[\hat{H}(t), \hat{\rho}] + \kappa \mathcal{D}[\hat{a}]\hat{\rho}, \quad (2) $$

with the dissipator

$$ \mathcal{D}[\hat{O}]\hat{\rho} = \hat{O}\hat{\rho}\hat{O}^{\dagger} - \frac{1}{2}\{\hat{O}^{\dagger}\hat{O}, \hat{\rho}\}. \quad (3) $$

Because we are interested in capturing the dynamics of the system in the presence of a large-amplitude nearly resonant drive on the resonator, leading to several hundred of photons and highly excited states of the transmon, we keep up to 32 states in the transmon and 1024 states in the resonator for the most demanding simulations, corresponding to a total Hilbert-space dimension $2^{15} = 32768$. That is, the joint density matrix $\hat{\rho}$ for the transmon-resonator system is a Hermitian matrix of size $2^{15} \times 2^{15}$, which thus contains $2^{30} = 1073741824$ time-dependent complex coefficients. Furthermore, given that the unbound transmon states are approximately eigenstates of the charge operator, their eigenvalues increase quadratically rather than linearly, as is the case for the bound states. This increases the complexity of numerical simulations, as larger eigenvalues require smaller integration step sizes for convergence. Additional details, including a prescription for adapting the resonator Hilbert-space size, are given in Appendix A.

In order to perform these challenging numerical simulations, we resort to TPUs. Google’s TPUs are application-specific integrated circuits designed exclusively to accelerate large-scale machine-learning workloads [15]. Recently, they have been repurposed for other high-performance computational tasks [16–28], including simulations of quantum systems in large Hilbert spaces [24–28]. Due to their (i) matrix-multiply units capable of accelerating matrix multiplication, (ii) large amounts of high-bandwidth memory, and (iii) fast intercore interconnects directly connecting up to thousands of cores, TPUs are particularly fast at performing large-scale dense linear-algebra operations, which are required, e.g., to numerically integrate the above Lindblad master equation. More details on TPUs can be found in Appendix B.
For a given drive amplitude $\mathcal{E}$, simulations begin with the system initialized in either the dressed ground $|0,0\rangle$ or excited $|1,0\rangle$ states. To reduce the simulation time while capturing the transient dynamics including the ionization of the transmon, we simulate the evolution over a time of at least $\kappa^{-1}$. After each period of the drive, the reduced transmon, $\hat{\rho}_t$, and resonator, $\hat{\rho}_r$, density matrices are recorded.

Unless otherwise stated, in all of our simulations as well as in the semiclassical model discussed in Sec. III, we use the parameters $E_J/E_C = 50$ and $E_C/h = 280$ MHz, corresponding to a transmon of frequency $\omega_t/2\pi = 5.304$ GHz and anharmonicity $\alpha/2\pi = 322$ MHz. These parameters result in approximately six bound transmon states within the cosine potential (see Fig. 1). The resonator frequency is $\omega_r/2\pi = 7.5$ GHz and $g/2\pi = 250$ MHz. A resonator loss rate of $\kappa/2\pi = 20$ MHz is chosen to ensure fast resonator dynamics and the drive amplitude takes values in the range $\mathcal{E}/2\pi \in [0, 440]$ MHz. We fix the drive frequency to the bare-resonator frequency, $\omega_d = \omega_r$. The above parameters place the system in the dispersive regime, with $\chi/2\pi = -5.6$ MHz and $\chi/\kappa = 0.28$, where the dispersive shift is defined here as $\chi = (E_{11} - E_{10} - E_{01} + E_{00})/2$. Note that the resonator population is much smaller than the critical photon number $n_{\text{crit}} = (\Delta/2g')^2 \approx 15$, with

$g' = (E_J/32E_C)^{1/4}g$ [2]. Finally, to minimize the number of simulations that are required, we note that we ignore a possible offset gate charge on the transmon. The effect of the gate charge on the ionization will be studied in a future work [29].

## B. Numerical results

Figures 2(a) and 2(b) illustrate the average photon number $N_r = \langle \hat{a}^\dagger \hat{a} \rangle$ (red lines) and the average transmon population $N_t = \sum_i i \langle i | \hat{\rho}_t | i \rangle$ (black lines) as a function of time for the drive amplitude $\mathcal{E}/2\pi = 280$ MHz. The system is initialized in the dressed ground state $|0,0\rangle$ in Fig. 2(a) and in the dressed excited state $|1,0\rangle$ in Fig. 2(b). Also shown is the instantaneous distribution $N_{i,t} = \langle i | \hat{\rho}_r | i \rangle$ of the transmon states (color scale). The difference between the two initial states is striking. For this drive amplitude, when starting in the ground state, the transmon leaks out of its initial state but the distribution largely remains confined within the cosine potential of the Josephson junction. In contrast, when the transmon is initialized in its excited state, we observe a sudden jump of the average transmon population with a distribution extending well above the top of the cosine potential. This is a clear illustration of the ionization of the transmon. For the drive amplitudes

<figure><img src="image_5.png" /><figcaption>FIG. 2. (a),(b) The transmon and photon population as a function of time with $\mathcal{E}/2\pi = 0.28$ GHz for (a) the ground state and (b) the excited state. The black line indicates the average transmon population $N_t$ (left axis) and the solid red line depicts the resonator average population $N_r$ (right axis). The full distribution $\log_{10}(N_{i,t})$ of the transmon population for each transmon level $|i\rangle$ is also plotted, with the color bar above providing the scale. (c),(d) The purity of the reduced transmon density matrix as a function of the drive amplitude and time for the transmon initialized (c) in the ground state and (d) in the excited state. The dashed red lines indicate when the resonator reaches a population of 105 and 42 photons for the ground and excited states, respectively. The white dashed lines indicate the drive amplitude $\mathcal{E}/2\pi = 0.28$ GHz of (a) and (b). Because the drive amplitude is not sufficiently large to reach a population of $N_r = 105$ photons when the qubit is the ground state, ionization is not observed in (a).</figcaption></figure>
considered here, the simulated measurement is far from QND and the dynamics are therefore not well described by a dispersive Hamiltonian [2].

As illustrated in Figs. 2(c) and 2(d), transmon ionization is accompanied by a sudden drop in the purity $P = \text{Tr}[\hat{\rho}_t^2]$ of the reduced transmon density matrix. The rapid decline of the purity is observed at specific resonator-photon populations indicated by the dashed red line: $N_r = 105$ when initialized in the ground state in Fig. 2(c) and $N_r = 42$ for the excited state in Fig. 2(d). The dashed red lines terminate whenever the drive is too weak for the resonator to reach those populations. That transmon ionization occurs at specific photon populations suggests that the phenomenon is due to resonances, which we discuss in more detail below. These observations are compatible with Ref. [9], where steady-state calculations with an off-resonant drive have also shown drops in transmon purity in steady-state numerical calculations. Resonances at large photon numbers have also been observed in Ref. [14].

As a further illustration of transmon ionization at distinct photon numbers, we plot in Fig. 3(a) the transmon population $N_t$ as a function of the resonator population $N_r$, both values taken from time traces such as shown in Figs. 2(a) and 2(b). The different curves correspond to different drive amplitudes and the initial state of the transmon is identifiable from the starting transmon population. Remarkably, the responses essentially collapse to single curves for each of the initial states. Below the ionizing photon population, the transmon population remains close to the computational manifold and the resonator population exhibits transient oscillations due to the drive being off resonant with the Lamb-shifted resonator frequency. Together with Purcell decay, these transient oscillations are responsible for the features observed at small $N_r$. Above the ionizing photon populations, the transmon population rapidly increases. This coincides with the sudden drop in the purity observed in Figs. 2(c) and 2(d).

Focusing now on the dynamics of the resonator, Fig. 4 shows the Wigner function of the resonator at different times with the transmon starting in its excited state and with the same drive amplitude $\mathcal{E}/2\pi = 0.28$ GHz as used in Figs. 2(a) and 2(b). For small features to appear more clearly, we plot the logarithm of the absolute value of the Wigner function. As expected, the cavity is initially well described by the single coherent states $|\alpha_1\rangle$ associated with the first excited state of the transmon [see Fig. 2(a)]. As the transmon population in other levels increases, additional features emerge corresponding to the coherent states associated with the now occupied transmon levels. Additionally, “bananization” caused by transmon-induced nonlinearities becomes apparent [30,31]. The dashed lines correspond to fixed photon numbers and are used as a guide to the eye. The full colored lines overlying the Wigner

<figure><img src="image_6.png" /><figcaption>FIG. 3. (a) A parametric plot of the average transmon population versus the average resonator population for an evolution time of $\kappa^{-1}$. Each line represents a unique simulation with a different drive amplitude $\mathcal{E}$ and starting either in the ground or excited state. These metrics are calculated in the bare basis and are defined as in Sec. II B. (b) The average transmon population as a function of the bare-resonator population for each resonator branch $\{i\}$. Each branch is represented by a different color, with dashed gray lines indicating branches that do not feature significantly in the dynamics. The nearly vertical dark lines represent the overlap $C_i(n)$ between the different branches, darker lines indicating stronger overlap.</figcaption></figure>

functions are obtained from a semiclassical approximation that we now introduce.

### III. STATE IDENTIFICATION AND SEMICLASSICAL INTERPRETATION

The TPU-based large-scale numerical simulations presented in Sec. clearly illustrate the breakdown of the dispersive approximation. Notably, at specific resonator-photon populations, we observe a sudden jump in transmon population above the cosine-potential well, which is associated with a sharp drop in transmon purity. This results in complex dynamics of the resonator field, as illustrated by its Wigner function. In this section, we develop a semiclassical model to understand the main features of these observations.
<figure><img src="image_7.png" /><figcaption>FIG. 4. The logarithm of the absolute value of the Wigner function $log_{10}(|W|)$ of the resonator at different times for a drive amplitude $E/2\pi = 0.28$ GHz and the transmon initialized in its dressed excited state. The logarithm is chosen to reveal the smaller features. Overlaid are the semiclassical flows of the relevant branches following the color convention of Fig. 3(b). The dashed circles are a guide to the eye at constant photon numbers (50, 100, 150, and 200).</figcaption></figure>

A. Identification of dressed states and resonator branches

Our semiclassical approach is based on the dressed states $|i, n\rangle$ and energies $E_{i, n}$ of the transmon-resonator Hamiltonian in Eq. (1) in the absence of the drive. More precisely, at arbitrary photon number, we identify the dressed states that are closest to transmon eigenstates. At low photon number, this identification is simple in the dispersive regime. However, as the photon population approaches and even exceeds $n_{crit}$, the dressed states are highly entangled and identification becomes difficult.

Building on Ref. [13], our approach relies on first identifying the eigenstates of the Hamiltonian of Eq. (1) for $E = 0$ obtained from numerical diagonalization with the largest overlap with $|i, 0\rangle$, the bare transmon states at zero photon population. Then, given an identified state $|i, n\rangle$ for $n \ge 0$, the next state $|i, n+1\rangle$ is chosen from the subset of remaining eigenstates $\{|\lambda\}\}$ such that the overlap

$$C_i(n) = \left| \frac{\langle \lambda | \hat{a}^\dagger | i, n \rangle}{\langle i, n | \hat{a} \hat{a}^\dagger | i, n \rangle} \right|^2 \quad (4)$$

is maximized. Following this procedure recursively, we obtain a set of states $\{|i, n\}$, where the bare transmon label $i$ is fixed and $n$ spans a desired range of resonator population. We refer to each such set of states as a branch $\{i\}$ of

the resonator. In Fig. 3(b), we plot the average transmon population of the first 16 of these resonator branches as a function of the photon number. The full colored lines are branches that play an important role in the understanding of the numerical results. Branches playing a more minor role for our particular choice of parameters are illustrated as dashed gray lines. At finite transmon-resonator coupling and because we include non-RWA terms in the system Hamiltonian, the different branches have a nonzero overlap $C_i(n)$.

To illustrate this, in Fig. 3(b) we plot, with dark and nearly vertical lines, the overlap $C_i(n)$ between branches whenever it rises above 0.01. Darker lines indicate a stronger overlap between the states. At small photon number, the overlap between $\{0\}$ and $\{1\}$ corresponds to Purcell decay, which is seen to decrease with increasing photon number [32]. More importantly, we note a strong overlap between the branches $\{0\}$ and $\{9\}$ for $N_r \sim 110$ and between $\{1\}$ and $\{5\}$ for $N_r \sim 50$, in agreement with the photon numbers at which sudden drops of transmon purity are observed in Figs. 2(c) and 2(d) (red dashed lines). To further emphasize the link between the large overlaps and transmon ionization, we reproduce as solid blue and red lines the first two branches $\{0\}$ and $\{1\}$ of Fig. 3(b) together with the numerical data in Fig. 3(a). The sharp elbows, which are indicative of the fact that transmon
ionization results from a resonance, align well between the semiclassical results and the numerical data.

A strong overlap is also observed between {0} and {8} close to $N_r = 30$ [see Figs. 3(a)and 3(b)]. In contrast to what is observed at $N_r \sim 50$ and 110, this sharper resonance results in little population transfer, seen as a small drop of purity for $\mathcal{E} \in (0.2 - 0.3)/2\pi$ GHz in Fig. 2(c). This can be understood from the framework of Landau-Zener transitions. In our semiclassical picture, as the resonator rings up toward its transmon-state-dependent steady state, the photon number is swept at a rate related to $\kappa$ and to the drive amplitude. A large drive amplitude leads to a rapid sweep through the narrow feature at $N_r \sim 30$, resulting in a diabatic passage with little leakage [see the darker-blue numerical lines in Fig. 3(b)]. We note that similar effects have been observed in the study of a resonator-induced phase gate, in which qubit-resonator leakage is maximized near avoided crossings [33]. On the other hand, for our choice of parameters, the features at $N_r \sim 50$ and 110 are broader and the sweep rate therefore comparatively slower. This slower sweep rate results in a nondiabatic process, leading to leakage out of the transmon computational subspace and to the observed drop in purity. This interpretation is further confirmed in Appendix C for numerical experiments with a different drive frequency and loss rates. In the steady state, the above concepts are not applicable and ionization is determined only by the chosen parameters and drive amplitude [9].

## B. Semiclassical dynamics

Following Ref. [13], we assign to each branch {$i$} an effective oscillator of photon-number dependent frequency

$$ \omega_i(n) = E_{\overline{i,n+1}} - E_{\overline{i,n}} \quad (5) $$

As seen in Fig. 5, where we plot $\omega_i(n)$ versus $n$, Eq. (5) accounts for the photon-number dependence of the resonator frequency, including that the resonator responds at its bare frequency at large photon numbers. Assuming that the transmon remains in a given branch {$i$}, the dynamics of this effective oscillator approximately obey the classical equation of motion of a driven damped oscillator,

$$ \dot{\alpha}_i = -i[\omega_i(|\alpha_i|^2) - \omega_d]\alpha_i - i\mathcal{E}/2 - \kappa\alpha_i/2, \quad (6) $$

where nonlinear effects are encapsulated in the dependence of the branch frequency on the photon number $\omega_i(|\alpha_i|^2)$ according to Eq. (5). Because the quantity $|\alpha_i|^2$ takes arbitrary real values, we generalize Eq. (5) from discrete values $n$ to a continuous function. To do this, we smooth $\omega_i(n)$ with a first-order Savitzky-Golay filter and linearly interpolate between each $n$. Doing this additionally removes large discontinuities in the photon-number dependence of the effective frequencies $\omega_i(n)$ caused by

<figure><img src="image_8.png" /><figcaption>FIG. 5. The effective frequency $\omega_i(N_r)$ of the resonator branches as a function of the resonator population. The frequency is smoothed using a Savitzky-Golay filter to remove discontinuities caused by strong interactions with other branches. The dashed line indicates the bare-cavity frequency $\omega_r/2\pi = 7.5$ GHz.</figcaption></figure>

strong interactions with other branches. This semiclassical approximation is expected to accurately describe the system for large photon number $|\alpha_i|^2 \gg 1$.

Numerically integrating Eq. (6), we plot in Fig. 4(a) the time dependence of $\alpha_1(t)$ (red lines) together with a snapshot of the Wigner function at time $t = 0.13\kappa^{-1}$ assuming that the effective oscillator starts in vacuum. The arrows indicate the flow of time and the different lines are obtained from Eq. (6) for a set of initial conditions in an area close to vacuum representing zero point fluctuations. The red dot corresponds to the steady-state value $\alpha_1^s$.

To account for transitions between transmon states, the evolution of $\alpha_i(t)$ associated with other branches is illustrated whenever the semiclassical evolution reaches photon numbers at which we anticipate non-negligible rates for transitions into those branches according to the overlaps observed in Fig. 3(b). To distinguish between the different branches, we use the color scheme of Fig. 3(b). For example, the blue lines in Figs. 4(a)–4(d) correspond to branch {0}, which appears as a result of Purcell decay. In the same way, starting in Fig. 4(b), we include the flow of $\alpha_5(t)$ associated with branch {5} (brown lines), which has a strong overlap with {1} when $N_r \sim 50$. Following the evolution of the Wigner function from panel to panel, it is possible to see a feature following the flow of $\alpha_5(t)$ and settling at the expected steady-state value $\alpha_5^s$ in Fig. 4(d) (brown dot). Applying this procedure whenever the resonator number is such that an occupied branch has a strong overlap $C_i(n)$ with another branch, the vast majority of the Wigner-function transient behavior can be understood. In particular, starting in Fig. 4(d), we see the bistable behavior of {$\alpha_1$} in both the semiclassical results (red lines) and the numerical data.

Rather than focusing on the steady-state response, in Fig. 6(a) we plot, for different drive amplitudes and as a
function of time, the difference in the resonator populations $\Delta N_r$ obtained by solving the semiclassical expression in Eq. (6), given that the transmon is initialized in its ground or excited state. As expected from Fig. 5, at some threshold power, the photon number rapidly increases if the qubit is initially in its excited state, while the increase is not as pronounced for ground state. This results in the observed large $\Delta N_r$. Figure 6(b) shows the same quantity obtained from numerical simulations. Given that the population of the resonator $N_r \ge 450$ for the strongest drive amplitude $\mathcal{E}/2\pi = 440$ MHz, these simulations are particularly demanding and require Hilbert-space sizes up to $2^{15}$. The agreement between the full numerical simulations and the simple semiclassical model is remarkable. In particular, in both approaches, we observe that the photon-number difference $\Delta N_r$ goes to zero at the strongest drive amplitudes. This is expected from Fig. 5, where at very large photon numbers, the frequencies $\omega_i(|\alpha_i|^2)$ eventually collapse to the bare value $\omega_r$. The different drive amplitude at which this collapse occurs depend on the initial transmon state and the resulting large $\Delta N_r$ is exploited in the high-power qubit readout [11–13]. The choice of drive amplitude range $\mathcal{E}/2\pi \in [0, 440]$ MHz for the TPU simulations is made to demonstrate the full range of the high-power readout behavior and demonstrate the decrease in $\Delta N_r$ at the largest amplitudes.

A possible interpretation for the observed large response at the bare-resonator frequency is that once the transmon is ionized, mostly chargelike states are occupied. Because these states couple longitudinally to the resonator, they do not lead to a resonator frequency pull [10]. However, the numerical results show that even for the highest drive powers, a significant distribution of states inside the well remain populated (not shown). As a result, the collapse of the resonator to its bare frequency cannot be explained as resulting alone from the longitudinal-type coupling of the

<figure><img src="image_9.png" /><figcaption>FIG. 6. The difference in the resonator population $\Delta N_r$ given that the transmon is initialized in its ground or excited state as a function of the drive amplitudes and time as obtained from (a) the semiclassical model and (b) the numerical TPU data.</figcaption></figure>

unbound states. We leave a more detailed analysis of this effect to future work [29].

### C. Resonances at low photon number

For the parameters used above, the first resonance leading to non-QNDness occurs at a relatively large photon number of approximately 30, followed by resonances at even larger photon numbers. Such resonances can, however, occur at much lower photon numbers. To demonstrate this, we now use parameters based on the experiment of Ref. [4]: $E_C/h = 314$ MHz, $E_J/E_C = 55.47$, $g/2\pi = 211$ MHz, and $\omega_r/2\pi = 4.804$ GHz, corresponding to $n_{crit} \approx 10$. The cavity damping rate is set to $\kappa/2\pi = 40$ MHz. The drive amplitude is varied in the range $\mathcal{E}/2\pi \in [0.02, 0.14]$ MHz and the evolution time is set to 48 ns, corresponding to the smallest measurement time used in Ref. [4]. We note that these are bare system parameters chosen such as to approach the dressed parameters reported in Ref. [4].

In a similar fashion to Fig. 3(a), in Fig. 7(a) we parametrically plot the transmon population $N_t$ against the resonator population $N_r$ for the transmon initialized in the excited state. For this choice of parameters and integration time, Purcell decay is apparent for the three lowest drive amplitudes (see the sharp decrease of $N_t$). For larger

<figure><img src="image_10.png" /><figcaption>FIG. 7. (a) A parametric plot of the average transmon population versus the average resonator population for an evolution time of 48 ns. Each line represents a unique simulation with a different drive amplitude $\mathcal{E}$ and starting in the excited state. Because of the small Hilbert-space size, these results are obtained using CPUs. (b) The average transmon population as a function of the bare-resonator population for each resonator branch $\{i\}$. The dashed lines indicate the overlap between different branches when $C_i(n) > 0.01$.</figcaption></figure>
measurement amplitudes corresponding to an average resonator population of approximately 2.5 photons, upward transition of the transmon population is clearly observed. Interestingly, this corresponds to the average photon number at which a decrease in measurement fidelity is observed in Ref. [4] [see Fig. 3(b) of that reference].

To understand the origin of this population leakage under measurement, we show in Fig. 7(b) the average photon population for the resonator branches {i} as obtained from our semiclassical model. At $N_r \approx 5$, the transmon population rapidly rises for branch {1} associated with the first excited state of the transmon (red line). In the numerical simulation shown in Fig. 7(a) as well as in experiments, the resonator field is in a coherent state. As a result, because of the $\sqrt{n}$ fluctuations of coherent states, the resonance at $N_r \approx 5$ already plays a role at $N_r \approx 2.5$, leading to the observed non-QNDness. While care must be taken when using a semiclassical model at such a small photon population, these observations are suggestive of the fact that the resonances that are observed here could play a role in the drop of measurement fidelity and QNDness that is experimentally observed at low photon number.

## IV. CONCLUSIONS

We leverage the computational power of TPUs to perform large Hilbert-space size time-dependent master-equation simulations of transmon readout in circuit QED. From these simulations, we identify resonances occurring at a distinct resonator-photon number where the transmon population escapes to states above the Josephson cosine-potential well. To interpret these results, we develop a semiclassical model capturing the nonlinear transmon-state-dependent change of the resonator frequency with photon number. In particular, this model correctly captures which states play a significant role in transmon ionization. Using a different set of parameters, we show that these resonances can occur at small photon population. These results suggest that the non-QND nature of the dispersive readout experimentally observed at small drive amplitudes could be due to these resonances. The location of ionization events can be predicted from the semiclassical branch analysis, thus providing an upper bound on the maximum population of the measurement photons to maintain the QND character of the readout. Moreover, because the location of these ionizing events is strongly parameter dependent, this work suggests a pathway to avoid these spurious effects when optimizing device parameters for readout.

## ACKNOWLEDGMENTS

We thank Cristóbal Lledó and Adrien Moulinas for insightful discussions. This work was undertaken in part thanks to funding from the Natural Sciences and Engineering Research Council (NSERC), the Canada First

Research Excellence Fund, the Ministère de l'Économie et de l'Innovation du Québec, the U.S. Army Research Office under Grant No. W911NF-18-1-0411, and the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers and Quantum Systems Accelerator. This research was also supported by Cloud TPUs from Google's TPU Research Cloud (TRC). G.V. is a Distinguished Invited Professor at the Institute of Photonic Sciences (ICFO) and a Distinguished Visiting Research Chair at Perimeter Institute. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research, Innovation and Science.

## APPENDIX A: NUMERICAL IMPLEMENTATION

To simulate the master equation in Eq. (2), we move to a rotating frame defined by

$$ \hat{U}_{\text{rf}}(t) = \exp(i\omega_d \hat{a}^{\dagger}\hat{a}) \quad (\text{A1}) $$

In that frame, the system Hamiltonian takes the form

$$ \begin{aligned} \hat{H}_I(t) &= \hat{U}_{\text{rf}}(t)\hat{H}(t)\hat{U}_{\text{rf}}^{\dagger}(t) - i\hat{U}_{\text{rf}}(t)\hat{U}_{\text{rf}}^{\dagger}(t), \\ &= \hat{H}_{I0} + \hat{X}_I(t), \end{aligned} \quad (\text{A2}) $$

with

$$ \hat{H}_{I0} = 4E_c \hat{n}_t^2 - E_J \cos(\hat{\phi}_t) + (\omega_r - \omega_d) \hat{a}^\dagger \hat{a} - \frac{\varepsilon}{2} (\hat{a}^\dagger + \hat{a}) \quad (\text{A3}) $$

and

$$ \hat{X}_I(t) = -ig\hat{n}_t(\hat{a}^\dagger e^{i\omega_d t} - \hat{a}e^{-i\omega_d t}) + \frac{\varepsilon}{2}(\hat{a}^\dagger e^{2i\omega_d t} + \hat{a}e^{-2i\omega_d t}). \quad (\text{A4}) $$

Moreover, in this rotating frame the master equation reads

$$ \dot{\hat{\rho}}_I(t) = -i[\hat{H}_I(t), \hat{\rho}_I] + \kappa \mathcal{D}[\hat{a}] \hat{\rho}_I \equiv \hat{\mathcal{L}}\hat{\rho}_I, \quad (\text{A5}) $$

with $\hat{\rho}_I(t) = \hat{U}_{\text{rf}}(t)\hat{\rho}(t)\hat{U}_{\text{rf}}^{\dagger}(t)$ and where the dissipator $\mathcal{D}[\hat{a}]$ is unaffected by the transformation.

To solve this master equation on TPUs, we approximate the action of the time-ordered Lindbladian exponential $\mathcal{T} \exp(\int_0^t \hat{\mathcal{L}}(t')dt')$. To do this, we first solve for the roots $\{z_i\}$ of the Taylor approximation of the exponential to nth order and rearrange the terms:

$$ e^x \approx \sum_{i=0}^{n} \frac{x^i}{i!} = \frac{1}{n!} \prod_{i=1}^{n} (x - z_i) = \prod_{i=1}^{n} (1 - x/z_i), \quad (\text{A6}) $$
where $\prod_i z_i = n!$. We then evolve the system by time $\delta t$ by iterating through the $n$ roots:

$$ \hat{\rho}(t + \delta t) = \hat{\rho}_n(t), \quad \hat{\rho}_{i=0}(t) = \hat{\rho}(t) \\ \hat{\rho}_i(t) = \left[ 1 - (\delta t/z_i) \overline{\hat{\mathcal{L}}(t, t + \delta t)} \right] \hat{\rho}_{i-1}(t), \quad (A7) $$

where $\overline{\hat{\mathcal{L}}(t, t + \delta t)}$ is calculated using a second-order Magnus expansion. The order $n$ and step size $\delta t$ required for convergence depend on the number of transmon and resonator states, with values of $n = 8-15$ and $1/(\omega_d \delta t) \approx 50-100$ found to be optimal for our choice of parameters.

The Hilbert space of the cavity is adapted throughout the evolution. At each step, the quantity $E = |1 - \text{Tr}\{\rho(N\delta t)[\hat{a}, \hat{a}^\dagger]\}|$ is calculated to ensure minimal occupation of the highest cavity excited states. The cavity size is doubled if $E > 10^{-6}$, with the previous step recalculated with the larger Hilbert-space size once this condition is met.

In the case of a transmon-resonator system driven off resonantly, a simple frame change can be used to follow the linear response of the oscillator, subsequently reducing the required Hilbert space [9]. In the case of a resonant drive, such a transformation does not exist.

## APPENDIX B: TENSOR PROCESSING UNITS

Google's TPUs are application-specific integrated circuits (ASICs) originally designed to accelerate and scale up machine-learning workloads [15]. By leveraging the JAX library [34], it is possible to repurpose TPUs to also accelerate other large-scale computational tasks [16–28]. For instance, in Refs. [24–27], TPUs are used to simulate the wave function of up to 36–40 (two-level) qubits. In this work, we use the power of TPUs to simulate, instead, the time evolution of the joint density matrix of a transmon-resonator system, as given by a Lindblad master equation.

We employ TPUs of the third generation, denoted v3. Each single TPU v3 core is equipped with two matrix-multiply units (MXUs) to formidably accelerate matrix-matrix multiplication (matmul), resulting in about ten teraflops (TFLOPs) of measured single-core matmul performance in single precision.

The smallest available TPU configuration consists of eight TPU v3 cores with a total of 128 GB of dedicated high-bandwidth memory (HBM), controlled by a single host with 48 CPU cores. The largest configuration is a pod with 2048 TPU v3 cores and 32 TB of HBM, controlled by 256 hosts. Given a choice of configuration, the available TPU cores are directly connected to nearest neighbors in a two-dimensional torus network through fast inter-core interconnects (ICIs). The ICIs play an essential role in the ability of the TPUs to maintain high performance when distributing matmuls and other dense linear-algebra

operations over all available TPU cores. In this work, we use the JAX library [34] to write single-program multiple data (SPMD) code and execute it on configurations made of multiple TPU cores. Specifically, for the largest density matrix under consideration, of size $2^{15} \times 2^{15}$, we use 128 TPU v3 cores. The density matrix is distributed over all available cores and updated according to the different terms in the Lindblad operator, with a typical update $\hat{\rho}(t) \rightarrow \hat{\rho}(t + \delta t)$ for $n = 15$ in Eq. (A7) taking on the order of seconds.

## APPENDIX C: LANDAU-ZENER TRANSITIONS FOR DIFFERENT PARAMETER CHOICES

Figure 8(a) shows the average transmon population versus the average resonator population for the same

<figure><img src="image_11.png" /><figcaption>FIG. 8. Parametric plots of the average transmon population versus the average resonator population for an evolution time of $\kappa^{-1}$ for (a) $\kappa/2\pi = 20$ MHz as in the main text and (b) $\kappa/2\pi = 80$ MHz. Each line represents a unique simulation with a different drive amplitude $\mathcal{E}$ starting in the ground state. The frequency of the drive is set to the mean of the pulled resonator frequencies associated with the ground and excited states of the transmon, as is often the case in dispersive readout. For the larger $\kappa$, the system sweeps through the resonance $N_r \sim 30$ at a faster rate. The transition becomes more diabatic and the leakage through other states is reduced, as expected from the Landau-Zener theory.</figcaption></figure>
parameters as in the main text, except for a slightly dif-
ferent drive frequency corresponding to the mean of the
pulled resonator frequencies associated with the ground
and excited states of the transmon. This different drive
frequency only leads to a small quantitative change with
respect to Fig. 3. More importantly, Fig. 8(b) is obtained
with that same drive frequency but now a larger decay
rate $\kappa/2\pi = 80$ MHz. Here, the resonance around $N_r \sim 30$
mentioned in Sec. III A leads to smaller leakage. This is
compatible with Landau-Zener theory, since with a larger
$\kappa$ the photon number rises more rapidly, corresponding to
a faster sweep through this resonance and therefore a more
diabatic transition.

In addition, we note that the numerically observed slow
increase of the transmon population in Fig. 8(b) is over-
estimated by the semiclassical model. This is not the
case for Fig. 3(a), where the drive frequency is set to
the bare-resonator frequency. We suspect that the follow-
ing mechanism is responsible for this discrepancy: for the
drive frequency used in Fig. 8, we expect the cavity state
to evolve to two very distinct coherent states depending
on whether one initializes the transmon in the ground or
the excited states. As has been highlighted in Leroux et
al. [35], the distance between the resulting polaronic states
is likely to diminish mixing of the states. This effect is,
however, not accounted for by the semiclassical model,
since the branches {$i$} are derived from the static spectrum,
where the drive frequency is not involved.

[1] A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Approaching Unit Visibility for Control of a Superconducting Qubit with Dispersive Readout, Phys. Rev. Lett. 95, 060501 (2005).

[2] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).

[3] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Charge-insensitive qubit design derived from the Cooper pair box, Phys. Rev. A 76, 042319 (2007).

[4] T. Walter, P. Kurpiers, S. Gasparinetti, P. Magnard, A. Potočnik, Y. Salathé, M. Pechal, M. Mondal, M. Oppliger, C. Eichler, and A. Wallraff, Rapid High-Fidelity Single-Shot Dispersive Readout of Superconducting Qubits, Phys. Rev. Appl. 7, 054020 (2017).

[5] Z. K. Minev, S. O. Mundhada, S. Shankar, P. Reinhold, R. Gutiérrez-Jáuregui, R. J. Schoelkopf, M. Mirrahimi, H. J. Carmichael, and M. H. Devoret, To catch and reverse a quantum jump mid-flight, Nature 570, 200 (2019).

[6] M. Boissonneault, J. M. Gambetta, and A. Blais, Dispersive regime of circuit QED: Photon-dependent qubit dephasing and relaxation rates, Phys. Rev. A 79, 013819 (2009).

[7] M. Malekakhlagh, A. Petrescu, and H. E. Türeci, Lifetime renormalization of weakly anharmonic superconducting

qubits. I. Role of number nonconserving terms, Phys. Rev. B 101, 134509 (2020).

[8] A. Petrescu, M. Malekakhlagh, and H. E. Türeci, Lifetime renormalization of driven weakly anharmonic superconducting qubits. II. The readout problem, Phys. Rev. B 101, 134510 (2020).

[9] L. Verney, R. Lescanne, M. H. Devoret, Z. Leghtas, and M. Mirrahimi, Structural Instability of Driven Josephson Circuits Prevented by an Inductive Shunt, Phys. Rev. Appl. 11, 024003 (2019).

[10] R. Lescanne, L. Verney, Q. Ficheux, M. H. Devoret, B. Huard, M. Mirrahimi, and Z. Leghtas, Escape of a Driven Quantum Josephson Circuit into Unconfined States, Phys. Rev. Appl. 11, 014030 (2019).

[11] M. D. Reed, L. DiCarlo, B. R. Johnson, L. Sun, D. I. Schuster, L. Frunzio, and R. J. Schoelkopf, High-Fidelity Readout in Circuit Quantum Electrodynamics Using the Jaynes-Cummings Nonlinearity, Phys. Rev. Lett. 105, 173601 (2010).

[12] L. S. Bishop, E. Ginossar, and S. M. Girvin, Response of the Strongly Driven Jaynes-Cummings Oscillator, Phys. Rev. Lett. 105, 100505 (2010).

[13] M. Boissonneault, J. M. Gambetta, and A. Blais, Improved Superconducting Qubit Readout by Qubit-Induced Nonlineabilities, Phys. Rev. Lett. 105, 100504 (2010).

[14] D. Sank et al., Measurement-Induced State Transitions in a Superconducting Qubit: Beyond the Rotating Wave Approximation, Phys. Rev. Lett. 117, 190503 (2016).

[15] N. Jouppi, D. Yoon, G. Kurian, S. Li, N. Patil, J. Laudon, C. Young, and D. Patterson, A domain-specific supercomputer for training deep neural networks, Commun. ACM 63, 67 (2020).

[16] T. Lu, T. Marin, Y. Zhuo, Y.-F. Chen, and C. Ma, in 2021 IEEE 18th International Symposium on Biomedical Imaging (ISBI) (IEEE, Nice, France, 2021), p. 783.

[17] T. Lu, Y.-F. Chen, B. Hechtman, T. Wang, and J. Anderson, Large-scale discrete Fourier transform on TPUs, IEEE Access 9, 93422 (2021).

[18] F. Huot, Y.-F. Chen, R. Clapp, C. Boneti, and J. Anderson, High-resolution imaging on TPUs, arXiv:1912.08063 [cs.CE] (2019).

[19] F. Belletti, D. King, K. Yang, R. Nelet, Y. Shafi, Y.-F. Shen, and J. Anderson, in Proceedings of the 2020 SIAM Conference on Parallel Processing for Scientific Computing (PP) (SIAM, Seattle, Washington, USA, 2020), p. 12.

[20] T. Lu, T. Marin, Y. Zhuo, Y.-F. Chen, and C. Ma, in 2020 IEEE High Performance Extreme Computing Conference (HPEC) (IEEE, Waltham, MA, USA, 2020), p. 1.

[21] Q. Wang, M. Ihme, Y.-F. Chen, and J. Anderson, A TensorFlow simulation framework for scientific computing of fluid flows on tensor processing units, Comput. Phys. Commun. 274, 108292 (2022).

[22] A. G. M. Lewis, J. Beall, M. Ganahl, M. Hauru, S. B. Mallick, and G. Vidal, Large-scale distributed linear algebra with tensor processing units, Proc. Natl. Acad. Sci. 119, e2122762119 (2022).

[23] R. Pederson, J. Kozlowski, R. Song, J. Beall, M. Ganahl, M. Hauru, A. G. M. Lewis, S. B. Mallick, V. Blum, and G. Vidal, Tensor processing units as quantum chemistry supercomputers, arXiv:2202.01255 [physics] (2022).
[24] A. Morningstar, M. Hauru, J. Beall, M. Ganahl, A. G. M. Lewis, V. Khemani, and G. Vidal, Simulation of Quantum Many-Body Dynamics with Tensor Processing Units: Floquet Prethermalization, PRX Quantum **3**, 020331 (2022).

[25] M. Hauru, A. Morningstar, J. Beall, M. Ganahl, A. Lewis, and G. Vidal, Simulation of quantum physics with tensor processing units: Brute-force computation of ground states and time evolution, arXiv:2111.10466 [quant-ph] (2021).

[26] E. Gustafson, B. Holzman, J. Kowalkowski, H. Lamm, A. C. Y. Li, G. Perdue, S. V. Isakov, O. Martin, R. Thompson, J. Beall, M. Ganahl, G. Vidal, and E. Peters, in *2021 IEEE/ACM Second International Workshop on Quantum Computing Software (QCS)* (2021), p. 72.

[27] M. Ganahl et al., Tensor processing units for simulating quantum circuits (to be published).

[28] M. Ganahl, J. Beall, M. Hauru, A. G. M. Lewis, J. H. Yoo, Y. Zou, and G. Vidal, Density matrix renormalization group with tensor processing units, arXiv:2204.05693 [cond-mat.str-el] (2022).

[29] J. Cohen, A. Petrescu, R. Shillito, and A. Blais, Reminiscence of classical chaos in driven transmons, arXiv:2207.09361 [quant-ph] (2022).

[30] S. Boutin, D. M. Toyli, A. V. Venkatramani, A. W. Eddins, I. Siddiqi, and A. Blais, Effect of Higher-Order Nonlinearities on Amplification and Squeezing in Josephson Parametric Amplifiers, Phys. Rev. Appl. **8**, 054030 (2017).

[31] V. V. Sivak, N. E. Frattini, V. R. Joshi, A. Lingenfelter, S. Shankar, and M. H. Devoret, Kerr-Free Three-Wave Mixing in Superconducting Quantum Circuits, Phys. Rev. Appl. **11**, 054060 (2019).

[32] E. A. Sete, J. M. Gambetta, and A. N. Korotkov, Purcell effect with microwave drive: Suppression of qubit relaxation rate, Phys. Rev. B **89**, 104516 (2014).

[33] M. Malekakhlagh, W. Shanks, and H. Paik, Optimization of the resonator-induced phase gate for superconducting qubits, Phys. Rev. A **105**, 022607 (2022).

[34] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: Composable transformations of PYTHON+NumPy programs, (2018).

[35] C. Leroux, A. Di Paolo, and A. Blais, Superconducting Coupler with Exponentially Large On:Off Ratio, Phys. Rev. Appl. **16**, 064062 (2021).