Dynamics of Transmon Ionization

Qubit measurement and control in circuit 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 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 which may explain the experimentally observed degradation in measurement fidelity.

Qubit measurement and control in circuit 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 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 which may explain the experimentally observed degradation in measurement fidelity.

I. INTRODUCTION
Dispersive readout in circuit 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 which is expected to lead to fast, high-fidelity and quantum non-demolition (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][7][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 ionization 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 is consequently modified.
Accurately simulating the dynamics of transmon ionization requires describing the system's density matrix on a truncated Hilbert space of very large dimension, something that is made possible here by the use of largescale computational accelerators known as Tensor Processing Units (TPUs), introduced below. While other studies were 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 measurement. Accounting for tens of transmon levels and hundreds of resonator states, we moreover see signatures of the high-power readout [11][12][13]. We interpret these numerical results using a semiclassical model capturing the nonlinear 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 which are in qualitative agreement with experiments [4]. This paper is organized as follows. In Sec. II 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 which 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 illustrating 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 E and frequency ω d on the resonator, the Hamil- (1) The first two terms correspond to the free transmon Hamiltonian with charging energy E c , Josephson energy E J , charge operatorn t and phase operatorφ t . We denote the eigenenergies and eigenstates of the free transmon Hamiltonian E i and |i , respectively. The first two of those eigenstates, labelled {|0 , |1 }, span the computational basis of the qubit. Of the higher excited states, a number ∼ 2E J /ω p are bound states lying within the cosine potential well illustrated in Fig. 1(b). Here, ω p = √ 8E C E J is the plasma frequency which is approximately the transmon's 0-1 transition frequency [3]. Moreover, we label |n the eigenstates of the free resonator Hamiltonian of frequency ω r , corresponding to the bosonic annihilation operatorâ. The transmonresonator coupling of amplitude g in the second line of Eq. (1) includes fast-rotating terms that 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 i,n and i, n .
Including cavity loss at a rate κ, the driven transmonresonator system is described by the usual Lindblad master equation [ with the dissipator Because we are interested in capturing the dynamics of the system in the presence of a large amplitude 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 1, 024 states in the resonator, corresponding to a total Hilbert space dimension 2 15 = 32,768. That is, the joint density matrixρ for the transmon-resonator system is a Hermitian matrix of size 2 15 × 2 15 , which thus contains 2 30 = 1,073,741,824 timedependent 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. More details are provided in Appendix A.
For a given drive amplitude E, simulations begin with the system initialized in either the dressed ground 0, 0 or excited 1, 0 states. To reduce the simulation time while capturing the transient dynamics including the ionization of the transmon, we simulate the evolution over a time at least κ −1 . After each period of the drive, the reduced transmon,ρ t , and resonator,ρ r , density matrices are recorded.
Except otherwise stated, in all of our simulation as well as in the semiclassical model discussed in Sec. III we use the parameters E J /E C = 50, E C /h = 280 MHz corresponding to approximately 6 bound transmon states within the cosine potential, see Fig. 1. The resonator frequency is ω r /2π = 7.5 GHz and g/2π = 250 MHz. A resonator loss rate of κ/2π = 20 MHz is chosen to ensure fast resonator dynamics and the drive amplitude takes values in the range E/2π ∈ [0, 450] MHz. We fix the drive frequency to the bare resonator frequency, ω d = ω r . The above parameters place the system in the dispersive regime with χ/κ = 0. 28 [29] when the resonator population is much smaller than the critical photon number Fig. 2(a,b) illustrate the average photon number N r = â †â (red lines) and the average transmon population N t = i i i|ρ t |i (black line) as a function of time for the drive amplitude E/2π = 280 MHz. The system is initialized in the dressed ground state 0, 0 in panel (a) and in the dressed excited state 1, 0 in panel (b). Also shown is the instantaneous distribution N i,t = i|ρ t |i 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 amplitude used here, the simulated measurement is far from QND and the dynamics therefore not well described by a dispersive Hamiltonian. As illustrated in Fig. 2(c,d), transmon ionization is accompanied by a sudden drop in the purity P = Tr[ρ 2 t ] 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 panel (c) and N r = 42 for the excited state in panel (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 also showed drops in transmon purity in steady-state numerical calculations. Resonances at large photon numbers have also been observed in Ref. [14].

B. Numerical Results
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 Fig. 2(a,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 exhibit 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 of purity observed in Fig. 2(c,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 E/2π = 0.28 GHz as used in Fig. 2(a,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 |α 1 associated to the transmon's first excited state, see panel (a). As the transmon population in other levels increases, additional features emerge corresponding to the coherent states associated to the now occupied transmon levels. Additionally, "bananization" caused by transmon-induced nonlinearities become apparent [30,31]. The dashed lines correspond to fixed photon numbers and are used as guide to the eye. The full colored lines overlaying the Wigner functions are obtained from a semiclassical approximation which we now introduce.

III. STATE IDENTIFICATION AND SEMICLASSICAL INTERPRETATION
The TPU-based, large-scale numerical simulations presented in the previous section clearly illustrate the breakdown of the dispersive approximation. Notably, we observe at specific resonator photon populations 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.

A. Identification of dressed states and resonator branches
Our semiclassical approach is based on the dressed states i, n and energies E i,n of the transmon-resonator Hamiltonian Eq. (1) in the absence of the drive. More precisely, at arbitrary photon number we identify the dressed states which 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 , the bare transmon states at zero photon population. Then, given an identified state i, n for n ≥ 0, the next state i, n + 1 is chosen from the subset of remaining eigenstates {|λ } such that the overlap with the state i, n plus one additional resonator excitation is maximum. 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 photon number. The full coloured lines are branches which 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 non-zero 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 ∼ 110, and between {1} and {5} for N r ∼ 50, in agreement with the photon numbers at which sudden drops of transmon purity is observed in Fig. 2(c,d). 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 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 Fig. 3(a,b). In contrast to what is observed at N r ∼ 50 and 110, this sharper resonance results in very little population transfer. This can be understood from the framework of Landau-Zener transitions. In our semiclassical picture, as the resonator rings up towards its transmon-state-dependent steady state, the photon number is swept at a rate related to κ and to the drive amplitude. A large drive amplitude leads to a rapid sweep through the narrow feature at N r ∼ 30 resulting in a diabatic passage with little leakage, see the darker blue numerical lines in Fig. 3(b). Imperfect diabatic transitions at N r ∼ 30 can be observed for smaller drive amplitudes. On the other hand, for our choice of parameters, the features at N r ∼ 50 and 110 are broader and the sweep rate therefore comparatively slower. This results in a non-diabatic 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 where numerical experiments with a different drive frequency and loss rates are presented.

B. Semiclassical dynamics
Following Ref. [13], to each branch {i} we assign an effective oscillator of photon-number dependent frequency As seen in Fig. 5 where we plot ω 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 obeys the classical equation of motion of a driven damped oscillatorα where nonlinear effects are encapsulated in the dependence of the branch frequency on photon number ω i (|α i | 2 ) according to Eq. (5). Because the quantity |α i | 2 takes arbitrary real values, we generalize Eq. (5) from discrete values n to a continuous function. To do this, we smooth ω 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 ω i (n) caused by strong interactions with other branches. This semiclassical approximation is expected to accurately describe the system for large photon number |α i | 2 1. Numerically integrating Eq. (6), we plot in Fig. 4(a) the time dependence of α 1 (t) (red lines) together with a snapshot of the Wigner function at time t = 0.13κ −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 α s 1 . To account for transitions between transmon states, the evolution of α i (t) associated to 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 Fig. 4(a-d) correspond to branch {0} which appears as a result of Purcell decay. In the same way, starting in panel (b) we include the flow of α 5 (t) associated to branch {5} (brown lines) which has a strong overlap with {1} when N r ∼ 50. Following the evolution of the Wigner function from panel to panel, it is possible to see a feature following the flow of α 5 (t) and settling at the expected steady-state value α s 5 in the last panel (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 behaviour can be understood. In particular, starting in panel (d) we see the bistable behavior of {α 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 ∆N r obtained by solving the semiclassical expression Eq. (6) given that the transmon is initialized in its ground or excited state. As expected from   5, at some threshold power the photon number rapidly increase 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 ∆N r . Panel (b) shows the same quantity obtained from numerical simulations. Given the large photon numbers that are involved, those simulations are particularly demanding and required 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 ∆N r goes to zero at the strongest drive amplitudes. This is expected from Fig. 5 where at very large photon numbers the frequencies ω i (|α i | 2 ) eventually collapse to the bare value ω r . The different drive amplitude at which this collapse occurs depend on the initial transmon state and the resulting large ∆N r is exploited in the high-power qubit readout [11][12][13].
A possible interpretation for the observed large response at the bare resonator frequency is that once the transmon is ionized, mostly charge-like states are occupied. Because these states couple longitudinally to the resonator, they do not lead to a resonator frequency pull [10]. However, 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 unbound states. We leave a more detailed analysis of this effect to future work [33].

C. Resonances at low photon number
For the parameters used above, the first resonance leading to non-QNDness occurs at a relatively large photon numbers of ∼ 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 experi-  [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 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 ob-tained from our semiclassical model. At N r ≈ 5, the transmon population rapidly rises for branch {1} associated with the transmon's first excited state (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 √ n fluctuations of coherent states, the resonance at N r ≈ 5 already plays a role at N r ≈ 2.5 leading to the observed non-QNDness. While care must be taken when using a semiclassical model at such small photon population, these observations are suggestive 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. CONCLUSION
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 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 dispersive readout which is experimentally observed at small drive amplitudes could be due these resonances. The semiclassical model developed here can be used as a tool to avoid these spurious effects when optimizing device parameters for readout.  The frequency of the drive is set to mean of the pulled resonator frequencies associated to the ground and excited states of the transmon, as is often the case in dispersive readout. For the larger κ, the system sweeps through the resonance Nr ∼ 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.
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 parameters as in the main text except for a slightly different drive frequency corresponding to the mean of the pulled resonator frequencies associated to the transmon's ground and excited states. 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 rates κ/2π = 80 MHz. The resonance around N r ∼ 30 mentioned in Sec. III A leads here to smaller leakage. This is compatible with Landau-Zener theory since with a larger κ the photon number rises more rapidly corresponding to a faster sweeps 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 overestimated 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 following 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 whether one initializes the transmon in the ground or the excited states. As it was highlighted in Leroux et al. [35], the distance between the resulting polaronic states are 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.