# Unified picture of measurement-induced ionization in the transmon

Marie Frédérique Dumas  Institut Quantique and Département de Physique,
Université de Sherbrooke, Sherbrooke J1K 2R1 QC, Canada  Benjamin Groleau-Paré
Institut Quantique and Département de Physique, Université de Sherbrooke,
Sherbrooke J1K 2R1 QC, Canada  Alexander McDonald  Institut Quantique and
Département de Physique, Université de Sherbrooke, Sherbrooke J1K 2R1 QC,
Canada  Manuel H. Muñoz-Arias  Institut Quantique and Département de Physique,
Université de Sherbrooke, Sherbrooke J1K 2R1 QC, Canada  Cristóbal Lledó
Institut Quantique and Département de Physique, Université de Sherbrooke,
Sherbrooke J1K 2R1 QC, Canada  Benjamin D’Anjou  Institut Quantique and
Département de Physique, Université de Sherbrooke, Sherbrooke J1K 2R1 QC,
Canada  Alexandre Blais  Institut Quantique and Département de Physique,
Université de Sherbrooke, Sherbrooke J1K 2R1 QC, Canada  Canadian Institute
for Advanced Research, Toronto, M5G 1M1 Ontario, Canada

###### Abstract

Despite the high measurement fidelity that can now be reached, the dispersive
qubit readout of circuit quantum electrodynamics is plagued by a loss of its
quantum nondemolition character and a decrease in fidelity with increased
measurement strength. In this work we elucidate the nature of this dynamical
process, which we refer to as transmon ionization. We develop a comprehensive
framework which provides a unified physical picture of the origin of transmon
ionization. This framework consists of three complementary levels of
descriptions: a fully quantized transmon-resonator model, a semiclassical
model where the resonator is treated as a classical drive on the transmon, and
a fully classical model. Crucially, all three approaches preserve the full
cosine potential of the transmon, and lead to similar predictions. This
framework identifies the multiphoton resonances responsible for transmon
ionization. It also allows us to efficiently compute numerical estimates of
the photon number threshold for ionization, which are in remarkable agreement
with recent experimental results. The set of tools developed within this work
are both conceptually and computationally simple, and we expect them to become
an integral part of the theoretical support of all circuit QED experiments.

##  I Introduction

Circuit quantum electrodynamics with transmon qubits is a leading platform for
quantum computation Blais  _et al._ (2021); Koch  _et al._ (2007). A key
feature of this architecture is the dispersive readout where the qubit state
is inferred from quadrature measurement of a microwave signal transmitted or
reflected from a resonator coupled to the qubit Wallraff  _et al._ (2005). In
principle, increasing the amplitude of the measurement drive, or equivalently
the average resonator photon population, should lead to an increase in the
measurement’s signal-to-noise ratio (SNR). Moreover, because the qubit is
strongly detuned from the resonator and thus from the readout drive, the
process is expected to be quantum nondemolition (QND) Blais  _et al._ (2004).

In practice, however, increasing the readout drive amplitude results in a
reduction of the measurement fidelity and to measurement-induced transitions
spoiling the QND character of the readout. This is a widely observed
phenomenon Walter  _et al._ (2017); Sank  _et al._ (2016); Khezri  _et al._
(2023); Minev  _et al._ (2019); Swiadek  _et al._ (2023); Johnson  _et al._
(2012). Although it has been experimentally possible to optimize the
dispersive readout to obtain high-fidelity state discrimination Walter  _et
al._ (2017); Touzard  _et al._ (2019); Sunada  _et al._ (2022); Lledó  _et
al._ (2023); Chen  _et al._ (2023); Swiadek  _et al._ (2023), single-shot
qubit readout now lags in fidelity behind the best one- and two-qubit gates
Kjaergaard  _et al._ (2020). Improving qubit readout is thus an outstanding
challenge for the field. Attempts at understanding the origin of these
observations have been made using perturbative expansions in the qubit-
resonator coupling, qubit anharmonicity or drive amplitudes, while treating
the qubit as a two-level system or a Kerr nonlinear oscillator Boissonneault
_et al._ (2008, 2009); Petrescu  _et al._ (2020); Hanai  _et al._ (2021).
While these approaches highlight important mechanisms for qubit decay and
excitation in circuit QED, they do not explain the observed non-QNDness of the
dispersive readout. This situation is exacerbated by the fact that no clear
parameter dependence for the onset of non-QNDness emerges from experimental
observations.

There has recently been a flurry of activity aimed at understanding the origin
of the failure of the dispersive readout Shillito  _et al._ (2022); Cohen  _et
al._ (2023); Khezri  _et al._ (2023); Sank  _et al._ (2016); Xiao  _et al._
(2023). In particular, numerical simulations of the dynamics of the dispersive
readout have shown that the measurement drive, although strongly detuned from
the qubit in the dispersive regime, can cause measurement-induced transitions
to highly excited states of the transmon at specific resonator photon number
Shillito  _et al._ (2022). In some cases, this leakage extends to states
beyond the transmon’s cosine potential well, a phenomenon referred to as
transmon ionization. Cohen  _et al._ (2023) moreover showed that the chaotic
behavior of the driven nonlinear classical pendulum can be used to gain
qualitative understanding of transmon ionization, an observation which
suggests that perturbative methods may not be sufficient to understand
ionization. In view of the transmon’s negative anharmonicity, the importance
of the qubit-resonator frequency detuning, i.e., placing the transmon
frequency above or below that of the resonator, was highlighted in Refs.
Khezri  _et al._ (2023); Cohen  _et al._ (2023).

However, it is unclear how these seemingly disparate approaches relate to each
other and whether they make the same quantitative predictions. In this work,
we present a unified picture of ionization in the transmon qubit. To do so, we
divide the analysis into three main levels of approximation, in all cases
considering both positive and negative detuning. We first consider a fully
quantized transmon-resonator model which enables us to investigate the
mechanisms leading to ionization. We then argue that under typical conditions,
a simplified semiclassical model where the resonator is treated as a classical
field within Floquet theory is able to accurately predict ionization. Finally,
we demote the transmon to a fully classical system and use the Bohr-Sommerfeld
quantization rule to make predictions about ionization. Remarkably, we show
that all three models lead to similar predictions for the onset of ionization,
and that they agree with experimental results. We also highlight the
importance of gate charge in the onset of ionization of the computational
states, even deep in the transmon regime Cohen  _et al._ (2023); Khezri  _et
al._ (2023). Common to all three approaches is that they consider the full
cosine potential of the transmon rather than the common Kerr nonlinearity
approximation.

An important feature emerging from all three methods is that there in general
exist several critical photon numbers at which ionization can in principle
occur. These critical points correspond to resonances that are activated by
the strong resonator field due to the Josephson junction’s nonlinearity. We
show that not all these resonances play an equally important role and that
some can be safely ignored in the appropriate conditions. Moreover, we find
that the critical photon numbers relevant to ionization are not simply related
to the critical photon number $n_{\rm crit}$ of the Jaynes-Cummings
Hamiltonian at which the dispersive approximation breaks down Blais  _et al._
(2004). Indeed, while the latter indicates when the qubit-induced nonlinearity
becomes important, it is not a good predictor for the onset of ionization.
Throughout this work we use the term “ionization” loosely, that is, even if in
some instances the measurement-induced transitions do not involve states above
the cosine potential well.

##  II Fully quantum model

Figure 1: (a) Schematic of a transmon qubit capacitively coupled to a driven
resonator. (b) Cosine potential well of the transmon with $E_{J}/E_{C}=110$
and $n_{g}=0$ together with the first 13 eigenstates represented in the phase
basis. There are about 9 levels in the well. Transmon states at the bottom of
the well are close to eigenstates of the harmonic oscillator while states
above the well are close to charge states, i.e., plane waves in the phase
basis.

###  II.1 Model

As illustrated in Fig. 1(a) we consider a transmon qubit capacitively coupled
to a resonator described by the Hamiltonian ($\hbar=1$) Koch  _et al._ (2007);
Blais  _et al._ (2021)

$\begin{split}\hat{H}_{tr}&=\omega_{r}\hat{a}^{\dagger}\hat{a}+\hat{H}_{t}-ig(\hat{n}_{t}-n_{g})(\hat{a}-\hat{a}^{\dagger}),\end{split}$
(1)

where $\omega_{r}$ is the bare resonator frequency, $g$ the transmon-resonator
coupling strength, and $\hat{a}$ the annihilation operator of the resonator.
Moreover,
$\hat{H}_{t}=4E_{C}(\hat{n}_{t}-n_{g})^{2}-E_{J}\cos\hat{\varphi}_{t}$ is the
transmon Hamiltonian with $E_{C}$ the charging energy, $E_{J}$ the Josephson
energy, and $n_{g}$ the gate charge. The operators $\hat{\varphi}_{t}$ and
$\hat{n}_{t}$ are the transmon’s canonically-conjugate phase and charge
operators, respectively. Throughout this work, we assume that the qubit is
operated in the transmon regime with $E_{J}/E_{C}\gg 1$. Note that we include
the often-dropped contribution of the gate charge to the last term of
$\hat{H}_{tr}$ to preserve the correct gauge freedom of the Hamiltonian.

The resonator drive used to measure the qubit takes the form

$\displaystyle\hat{H}_{d}(t)=-i\varepsilon_{d}\sin(\omega_{d}t)\left(\hat{a}-\hat{a}^{\dagger}\right),$
(2)

with $\varepsilon_{d}$ the drive amplitude and $\omega_{d}$ the drive
frequency. Including single-photon resonator losses, the master equation for
the transmon-resonator density matrix is then Blais  _et al._ (2021)

$\displaystyle\partial_{t}\hat{\rho}=-i[\hat{H}_{tr}+\hat{H}_{d}(t),\hat{\rho}]+\kappa\mathcal{D}[\hat{a}]\hat{\rho},$
(3)

where $\kappa$ is the resonator decay rate and
$\mathcal{D}[\hat{a}]\hat{\rho}=\hat{a}\hat{\rho}\hat{a}^{\dagger}-\\{\hat{a}^{\dagger}\hat{a},\hat{\rho}\\}/2$
the Lindblad dissipator.

In the absence of drive and qubit-resonator coupling, we denote the
eigenstates of the Hamiltonian by $|i_{t},n_{r}\rangle$ where $i_{t}$ labels
the bare transmon eigenstates and $n_{r}$ labels the bare resonator photon
number. The bare energy associated to $|i_{t},n_{r}\rangle$ is
$E_{i_{t}}+n_{r}\omega_{r}$, and the transition frequencies between different
states of the transmon are denoted $\omega_{i_{t}j_{t}}\equiv
E_{j_{t}}-E_{i_{t}}$. With this notation, the bare qubit transition frequency
is $\omega_{q}=E_{1_{t}}-E_{0_{t}}$. The matrix elements of the qubit-
resonator coupling are $g_{i_{t}j_{t}}=g\langle
i_{t}|\hat{n}_{t}|j_{t}\rangle$.

As is usual for quantum information processing in circuit QED, throughout this
work we assume that the system is operated in the dispersive limit where
$|g/\Delta|\ll 1$ with $\Delta=\omega_{q}-\omega_{r}$ the qubit resonator
detuning. In that regime, the computational subspace of the qubit is spanned
by the dressed ground $|\overline{0_{t},0_{r}}\rangle$ and excited
$|\overline{1_{t},0_{r}}\rangle$ states of $\hat{H}_{0}$. Logical operations
act on these two dressed states, and it is those states that are resolved in a
dispersive qubit measurement. A more precise condition for the dispersive
regime to hold is that $\bar{n}_{r}/n_{{\rm crit},i_{t}}\ll 1$, where
$\bar{n}_{r}$ is the average resonator photon number and $n_{{\rm
crit},i_{t}}$ is a photon number at which the transmon-resonator interaction
stops being dispersive for state $i_{t}$ (see Appendix A). With the above
notation, the commonly used Jaynes-Cummings critical photon number is $n_{{\rm
crit},0_{t}}=|\Delta/2g_{0_{t}1_{t}}|^{2}$ Blais  _et al._ (2004, 2021). As
will become clear below, the breakdown of the dispersive approximation during
a qubit readout does not necessarily mean that there will be ionization.

Unless otherwise stated, throughout this work we fix the transmon parameters
to be $E_{J}/E_{C}=110$ with
$E_{C}/2\pi=$220\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$. This corresponds to
a qubit transition frequency
$\omega_{q}/2\pi=$6.298\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$,
anharmonicity of magnitude
$\alpha/2\pi=$240\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$, and from 9 to 10
levels in the cosine potential well depending on the value of $n_{g}$, see
Fig. 1(b). Moreover, we take the coupling to be
$g/2\pi=$120\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$. Due to the negative
sign of the transmon’s anharmonicity, there is a striking difference in the
onset of ionization depending on the sign of the qubit-resonator detuning
$\Delta$ Khezri  _et al._ (2023); Cohen  _et al._ (2023). In this work, we
study both negative and positive detunings. When considering negative
detunings, the resonator frequency is
$\omega_{r}/2\pi=$7.5\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$
($\Delta/2\pi=$-1.202\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$, $n_{{\rm
crit},0_{t}}=14$, $n_{{\rm crit},1_{t}}=10.5$) while for positive detuning we
take $\omega_{r}/2\pi=$5.3\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$
($\Delta/2\pi=$0.998\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$, $n_{{\rm
crit},0_{t}}=9.7$, $n_{{\rm crit},1_{t}}=2.9$). In Sects. III and IV the
driven resonator is substituted by a direct drive on the transmon. There, we
take $\omega_{d}/2\pi=$7.515\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ at
negative detuning and
$\omega_{d}/2\pi=$5.267\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ at positive
detuning. In both cases, the chosen drive frequency is positioned between the
two pulled resonator frequencies corresponding to the qubit’s ground and
excited states.

Crucially, although they are discussed for those specific sets of parameters,
the phenomena we discuss in this article are generic. Once the sign of the
detuning $\Delta$ has been fixed, the phenomenology for a given choice of
parameters is easily understood.

###  II.2 Branch Analysis

To understand how the presence of photons in the resonator leads to transmon
ionization, we first label dressed states $|\overline{i_{t},n_{r}}\rangle$ by
considering which is the closest to the bare transmon state $\ket{i_{t}}$ at
an arbitrary photon number $\ket{n_{r}}$. Each such sets of states is referred
to as a branch. While this identification can easily be done at low photon
number, there is no clear such closest state when the photon number approaches
or exceeds $n_{\rm crit}$, where the dressed states are highly entangled
qubit-resonator states. To build branches that allow us to gain insight into
the dynamics of the driven system at high photon numbers, we instead follow
Refs. Shillito  _et al._ (2022); Boissonneault  _et al._ (2010). In this
approach, we use as our starting point the eigenstates $\\{|\lambda\rangle\\}$
of $\hat{H}_{tr}$ obtained from numerical diagonalization. For each $i_{t}$,
we identify $|\overline{i_{t},0_{r}}\rangle$ to be the low-energy eigenstate
with the largest overlap with $|i_{t},0_{r}\rangle$. The branches are built
recursively from those starting points and in parallel: given
$|\overline{i_{t},n_{r}}\rangle$, the next eigenstate
$|\overline{i_{t},n_{r}+1}\rangle$ added to each branch is the one which
maximizes the overlap

$\displaystyle
C_{\overline{i_{t},n_{r}}}(\lambda)\equiv|\langle\lambda|\hat{a}^{\dagger}|\overline{i_{t},n_{r}}\rangle|^{2}.$
(4)

At each step, $|\lambda\rangle$ is taken from the set of unassigned states. In
this way, we obtain branches labeled $B_{i_{t}}$ corresponding to the set of
eigenstates $\\{\ket{\overline{i_{t},n_{r}}}\\}$ with fixed $i_{t}$ and
$n_{r}$ spanning the full resonator Hilbert space. When the dispersive
approximation is valid, each such branch can thus be pictured as an effective
oscillator responding at the pulled resonator frequency associated with the
transmon state $i_{t}$ Shillito  _et al._ (2022).

The character of the branch eigenstates can be investigated by tracking their
average transmon population $N_{t}\equiv\sum_{i_{t},n_{r}}i_{t}|\langle
i_{t}|\overline{i_{t},n_{r}}\rangle|^{2}$ and their average resonator
population
$N_{r}\equiv\langle\overline{i_{t},n_{r}}|\hat{a}^{\dagger}\hat{a}|\overline{i_{t},n_{r}}\rangle$.
At very low photon numbers, the dispersive approximation holds and we expect
the states within a given branch $B_{i_{t}}$ to have a transmon population
$N_{t}$ close to the bare value $i_{t}$. However, plotting the average
transmon and resonator populations parametrically for each branch reveals that
this is not always the case, see Fig. 2(a). Instead, a drastic change in the
average transmon population $N_{t}$ is sometimes observed at specific branch-
dependent values of $N_{r}$, such as a “swapping” of the branches in which two
branches exchange their values of $N_{t}$. These features are a signature of
strong hybridization between transmon states and indicate a complete breakdown
of the dispersive approximation. When swapping occurs between a branch linked
to computational states and a branch that is near the top of the cosine
potential, there is a significant change in the character of the former. It
shifts from being predominantly made up of low-energy, localized bare transmon
states to being composed largely of high-energy, charge-like bare transmon
states. We thus generically refer to this process as ionization.

In Refs. Shillito  _et al._ (2022); Shillito (2023), these drastic changes
were shown to determine dynamical properties of the system in the presence of
a readout drive populating the resonator. The branch analysis thus serves as a
diagnosis of the onset of measurement-induced ionization of the transmon.
Below, we first focus on negative transmon-resonator detuning, followed by the
case where this detuning is positive.

Figure 2:  Negative detuning. (a) Average transmon population and (b) modular
branch eigenenergies versus $N_{r}$ for a negative detuning
$\Delta/2\pi=-$1.202\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$. Branches
$B_{0_{t}}$, $B_{1_{t}}$, $B_{7_{t}}$, $B_{10_{t}}$ and $B_{11_{t}}$ are
highlighted in color (see legend), the others are gray. Branch population
swappings occur at avoided crossings in the modular energy spectrum (vertical
dashed lines). (c) Dynamics when initializing in the excited state
$\ket{\overline{1_{t},0_{r}}}$. When the resonator population reaches
$\langle\hat{a}^{\dagger}\hat{a}\rangle\approx 84$ corresponding to the
$B_{1_{t}}-B_{7_{t}}$ swapping observed in panel (a), the population of the
state is mostly transferred to branch $B_{7_{t}}$, but a significant part of
the population ($\approx 7\%$) stays in branch $B_{1_{t}}$, indicating
ionization. At $\langle\hat{a}^{\dagger}\hat{a}\rangle\approx 140$, the
population is then mostly transferred to $B_{10_{t}}$, corresponding to the
$B_{7_{t}}-B_{10_{t}}$ swapping in (a), while residual population in
$B_{7_{t}}$ is ionized. Inset in (c) shows in log scale the bimodal Fock state
distribution of the resonator at the final time $\kappa t=2$. For dynamics
simulations, we use
$\varepsilon_{d}/2\pi=$180\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$,
$\omega_{d}/2\pi=$7.515\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$, and
$\kappa/2\pi=$7.95\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ (set to match the
full $\chi$ shift at zero photons). The full-time dynamics were performed
using quantum trajectories Carmichael (2007) with 16 transmon states, 300
resonator states, 200 trajectories, and time-evolving up to $\kappa t=2$.

###  II.3 Negative detuning $\omega_{q}<\omega_{r}$

We consider the case where the qubit frequency is below the resonator
frequeny, $\Delta=\omega_{q}-\omega_{r}<0$. To proceed, we note that Eq. 4 can
be rewritten as Cohen  _et al._ (2023)

$\displaystyle
C_{\overline{i_{t},n_{r}}}(\lambda)=g^{2}\frac{|\langle\lambda|\hat{n}_{t}|\overline{i_{t},n_{r}}\rangle|^{2}}{(E_{\lambda}-E_{\overline{i_{t},n_{r}}}-\omega_{r})^{2}}.$
(5)

The denominator of this expression conveys that, as expected, when two
eigenstates are nearly degenerate dramatic effects such as branch swapping can
be expected. But when do we expect dressed states of $\hat{H}_{tr}$ to be
nearly-resonant?

Perturbation theory provides an intuitive answer to this question. Using the
usual Schrieffer-Wolff transformation $e^{-\hat{S}}$ with an appropriately
chosen $\hat{S}$, one can approximately diagonalize $\hat{H}_{tr}$ as Koch
_et al._ (2007)

$\displaystyle e^{\hat{S}}\hat{H}_{tr}e^{-\hat{S}}$
$\displaystyle\approx\sum_{i_{t}=0}^{\infty}\left[E_{i_{t}}+\Lambda_{i_{t}}+(\omega_{r}+\chi_{i_{t}})\hat{a}^{\dagger}\hat{a}\right]|i_{t}\rangle\langle
i_{t}|,$ (6)

where $\Lambda_{i_{t}}$ are the Lamb shifts and

$\displaystyle\chi_{i_{t}}=\sum_{j_{t}=0}^{\infty}|g_{i_{t}j_{t}}|^{2}\left(\frac{1}{\omega_{r}-\omega_{i_{t}j_{t}}}-\frac{1}{\omega_{r}+\omega_{i_{t}j_{t}}}\right)$
(7)

are the $\chi$ shifts obtained here without the usual rotating-wave
approximation (RWA) on the transmon-resonator coupling term of Eq. 1 Blais
_et al._ (2021). Importantly, because of the transmon’s negative
anharmonicity, at negative detuning $\Delta<0$ the transitions between
nearest-neighbor higher-excited states are more off-resonant with the
resonator, e.g., $|\Delta-\alpha|>|\Delta|$. This trend continues for all
transmon eigenstates $i_{t}$ within the cosine potential and ensures the
validity of perturbation theory at low photon number. Given the small change
$\chi_{i_{t}}$ to the dispersion of the resonator, at low photon number we
then have a direct mapping between the bare $|i_{t},n_{r}\rangle$ and dressed
states $|\overline{i_{t},n_{r}}\rangle\approx
e^{-\hat{S}}|i_{t},n_{r}\rangle\equiv|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle i_{t},n_{r}$\cr} }}}\rangle$.

To understand how this standard analysis helps explain the presence of
degeneracies, it is instructive to consider the eigenvalues of $\hat{H}_{tr}$
modulo $\omega_{r}$. This modular spectrum is plotted versus resonator photon
number in Fig. 2(b). At low photon number, the slope of the modular branch
energy versus $N_{r}$ for states in the well is $\chi_{i_{t}}$, which is
positive for states deep in the well and negative for state near the top of
the well 111Emitting and absorbing a virtual photon — corresponding to the
first and second term in Eq. 7 respectively — leads to an oppositely-signed
ac-Stark shift for the transitions of interest. The former processes are less
detuned due to the negative anharmonicity and the corresponding contributions
are larger in magnitude, leading to a positive slope $\chi_{i_{t}}$ for states
at the bottom of the cosine potential. However $|g_{i_{t}j_{t}}|^{2}$ is
larger for higher-excited states, and the virtual emission processes become
dominant, leading to a negative $\chi_{i_{t}}$ for states at the top of the
well.. At large enough photon number, these opposite ac-Stark shifts can lead
to near collisions of the modular energies between states at the bottom and at
the top of the well, corresponding to near degeneracies
$E_{\lambda}-E_{\lambda^{\prime}}\approx 0\mod\omega_{r}$ of the
eigenenergies. Nonperturbative corrections to the eigenstates—which arise due
to dropped terms in Eq. 6—turn these near collisions into avoided crossings
(dashed vertical lines), see Appendix B. As expected in this situation, the
states before and after an avoided crossing swap character. That is, if the
eigenenergies of branches $B_{i_{t}}$ and $B_{j_{t}}$ collide modulo
$\omega_{r}$, we then have
$|\overline{i_{t},n_{r}}\rangle\approx|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle i_{t},n_{r}$\cr} }}}\rangle$,
$|\overline{j_{t},n_{r}}\rangle\approx|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle j_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle j_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle j_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle j_{t},n_{r}$\cr} }}}\rangle$ before the crossing
but $|\overline{i_{t},n_{r}}\rangle\approx|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle j_{t},n_{r}+k$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle j_{t},n_{r}+k$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle j_{t},n_{r}+k$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle j_{t},n_{r}+k$\cr} }}}\rangle$,
$|\overline{j_{t},n_{r}}\rangle\approx|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle i_{t},n_{r}-k$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle i_{t},n_{r}-k$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle i_{t},n_{r}-k$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle i_{t},n_{r}-k$\cr} }}}\rangle$ after the crossing
for some integer $k$.

These qualitative predictions capture the main features observed in Fig.
2(a,b). As expected, at low photon numbers the modular energies vary linearly
with photon number, and the slope is positive for branches corresponding to
states deep in the well such as $B_{1_{t}}$. On the other hand, the slope is
negative for states near the top such as $B_{7_{t}}$. An avoided crossing is
observed when these two branches eventually meet at $N_{r}\approx 84$. This is
precisely where the branches swap character, see panel (a). The same
phenomenology holds for the ground-state branch $B_{0_{t}}$ which swaps with
the charge-like branch $B_{11_{t}}$ after an avoided crossing at a higher
photon number $N_{r}\approx 166$. An avoided crossing between $B_{7_{t}}$ and
$B_{10_{t}}$ is also observed at $N_{r}\approx 140$. Importantly, not all
near-degeneracies lead to swapping, such as when $B_{0_{t}}$ and $B_{7_{t}}$
cross at around $N_{r}\approx 23$. This is simply a consequence of the
magnitude of the relevant matrix element: at such low photon number it is not
possible to efficiently couple the bare $0_{t}$ and $7_{t}$ states of the
transmon via the qubit-resonator coupling, see Appendix B for details.

We stress that although the above leading-order perturbative approach used to
interpret Fig. 2 is qualitatively correct, to obtain quantitative predictions
one must resort to full numerics. For instance, at large enough photon numbers
the self-Kerr of the dressed resonator becomes relevant and must be taken into
account when trying to predict the positions of the avoided crossings. In
addition, ionization of the ground and excited branches can occur with charge-
like states. The modular spectrum of these highly-excited states can vary
wildly as a function of photon number, and they can thus become resonant with
the computational states. Given that $\hat{n}_{t}$ can efficiently couple
states near the top of the well to each other and even to charge-like states
outside the cosine potential Cohen  _et al._ (2023), this observation is
relevant even at low photon number. Moreover, the charge dispersion of the
transmon states near and above the top of the cosine potential well implies
that $n_{g}$ plays a crucial role in determining at which photon number
ionization occurs Cohen  _et al._ (2023); Khezri  _et al._ (2023). See Sec.
II.5 for more details on these considerations.

It is worth emphasizing that written in the language of the branch analysis,
ionization occurs when the population remains in the ground or excited state
branch after the swapping; consider e.g. the red line in Fig. 2(a) which jumps
from $N_{t}\sim 1$ to $N_{t}\sim 7$ at the anticrossing. This corresponds to
an adiabatic transition to a highly-excited transmon state. By contrast,
purely diabatic-like transitions between branches do not cause ionization. In
practice, the character of the transition can be mixed between these two types
of processes. As an illustration, by starting in
$|\overline{1_{t},0_{r}}\rangle$, time-evolving via Eq. 3, and plotting the
average population of the branches
$P_{B_{i_{t}}}(t)\equiv\langle\sum_{n_{r}}|\overline{i_{t},n_{r}}\rangle\langle\overline{i_{t},n_{r}}|\rangle$
as in Fig. 2 (c), we see a nearly full diabatic population transfer from
$B_{1_{t}}$ to $B_{7_{t}}$ at $\langle\hat{a}^{\dagger}\hat{a}\rangle\approx
84$, yet the bare transmon population (not shown) does not exceed
$N_{t}\approx 2.75$. The approximately 7% of the population that remains in
$B_{1_{t}}$ after the swapping roughly translates to an equally large
probability of ionization. Ionizing by remaining in the ground or excited
state branch after the crossing implies a change in the resonator pull, which
will in this instance be closer to $\chi_{7_{t}}$ than to $\chi_{1_{t}}$. As
discussed in Ref. Shillito  _et al._ (2022), this in turn implies that the
resonator undergoes a different phase-space trajectory, leading to the bimodal
distribution of Fock states as seen in the inset of Fig. 2(c). Note that
ionization here happens at photon number that is much larger than the Jaynes-
Cummings critical photon number.

The importance of multiphoton resonances as a mechanism causing drive-induced
transitions in the transmon has previously been discussed in the literature.
Sank  _et al._ (2016) show that resonances between different excitation-
conserving subspaces of the Jaynes-Cummings Hamiltonian mediated by several
photons can lead to transmon transitions. That model was able to explain
experimentally-observed features. Interestingly, ionization was observed to
occur at a photon number significantly larger than the Jaynes-Cummings photon
number, as is expected from the above discussion. Moreover, that experiment
also shows that while readout-induced transitions occur at specific photon
numbers, it is possible to operate the readout at still larger photon numbers.
This highlights that, as discussed above, resonances that are traversed
diabaticaly do not lead to ionization. Moreover, Xiao  _et al._ (2023) use a
diagrammatic method to compute effective Hamiltonians for a driven transmon in
the absence of a resonator. This approach allows them to identify the
multiphoton processes that are responsible for drive-induced transmon
transitions. In the next section, we show that the model of a directly driven
transmon can capture ionization due to a readout drive. The branch analysis
unifies these pictures and provide a simple way to predict where resonances
occur. It also captures resonances that cannot be obtained via perturbative
arguments. This occurs for instance when invoking the RWA for higher-energy
states of the cosine potential, which would lead one to incorrectly predict at
what photon number ionization occurs, see Sec. II.5.

We can now succinctly summarize the three key ingredients leading to
ionization of the qubit ground and excited states at negative detuning:

  * •

At low photon number, $\chi_{i_{t}}$ is positive for states deep in the well,
leading to an increase in the average energy with photon number. In contrast,
transmon states near the top of the cosine potential are pushed down in energy
with increasing photon number. Charge-like states above the well can also be
pushed down in energy, in a manner which is highly-sensitive to the gate
charge.

  * •

These opposite behaviors result in near-degeneracies and thus avoided
crossings between these two sets of states at some large-enough photon number,
leading to swapping of the transmon branches.

  * •

The resulting hybridization of the transmon’s computational state with states
at the top or above the well leads to transmon ionization at specific
resonator photon number.

###  II.4 Positive detuning $\omega_{q}>\omega_{r}$

Figure 3: Positive detuning. (a) Average transmon population, (b) modular
eigenenergies of branches and (c)-(f) projection of branch eigenstates on the
transmon eigenbasis in logarithmic scale versus $N_{r}$ for a positive
detuning $\Delta/2\pi=$0.998\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$.
Branches $B_{0_{t}}$, $B_{1_{t}}$, $B_{4_{t}}$, $B_{7_{t}}$ and $B_{8_{t}}$
are highlighted in color (see legend), the others are gray. Branch population
bunching is observed about state $i_{t}^{\star}=4$. (g) Dynamically simulated
populations of branches $1_{t}$ (red line) and $7_{t}$ (light blue line) when
initializing in the excited state $\ket{\overline{1_{t},0_{r}}}$. Dynamics
simulations are done similarly to Fig. 2(c), but with
$\varepsilon_{d}/2\pi=$93\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$,
$\omega_{d}/2\pi=$5.267\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$,
$\kappa/2\pi=$11.98\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$, and 120
resonator states.

In contrast to negative detuning where branch swapping is observed, at
moderate positive qubit-resonator detuning we instead observe that the
transmon populations of the branches within the cosine potential well coalesce
as the resonator photon number is increased, see Fig. 3(a). Much of the
phenomenology associated with this observation stems from the transmon’s
negative anharmonicity, leading to transitions between neighboring states in
the well becoming more resonant with the resonator frequency when climbing the
cosine potential. As a result, because of the small anharmonicity of the
transmon and for typical large values of $g$ in circuit QED, generically one
of the transmon transitions within the well is such that
$|E_{i_{t}+1}-E_{i_{t}}-\omega_{r}|\sim g$. For that transition, the
dispersive approximation breaks down even at zero resonator photon population.
In the Kerr oscillator approximation for the transmon where
$\omega_{i_{t},i_{t}+1}\sim\omega_{q}-i_{t}\alpha$, this 1-photon resonance
occurs for the transmon level
$i_{t}^{\star}\sim\mathrm{round}[(\omega_{q}-\omega_{r})/\alpha]$ Khezri  _et
al._ (2023).

For example, with the parameters of Fig. 3 we have that

$\begin{cases}E_{i_{t}+1}-E_{i_{t}}>\omega_{r}&\text{if $i_{t}\leq 3$},\\\
E_{i_{t}+1}-E_{i_{t}}\approx\omega_{r}&\text{if $i_{t}=4$},\\\
E_{i_{t}+1}-E_{i_{t}}<\omega_{r}&\text{if $i_{t}\geq 5$},\end{cases}$ (8)

such that $i_{t}^{\star}=4$ is the state for which the dispersive
approximation first fails. This is confirmed by examining the modular energies
of the branches, see Fig. 3(b). Unlike the prediction made by the dispersive
approximation of Eq. 6, $E_{\overline{4_{t},n_{r}}}$ does not vary linearly at
small photon number (pink line), indicating a strong mixing of the bare state
$|4_{t},n_{r}\rangle$ with other neighboring bare transmon states. Increasing
the photon number, hybridization of $i_{t}^{\star}$ with neighboring states
leads to level repulsion, and in turn to collisions with other branches.

At the moderate positive detuning considered here, this hybridization is
further accentuated by near resonances (modulo $\omega_{r}$) between pairs of
states in the well. For the parameters of Fig. 3, this is the case of the pair
of states $1_{t}$ and $7_{t}$ for which $E_{7_{t}}-E_{1_{t}}-6\omega_{r}\sim
g$, and $0_{t}$ and $8_{t}$ for which $E_{8_{t}}-E_{0_{t}}-8\omega_{r}\sim g$.
These transitions require multiple intermediate virtual transitions and are
thus suppressed at very low photon numbers, see in Fig. 3(b) the almost
resonant modular energies of branches $1_{t}$ and $7_{t}$ running parallel up
to $N_{r}\approx 7$. However, because the coupling matrix elements is
proportional to $g\sqrt{n_{r}}$, their hybridization becomes possible at
increased photon number, see Fig. 3(c-f). The early breakdown of the
dispersive approximation for $i_{t}^{\star}$ followed by the activation of
near resonances is the cause of the observed bunching of the branch population
at $N_{t}\sim i_{t}^{\star}$. This process starts with states closest to
$i_{t}^{\star}$ and, pair by pair, is followed by states connected by
multiphoton transition and that are symmetric about $i_{t}^{\star}$. This
phenomenology explains why branch bunching, and thus ionization, typically
occurs at lower photon number for the transmon’s excited than the ground
state, see Fig. 3(a). It also explains why, at moderate detuning, ionization
typically occurs at much lower photon number for positive detuning than for
negative detuning.

Figure 3(g) shows the average population of the branches when starting in
$|\overline{1_{t},0_{r}}\rangle$ and evolving using Eq. 3. Around
$\langle\hat{N}_{r}\rangle(t)\approx 7$, which coincides with branches $1_{t}$
and $7_{t}$ bunching around $N_{r}\approx 7$ in Fig. 3(a), a fraction of the
population is transferred from branch $1_{t}$ to $7_{t}$. However, even at the
relatively high speed at which the resonance is traversed, most of the
population remains in branch $1_{t}$ due to the wideness of the avoided
crossing. The resonator population drops after the crossing due to drive
becoming significantly off-resonant with the pulled resonator frequency
associated with the ionized state.

Figure 4: Positive detuning. Bare energy $\omega_{0_{t}i_{t}}-i_{t}\omega_{r}$
of the sites $\ket{i_{t},N-i_{t}}$ in the $N$-excitation subspace of
$\hat{H}_{0}$ in the RWA illustrated for $N=12$. (a) Parameters are as in Fig.
3 with $g=0$ and $n_{g}=0.25$. The site $\ket{4_{t},N-4_{t}}$ sits at the top
of the inverted parabola. The arrows represent multiphoton processes
connecting near-degenerate sites at the two extremities of the lattice. The
full vertical lines represent the charge dispersion of level $i_{t}$. The
dashed vertical line at $2E_{J}/\omega_{p}\sim 7.4$ indicates a lower bound
for the number of states in the cosine potential well. Here,
$\omega_{p}=\sqrt{8E_{C}E_{J}}$ is the plasma frequency. (b) Same parameters
except for the larger positive detuning
$\Delta/2\pi=$1.998\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$
($\omega_{r}/2\pi=$4.3\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$).

In Ref. Khezri  _et al._ (2023), this physics was shown to be captured within
a semiclassical model in the RWA. In the full quantum model we use here, this
can be understood by projecting the qubit-resonator Hamiltonian $\hat{H}_{0}$
of Eq. 1 on a subspace with $N$ excitations. After a RWA, dropping $n_{g}$
from the coupling term, and moving to a frame rotating at $\omega_{r}$ for the
qubit and the resonator, this leads to

$\begin{split}\hat{H}_{0}^{N}&=\sum_{i_{t}+n_{r}=N}(\omega_{i_{t}}-i_{t}\omega_{r})\ket{i_{t},n_{r}}\bra{i_{t},n_{r}}\\\
&-ig_{i_{t}+1,i_{t}}\sqrt{n_{r}}\ket{i_{t}+1,n_{r}}\bra{i_{t},n_{r}-1}+\mathrm{h.c.}\end{split}$
(9)

This is analogous to a 1D lattice with nearest-neighbor hopping amplitude
$g_{i_{t}+1,i_{t}}\sqrt{n_{r}}$. The bare energies of the sites
$\ket{i_{t},n_{r}}$ on this lattice form an inverted parabola with
$\ket{i_{t}^{\star},N-i_{t}^{\star}}$ at the top, see Fig. 4(a) Khezri  _et
al._ (2023). For the parameters of Fig. 3, the states $\ket{3_{t},N-3_{t}}$,
$\ket{4_{t},N-4_{t}}$, and $\ket{5_{t},N-5_{t}}$ near the top of the parabola
are near-degenerate and are connected by a single lattice hop. They therefore
rapidly hybridize leading to the almost immediate state bunching seen at low
photon number in Fig. 3(a). Hybridization between the other pairs of near
degenerate states located on each side of the lattice becomes relevant once
the photon number has increased enough to make the hopping amplitude
sufficiently large to allow virtual transitions between those two remote
states. At that point, the wavefunctions of that pair of states become
delocalized along the 1D lattice, and they join the layer of bunched states.
This occurs pair by pair until the photon number is large enough for the
ground state sitting at the bottom of the cosine potential well to finally
enter the bunching layer together with its partner state, see Fig. 3(a).

This 1D model, however, is not sufficient to explain the situation at larger
positive detuning, compare Fig. 3(a) to Fig. 5 obtained for
$\Delta/2\pi=$1.998\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$
($\omega_{r}/2\pi=$4.3\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$). The
discrepancy is easily understood. Indeed, because the index $i_{t}^{\star}$ of
the transmon level sitting at the top of the inverted parabola increases with
the detuning, at large detuning states $i_{t}<i_{t}^{\star}$ no longer have
quasidegenerate partners on the lattice, see Fig. 4(b). The absence of these
multiphoton resonances tends to push ionization to larger photon numbers.

Moreover, at small resonator frequency with respect to the qubit frequency,
more than $N$ photons may be needed to connect pairs of states
$\ket{i_{t},n_{r}}$ associated to transmon states deep in the well to states
near the top of the well. Considering only the projection on the
$N$-excitation subspace is therefore no longer sufficient, signaling the
breakdown of the RWA. More generally, Eq. 9 does not capture the full matrix
elements of the transmon charge operator when states near or above the top of
the cosine potential well are involved (see below).

As a result, while a bunching layer can still exist, at large positive
detuning branch swapping associated with multiphoton transitions of the type
discussed in Sec. II.3 is also observed. This is shown in Fig. 5 obtained for
$\Delta/2\pi=$1.998\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$, where this
occurs for both the ground and excited branches. At large photon numbers, the
increase with $\sqrt{n_{r}}$ of the coupling matrix elements is eventually
sufficient to overcome the large energy separations. This leads to significant
hybridization, causing the branches to merge with the bunching layer.

Figure 5: Positive detuning. Average transmon population versus $N_{r}$ for a
larger positive detuning
$\Delta/2\pi=$1.998\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ corresponding to
$i_{t}^{\star}=8$. Branches $0_{t}$ and $1_{t}$ are highlighted in color, the
others are gray.

To summarize this section, the key ingredients leading to ionization at
positive detuning are:

  * •

Because of the transmon’s negative anharmonicity, there generically exists a
transmon state $i_{t}^{\star}$ for which the dispersive approximation breaks
down at very low photon number. This results in rapid hybridization of the
transmon states and to branch bunching close to $N_{t}\sim i_{t}^{\star}$.

  * •

Branch bunching, and thus ionization, is precipitated by multiphoton
resonances becoming relevant with the increase in photon number. This is
captured within the RWA Khezri  _et al._ (2023).

  * •

At larger positive detuning, the RWA fails and weaker processes become
dominant. As a result, multiphoton resonances involving the computational
states are less likely. Generally, the photon number threshold for ionization
increases with the detuning.

###  II.5 Gate-charge dependence

An important feature of the transmon is the exponentially small charge
dispersion of its computational states with increasing $E_{J}/E_{C}$ Koch  _et
al._ (2007). This charge dispersion, however, rapidly increases for states
near and above the top of the cosine potential well, see the vertical full
lines in Fig. 4. Because of the strong delocalization of the branches over the
bare transmon states, see Fig. 3(c-f), even branches associated with states
deep in the well exhibit charge sensitivity. As a result, the exact position
of the resonances responsible for ionization depends on the gate charge.
Moreover, these gate-charge sensitive states can be efficiently coupled to
computational states given that matrix elements such as $\langle
i_{t}|\hat{n}_{t}|i_{t}+3_{t}\rangle$ can be large. When moving away from the
symmetry points $n_{g}=0$ and 0.5, other matrix elements such as $\langle
i_{t}|\hat{n}_{t}|i_{t}+2_{t}\rangle$ also become large. Both the large charge
dispersion and the influence of the gate charge on the matrix elements
therefore significantly impacts the photon number at which ionization of the
transmon computational states occurs, even deep in the transmon regime.

This is illustrated in Fig. 6 which shows the transmon branch population
versus resonator photon number at negative detuning and for two values of the
gate charge. For $n_{g}=0$ (a), branches $7_{t}$ and $1_{t}$ form an avoided
crossing at around $N_{r}\approx 84$ resulting in the observed branch
swapping. In contrast, for $n_{g}=0.3$ (b), because the charge operator can
now connect states of the same parity, branch $11_{t}$ hybridizes early with
$9_{t}$ and, in turn, branch $9_{t}$ hybridizes with $7_{t}$. The result is a
delayed ionization of $1_{t}$ at around $N_{r}\approx 108$. Dependence on gate
charge is observed at positive detuning for the same reasons. The importance
of gate charge on the driven transmon was also pointed out in Refs. Cohen  _et
al._ (2023); Khezri  _et al._ (2023).

Figure 6: Negative detuning. Average transmon population versus $N_{r}$ for
two different gates charge values (a) $n_{g}=0$ and (b) $n_{g}=0.3$. The
parameters are the same as in Fig. 2. Branches $B_{1_{t}}$, $B_{7_{t}}$,
$B_{9_{t}}$ and $B_{11_{t}}$ are highlighted (see legend), the others are
gray. Ionization of the excited state is pushed back from $N_{r}\approx 84$ at
$n_{g}=0$ to $N_{r}\approx 108$ at $n_{g}=0.3$. Figure 7: Gate charge averaged
critical photon number versus detuning extracted from the fully quantum model
(pink), and Floquet analysis (red). The top panel corresponds to the ground
state $\bar{n}_{{\rm crit},0_{t}}^{q}$ and $\bar{n}_{{\rm crit},0_{t}}^{f}$
for (a) negative detuning and (b) positive detuning. The bottom panel is for
the qubit excited state $\bar{n}_{{\rm crit},1_{t}}^{q}$ and $\bar{n}_{{\rm
crit},1_{t}}^{f}$ for (c) negative detuning and (d) positive detuning. Shaded
areas correspond to the distribution of $n_{{\rm crit},i_{t}}^{q}$ and
$n_{{\rm crit},i_{t}}^{f}$ between the 10th and 90th percentiles over 100
realizations of $n_{g}$ for both models. The black dashed lines are the
dispersive critical photon number $n_{{\rm crit},i_{t}}$ defined in Appendix A
for each state. The colored background shows the transmon population extracted
from the gate charge averaged dynamics of a driven transmon, with steady-state
photon number $\bar{n}_{r}(t\xrightarrow{}\infty)=400$,
$\kappa/2\pi=$10\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ (positive detuning)
and $\kappa/2\pi=$1\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ (negative
detuning), see Sec. III for details. Note the change in color scale between
the top and bottom panels.

###  II.6 Critical photon number

Using the fully-quantum branch analysis, we now introduce critical photon
numbers $n_{{\rm crit},i_{t}}^{q}$ at which ionization of the transmon state
$i_{t}$ is expected to occur. These critical numbers are defined as the
minimum photon number $N_{r}$ at which the average transmon population reaches
$N_{t}=2$ for the ground state branch and $N_{t}=3$ for the excited state
branch. These choices capture large changes caused by resonances while
avoiding to register the slow and smooth increase in $N_{t}$ with photon
number that is observed before ionization, see e.g. the red lines in Fig. 6.
Given the sensitivity of ionization to gate charge, the critical photon
numbers we report are averaged over 100 realizations of $n_{g}$ between 0 and
0.5.

Figure 7 shows the resulting $n_{g}$-averaged critical photon number
$\bar{n}^{q}_{\mathrm{crit},i_{t}}$ as a function of the qubit-resonator
detuning (pink line). The top panels correspond to the ground state and the
bottom panels to the excited state. We exclude a range of detunings around
$\Delta=0$ where the dispersive approximation does not hold (dashed region).
The shaded pink area shows the distribution of $n_{{\rm crit},0_{t}}^{q}$ and
$n_{{\rm crit},1_{t}}^{q}$ between the 10th and 90th percentiles for the
realization of $n_{g}$, emphasizing the strong dependence of ionization on
gate charge. The usual Jaynes-Cummings critical photon number $n_{{\rm
crit},i_{t}}$ related to the breakdown of the dispersive approximation (see
Appendix A) is shown as the black dashed line. Apart from a few values of
$n_{g}$, the onset of ionization is well above $n_{\rm crit}$. This
observation is in agreement with Ref. Sunada  _et al._ (2022) where a 99.1%
dispersive readout fidelity was obtained at an average photon number
$\bar{n}_{r}\sim 1.5n_{\rm crit}$ larger than $n_{\mathrm{crit}}$.

As a comparison, the background shows the average transmon population obtained
from numerical integration of Schrödinger’s equation for a driven transmon. We
show in Sec. III how with the appropriate choice of time-dependent drive
amplitude this simplified approach captures the essential physics of the full
quantum model. Because of the strong dependence on $n_{g}$, the results are
also averaged over 100 realization of the gate charge. The agreement between
the critical photon numbers extracted from the branch analysis is excellent,
confirming that the $n_{{\rm crit},i_{t}}^{q}$ are useful and numerically
simple to compute proxies for the onset of ionization.

The critical numbers $n_{{\rm crit},i_{t}}^{q}$ exhibit a complex behavior as
a function of detuning which can be understood from the phenomenology
discussed in the previous sections. Focusing first on the case of negative
detuning, we first observe the expected overall increase of $n_{{\rm
crit},i_{t}}^{q}$ with $|\Delta|$ due to the decrease of the $\chi_{i_{t}}$
shifts. As these shifts become smaller in magnitude, the slope of the modular
energies decreases, pushing the energy collisions responsible for branch
swapping to larger photon number and consequently increasing $n_{{\rm
crit},0_{t}}^{q}$ and $n_{{\rm crit},1_{t}}^{q}$. On the other hand, the large
detuning-dependent dips are due to branch swapping caused by anticrossings
between the ground state [panel (a)] or excited state [panel (c)] branches
with other branches. Increasing the detuning has the effect of sweeping
through these resonances, resulting in the observed dips if the photon number
at the anticrossing is large enough. For example, the three large dips in
panel (a) result from resonances in the modular spectrum between the ground
state branch and (from small to large negative detuning) branches $B_{7_{t}}$,
$B_{6_{t}}$ and $B_{5_{t}}$, see Appendix B for details.

At positive detuning, we also observe the expected overall increase of the
critical photon numbers with detuning. On top of this behavior are sharp dips
resulting from multiphoton resonances between the ground state [panel (b)] or
excited state [panel (d)] and a quasiresonant partner state in the inverted
potential of Fig. 4(a). At larger detuning, the ground and excited states no
longer have a quasiresonant partner, see Fig. 4(b). As a result, this
mechanism precipitating ionization is no longer active, and the increase of
the critical photon number with detuning becomes more monotonic.

##  III Floquet analysis of the driven transmon

During qubit readout, the drive $\varepsilon_{d}$ on the resonator displaces
the resonator field to a coherent state with amplitude $\alpha(t)$ which, in
turn, acts as a classical drive on the transmon Lledó  _et al._ (2023); Muñoz
Arias  _et al._ (2023). By making a displacement transformation
$\hat{a}\rightarrow\hat{a}+\alpha(t)$ on the master equation, Eq. 3, and
subsequently ignoring the quantum fluctuations of the resonator, we arrive at
the simplified semiclassical picture of a driven transmon with Hamiltonian
(see Appendix C)

$\begin{split}\hat{H}(t)=&4E_{C}(\hat{n}_{t}-n_{g})^{2}-E_{J}\cos(\hat{\varphi}_{t})+\mathcal{E}_{t}(t)\hat{n}_{t}\\\
=&\hat{H}_{t}+\mathcal{E}_{t}(t)\hat{n}_{t}.\end{split}$ (10)

Assuming the resonator drive frequency $\omega_{d}$ to be close to the
resonator frequency $\omega_{r}$, the classical drive on the transmon is
approximately $\mathcal{E}_{t}(t)\approx\varepsilon_{t}(t)\cos(\omega_{d}t)$,
with a time-dependent amplitude $\varepsilon_{t}(t)=2g\sqrt{\bar{n}_{r}(t)}$
where

$\bar{n}_{r}(t)=\left(\frac{\varepsilon_{d}}{\kappa}\right)^{2}(1-e^{-\kappa
t/2})^{2}$ (11)

is the average number of photons in the resonator, see Appendix C for details.
In this model the resonator frequency $\omega_{r}$ of the full quantum model
is replaced by the drive frequency $\omega_{d}$. As a result, we place the
qubit below or above $\omega_{d}$ when considering negative or positive
detunings.

The effective drive on the transmon originates from the displacement of the
transmon-resonator coupling and its amplitude is thus proportional to the
coupling strength $g$. With the usual large values of $g$ in circuit QED, this
can lead to a very large effective drive amplitude on the transmon. For
example, for the value $g/2\pi=$120\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$
used in this paper, we have
$\varepsilon_{t}/2\pi\sim$750\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ for
$\bar{n}_{r}=10$ photons in the resonator. This is significantly larger than
the transmon anharmonicity $\alpha$. Even though it is off resonant, this
effective drive cannot be treated perturbatively and can be expected to lead
to leakage out of the computational subspace. In this section, we show that
this simplified model leads to the same predictions for transmon ionization as
the full quantum model. The essence of the success of this simple approach is
that it preserves the full cosine potential of the transmon qubit.

Because the effective drive amplitude
$\varepsilon_{t}(t)\propto\sqrt{\bar{n}_{r}(t)}$ changes on a time scale
$1/\kappa$ that is much larger than the period of the drive
($T=2\pi/\omega_{d}$), our analysis relies on the instantaneous Floquet
spectrum Breuer and Holthaus (1989). Considering $\hat{H}(t)$ to be periodic,
in this approach the solution to Schrödinger’s equation can be expressed as
$|\psi(t)\rangle=\sum_{i_{t}}c_{i_{t}}e^{-i\epsilon_{i_{t}}t}|\phi_{i_{t}}(t)\rangle$
where $\epsilon_{i_{t}}$ is a Floquet quasienergy and
$|\phi_{i_{t}}(t)\rangle=|\phi_{i_{t}}(t+T)\rangle$ a Floquet mode Grifoni and
Hänggi (1998). The phases $\exp(-i\epsilon_{i_{t}}T)$ and the Floquet modes
$\ket{\phi_{i_{t}}(t)}$ are the eigenvalues and eigenvectors of the propagator
over one period of the drive,
$\hat{U}(t+T,t)|\phi_{i_{t}}(t)\rangle=e^{-i\epsilon_{i_{t}}T}|\phi_{i_{t}}(t)\rangle$.
These eigenvalues can also be obtained for a static Hamiltonian whose form is
reminiscent of the transmon-resonator Hamiltonian, see Appendix D.
Importantly, as the phases $\exp(-i\epsilon_{i_{t}}T)$ remain invariant under
a shift of $\epsilon_{i_{t}}$ by an integer multiple of $\omega_{d}$, the
quasienergies $\epsilon_{i_{t}}$ are only defined modulo $\omega_{d}$.
Quasienergies are thus analogous to quasimomenta in Bloch theory and they are
folded in a first ‘Brillouin zone’
$-\omega_{d}/2\leq\epsilon_{i_{t}}\leq\omega_{d}/2$. Importantly, the Floquet
modes and quasienergies are functions of the drive amplitude
$\varepsilon_{t}$.

Just as in the previous section, we sort the Floquet quasienergies
$\epsilon_{i_{t}}[\varepsilon_{t}]$ and modes
$|\phi_{i_{t}}[\varepsilon_{t}]\rangle$ into transmon branches. Our labeling
procedure at finite drive is performed by smoothly increasing the drive
amplitude to connect the eigenstates of the unitary operator to those at zero
drive. This is made possible by the fact that at zero drive, Floquet modes are
simply the bare transmon eigenstates
$\ket{\phi_{i_{t}}[\varepsilon_{t}=0]}=\ket{i_{t}}$ and the quasienergies the
bare transmon energies modulo the drive frequency
$\epsilon_{i_{t}}[\varepsilon_{t}=0]=E_{i_{t}}\mod\omega_{d}$. Given the above
observations, the Floquet quasienergy spectrum is thus akin to the modular
energy spectrum introduced above for the full quantum model. Both spectra
capture changes in (quasi)energy with photon number and, thanks to the folding
of the spectrum present in both cases, they can display avoided crossing
between states corresponding to low- and high-energy transmon states, see e.g.
Fig. 8(a). Avoided crossings in the quasienergy spectrum are thus linked to
resonances between ac-Stark shifted transmon states, up to $k$ drive photons
of energy $\omega_{d}$.

For unbounded Hamiltonians such as the transmon Hamiltonian, however, there
are avoided crossings with quasienergy gaps of arbitrarily small sizes Hone
_et al._ (1997). To avoid very small avoided crossings that are not relevant
to transmon ionization, we numerically obtain the Floquet spectrum as a
function of $\varepsilon_{t}$ by setting a finite increment
$\delta\varepsilon_{t}$ for tracking Floquet branches $B_{i_{t}}$. This
contrasts with the branch construction of the full quantum model where, by
default, one resonator photon is added at each step, corresponding to an
increment $\sim g/\sqrt{n_{r}}$. Here, we take
$\delta\varepsilon_{t}/2\pi=$10\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ which
allows to capture the avoided crossing seen in the modular spectrum of the
full quantum model. Another reason for this choice is that the effective drive
amplitude $\varepsilon(t)$ is linked to the photon population of the resonator
which in a qubit readout goes from 0 to a finite population $n_{r}(t)$ in a
time scale set by $\kappa$. As such, the size of the increment
$\delta\varepsilon_{t}$ can, to some extent, be linked to the system’s
dynamics. We have thus taken $\delta\varepsilon_{t}$ to be close to the value
of $\kappa/2\pi=$7.95\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ used in our
time-dependent simulations.

###  III.1 Negative detuning: $\omega_{q}<\omega_{r}$

Figure 8: Negative detuning. (a) Quasienergy spectrum and (b) period-averaged
mean excitation number of the Floquet modes as a function of the mean
resonator photon number $\bar{n}_{r}=(\varepsilon_{t}/2g)^{2}$ for a negative
detuning $\Delta/2\pi=-$1.217\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$.
Branches ${0_{t}}$, ${1_{t}}$, ${7_{t}}$, ${10_{t}}$ and ${11_{t}}$ are
highlighted in color (see legend), the others are gray. The transmon
parameters are the same as in Fig. 2. (c) Time evolution of the mean number of
excitation under a drive of amplitude
$\varepsilon_{t}(t)=2g\sqrt{100}(1-e^{-\kappa t/2})$ with
$\kappa/2\pi=$7.95\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ and frequency
$\omega_{d}/2\pi=$7.515\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$. The transmon
is initialized in $\ket{0_{t}}$ (blue line) or $\ket{1_{t}}$ (red line). (d)
Probability
$|\bra{\phi_{i_{t}}[\varepsilon_{t}(t_{f})]}\hat{\mathcal{U}}(t_{f})\ket{1_{t}}|^{2}$
of being in the Floquet mode $i_{t}=1_{t}$ (red line, adiabatic) or
$i_{t}=7_{t}$ (light blue line, diabatic) vs. $\kappa$ at the final time
$t_{f}\approx 10/\kappa$ of the evolution. The dashed black lines are the
Landau-Zener predictions, see Appendix E.

We saw in Sec. II.3 that the opposite signs of the ac-Stark shifts of pairs of
transmon branches (one deep in the cosine potential well and the other near
the top of the well) results in avoided crossings in the modular energies.
These avoided crossings — which can be precipitated by the early hybridization
of states at the top of the well — result in transmon ionization during
readout. We now show that the same phenomenology is present in the simplified
model of the classically driven transmon, Eq. 10, with its full cosine
potential.

Figure 8(a) shows the numerically-obtained quasienergy spectrum as a function
of the transmon drive amplitude here expressed in terms of the equivalent mean
resonator photon number $\bar{n}_{r}=(\varepsilon_{t}/2g)^{2}$. The similarity
with the modular energy spectrum of Fig. 2(b) obtained from the full quantum
model is remarkable. In correspondence with Fig. 2(b), the quasienergies of
states at the bottom and at the top of the cosine potential well move in
opposite directions with $\bar{n}_{r}$ (see Appendix D), leading to avoided
crossings at finite drive amplitudes. As an example, the quasienergies of
branches $1_{t}$ (red line) and $7_{t}$ (light blue line) approach each other
until an avoided crossing is formed at $\bar{n}_{r}\approx 82$ photons.

To illustrate the abrupt change in the character of the Floquet modes at the
avoided crossings, we show in Fig. 8(b) the period-averaged transmon
excitation,

$\langle\langle\hat{N}_{t}\rangle\rangle_{i_{t}}=\frac{1}{T}\int_{0}^{T}d\tau\,\sum_{j_{t}}j_{t}|\langle
j_{t}|\phi_{i_{t}}(\tau)\rangle|^{2},$ (12)

for each Floquet branch as a function of $\bar{n}_{r}$. The similarity with
the transmon population of the branches computed for the full transmon-
resonator system of Fig. 2(a) is again remarkable, with branch swapping
occurring at the resonances identified in Fig. 8(a). Note also the presence of
small features at $\bar{n}_{r}\approx 22$ of the period-averaged transmon
excitation corresponding to the branches $0_{t}$ (dark blue line) and $7_{t}$
(light blue line). This feature is due to a very small avoided crossing of the
same two branches also at $\bar{n}_{r}\approx 22$, which, with our choice of
increment $\delta\varepsilon_{t}$, is not resolved and is shown as a crossing,
see Fig. 8(a). We confirm below that with our choice of $\kappa$ this
resonance does not play a role in the dynamics of the system.

To confirm that the observed resonances are responsible for transmon
ionization in this semiclassical model, we now turn to time dynamics obtained
under the Hamiltonian of Eq. 10 together with resonator decay $\kappa$, see
Fig. 8(c). Here, $\varepsilon_{t}$ is no longer a static control parameter but
changes in time following $\varepsilon_{t}(t)=2g\sqrt{\bar{n}_{r}(t)}$ with
the average photon number given by Eq. 11. After preparing the transmon in the
bare state $\ket{\psi(0)}=\ket{i_{t}}$ with $\varepsilon_{t}(0)=0$, the state
adiabatically follows the instantaneous Floquet mode of the $i_{t}$-th branch,
$\ket{\phi_{i_{t}}[\varepsilon_{t}(t)](t)}$, as the drive amplitude is
increased and until an avoided crossing with branch $j_{t}$ is met. Following
the usual Landau-Zener argument Grifoni and Hänggi (1998), at that point the
state diabatically switches branch or adiabatically remains in the branch.
This transition occurs with the Landau-Zener probability,
$P_{\mathrm{LZ}}=\exp(-\pi\Delta_{\rm ac}^{2}/2v)$, where $\Delta_{\rm ac}$ is
the quasienergy gap at the avoided crossing and $v(\varepsilon_{t})$ the speed
of passage through the resonance Drese and Holthaus (1999); Oliver  _et al._
(2005); Shevchenko  _et al._ (2010); Ikeda  _et al._ (2022), see Appendix E.
Importantly, ionization is averted if the state transitions diabatically
through the crossing. By contrast, transitioning adiabatically means that the
transmon ionizes.

In Fig. 8(c) we show a parametric plot of the transmon population as a
function of time, here expressed in terms of $\bar{n}_{r}(t)$. When initially
prepared in $\ket{0_{t}}$ (blue line), the state of the transmon follows
branch $0_{t}$, going straight through the weak avoided crossing at
$\bar{n}_{r}\sim 22$ that our finite-step-size Floquet tracking purposely did
not capture. In this case, the drive does not induce transitions in the
transmon. On the other hand, when initially prepared in $\ket{1_{t}}$ (red
line), the state follows mostly adiabatically branch $1_{t}$ through the
avoided crossing at $\bar{n}_{r}\approx 82$ leading to a transition to high-
energy transmon states. Figure 8(d) shows the probability of transitioning
adiabatically or diabatically through the avoided crossing as a function of
the rate at which $\varepsilon_{t}(t)$ varies. More precisely, we show as a
function of $\kappa$ [which controls $\varepsilon_{t}(t)$, see Eq. 11] the
probability of the state being in branch $i_{t}=1_{t}$ (red line) or in branch
$i_{t}=7_{t}$ (light blue line) at the final time $t_{f}\approx 10/\kappa$ of
the evolution. By increasing the value of $\kappa$ and thus the speed at which
resonances are crossed, the Landau-Zener probability for a diabatic crossing
can be tuned from zero to almost one Shillito  _et al._ (2022). On top of
these curves, we show the transition probability obtained from the Floquet-
Landau-Zener formula (dashed black line) in perfect agreement with the time-
dynamics simulations, see Appendix E.

The onset of ionization therefore depends on the effective speed at which
resonances are crossed. According to Eq. 11, the drive ramp-up is faster at
low photon numbers and slows down when approaching the steady-state photon
number. We thus expect resonances that occur at the steady-state photon number
to be more easily resolved. Consequently, ionization can be pushed back by
driving the resonator stronger and thus diabatically crossing the resonance.
This is observed experimentally in Ref. Sank  _et al._ (2016), wherein the
probability of qubit ionization exhibits an initial rise at a certain drive
power, followed by a subsequent decline at higher powers.

###  III.2 Positive detuning : $\omega_{q}>\omega_{r}$

Because of the negative anharmonicity of the transmon, the branch swapping
observed at negative detuning is replaced by branch bunching at positive
detuning. We now show that this phenomenology is captured by the simplified
model of the classically driven transmon. Because ionization typically occurs
at smaller resonator photon numbers at positive detuning than at negative
detuning, larger discrepancies between the fully quantum model and the
semiclassical model are expected.

Figure 9:  Positive detuning. (a) Period-averaged mean excitation number of
the Floquet modes, (b) Quasienergy spectrum, and (c)-(f) projection of the
Floquet modes $7_{t}$, $1_{t}$, $8_{t}$ and $0_{t}$ on the transmon basis in
logarithmic scale as a function of the mean resonator photon number
$\bar{n}_{r}=(\varepsilon_{t}/2g)^{2}$ for a positive
$\Delta/2\pi=$1.031\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$. In (a) and (b),
branches ${0_{t}}$, ${1_{t}}$, ${4_{t}}$, ${7_{t}}$, and ${8_{t}}$ are
highlighted in color (see legend), the others are gray. The transmon
parameters are the same as in Fig. 3. (g) Time evolution of the mean number of
transmon excitations under a drive of amplitude
$\varepsilon_{t}(t)=2g\sqrt{30}(1-e^{-\kappa t/2})$ with
$\kappa/2\pi=$11.98\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ (set to match the
full $\chi$ shift at zero photons), and frequency
$\omega_{d}/2\pi=$5.267\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$. The transmon
is initialized in $\ket{0_{t}}$ (blue line) or $\ket{1_{t}}$ (red line).

Figure 9(a) and (b) shows the period-averaged transmon excitation of the
Floquet branches and the quasienergy spectrum, respectively, versus
$\bar{n}_{r}$ and obtained for the same parameters as in Sec. II.4. Although
there are some discrepancies especially at small photon numbers, the agreement
with the modular spectrum of the full quantum model is nevertheless excellent,
see Fig. 3. As in the full quantum model, rather than branch swapping we
observe in Fig. 9(a) the expected branch bunching at
$\langle\langle\hat{N}_{t}\rangle\rangle\sim 4$. The branch bunching also
correctly correlates with anticrossings in the quasienergy spectrum. For
example, at $\bar{n}_{r}\approx 6.4$, the quasienergies of branches $1_{t}$
and $7_{t}$ form a strong and wide avoided crossing, at which point we also
observe the bunching of these two branches.

Interestingly, a weak avoided crossing of the quasienergies associated to
$0_{t}$ and $8_{t}$ at $\bar{n}_{r}\approx 5.3$ is not captured with our
finite tracking increment of
$\delta\varepsilon_{r}/2\pi=$10\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$.
While the $1_{t}-7_{t}$ transition requires the absorption of six drive
photons, the $0_{t}-8_{t}$ transition requires eight. It is hence more
strongly suppressed and is not observed at this low drive amplitude. This
higher-order transition, however, becomes relevant at larger drive amplitudes.
Indeed, at $\bar{n}_{r}\approx 18$ a weak level repulsion of the quasienergies
associated to $0_{t}$ and $8_{t}$ are observed and which is correlated with
the start of the bunching of these two states.

Those correlations are confirmed by the projection of the Floquet modes onto
the bare transmon basis, see Fig. 9(c-f). There, we see that the Floquet modes
associated to $0_{t}$ and $8_{t}$ show a strong delocalization in the transmon
bare basis at $\bar{n}_{r}\approx 18$. On the other hand, for $1_{t}$ and
$7_{t}$ this occurs at $\bar{n}_{r}\approx 6.5$, matching the presence of the
avoided crossing in the quasienergy spectrum. To quantify the delocalization
of the modes swallowed into the bunching layer, we compute the inverse
participation ratio $\sum_{j_{t}}|\langle j_{t}|\phi_{i_{t}}\rangle|^{4}$
Kramer and MacKinnon (1993) which here measures the delocalization of a
Floquet mode $\ket{\phi_{i_{t}}}$ on the transmon basis. At
$\bar{n}_{r}\approx 11$, the inverse participation ratio is approximately
$0.14\approx 1/7$ for both $\ket{\phi_{7_{t}}}$ and $\ket{\phi_{1_{t}}}$,
corresponding to a delocalization of these modes between at least seven
transmon states.

To confirm that the presence of avoided crossings in the quasienergy spectrum
is an accurate predictor of the system’s dynamics, in Fig. 9(g) we show the
time-dependent transmon population as a function of time, here parameterized
by the instantaneous average photon number $\bar{n}_{r}(t)$. When the transmon
is initialized in $\ket{0_{t}}$ (blue line), the state follows branch $0_{t}$
avoiding an early measurement-induced ionization at $\bar{n}_{r}\approx 5.3$.
As a result, no feature is observed in $\langle\hat{N}_{r}\rangle(t)$ at that
average photon number, except for a slow overall increase of the transmon
population due to the gradual hybridization of branch $0_{t}$. At
$\bar{n}_{r}\approx 18$, the transmon population increases, as expected from
the delocalization of that branch, see Fig. 9(f). On the other hand, when
initialized in $\ket{1_{t}}$ (red line), the system mostly adiabatically
transitions at the large avoided crossing at $\bar{n}_{r}\approx 6.4$, thereby
following branch $1_{t}$ with a high probability. This results in measurement-
induced transitions and to the observed rapid increase of
$\langle\hat{N}_{t}\rangle(t)$.

###  III.3 Comparison to fully quantum model

We define the critical photon numbers in the Floquet analysis in a similar way
as in the branch analysis of the fully quantum model. The Floquet critical
photon numbers for the ground state $n_{{\rm crit},0_{t}}^{f}$ and for the
excited state $n_{{\rm crit},1_{t}}^{f}$ are thus respectively defined as the
smallest value of $\bar{n}_{r}$ for which
$\langle\langle\hat{N}_{t}\rangle\rangle_{0_{t}}=2$ and
$\langle\langle\hat{N}_{t}\rangle\rangle_{1_{t}}=3$. In Fig. 7, we compare the
extracted $n_{{\rm crit},i_{t}}^{f}$ (red line) with $n_{{\rm
crit},i_{t}}^{q}$ of the fully quantum model (pink line) for a range of
detunings. The agreement between the two methods is excellent for both the
average values and their fluctuations with charge noise (see the red and pink
shaded areas). This agreement confirms the validity of treating the driven
resonator as an effective classical field driving the transmon. Numerically,
this approach is advantageous because of it simplicity and efficiency.

As previously mentioned, in Fig. 7 we also compare the critical photon numbers
with the dynamics of the driven transmon (colored background). Because avoided
crossings are very narrow in the negative detuning regime, we take a small
$\kappa/2\pi=$1\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ to ensure a slow
ramp-up of the average photon number in the resonator. In this way, the system
crosses slowy the resonances responsible for $\bar{n}_{{\rm crit},i_{t}}^{q}$
and $\bar{n}_{{\rm crit},i_{t}}^{f}$, leading to ionization in the time
dynamics. At positive detunings, resonances are wider and a faster ramp-up
with $\kappa/2\pi=$10\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$ still results
in ionization around $\bar{n}_{{\rm crit},i_{t}}^{q}$ and $\bar{n}_{{\rm
crit},i_{t}}^{f}$. When comparing to experiments, the value of $\kappa$ and
$\varepsilon_{d}$ used in the time-dynamical simulations can be adjusted to
match the speed at which the resonator field builds up. For large speeds,
early weak resonances will be crossed rapidly and therefore not lead to
ionization. The step size $\delta\varepsilon_{t}$ in the Floquet tracking
should be adjusted accordingly.

##  IV Nonlinear dynamics of a driven classical transmon

We have shown in Sec. III that a simplified model of a transmon driven by a
classical monochromatic drive is enough to understand and predict the onset of
ionization. In this section, we go one step further and investigate a model
where the transmon itself is treated as a classical object. As in the previous
sections, we account for the full nonlinearity of the transmon cosine
potential and treat the readout drive nonperturbatively. This fully classical
description, combined with the Bohr-Sommerfeld quantization rule, was shown in
Ref. Cohen  _et al._ (2023) to capture the main features of ionizing
transitions. Here, we further clarify the classical mechanisms of ionization
and explicitly connect them to the quantum phenomenology discussed in Sec. II
and Sec. III.

###  IV.1 General properties of the classical model

Figure 10:  Poincaré sections for different values of the drive amplitude
$\tilde{\varepsilon}_{t}$ at negative (top row) and positive (bottom row)
detunings. For negative detuning (a-d) the drive frequency is
$\tilde{\omega}_{d}=1.152$
($\omega_{d}/2\pi=$7.515\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$,
$\Delta/2\pi=-$1.202\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$) and the drive
amplitudes are, from left to right,
$\tilde{\varepsilon}_{t}\in\\{0,0.167,0.501,0.624\\}$
($\bar{n}_{r}\in\\{0,20.6,185.6,287.6\\}$). For positive detuning (e-h) the
drive frequency is $\tilde{\omega}_{d}=0.807$
($\omega_{d}/2\pi=$5.267\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$,
$\Delta/2\pi=$0.998\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$) and the drive
amplitudes are, from left to right,
$\tilde{\varepsilon}_{t}\in\\{0,0.029,0.081,0.152\\}$
($\bar{n}_{r}\in\\{0,0.63,4.88,17.02\\}$). The Bohr-Sommerfeld orbits of the
qubit ground and excited states are shown for $E_{J}/E_{C}=110$ in dark blue
and dark red, respectively. At zero drive amplitude, (a) and (e), the
separatrix is shown in olive. In (b) $(5\\!:\\!3)$ and $(7\\!:\\!5)$
resonances are highlighted in purple and plum red, respectively. A pair of
$(3\\!:\\!2)$ resonances are also highlighted in orange and pale green. In (c)
the $(7\\!:\\!5)$ resonance collides with the Bohr-Sommerfeld orbit of the
excited state. This is illustrated by coloring the resonance in red. In (d),
merging of the excited state with the chaotic layer is assisted by a pair of
$(3\\!:\\!2)$ resonances that reduce the regular area available to the excited
states. In (f-h) a $(1\\!:\\!1)$ resonance occupies a large fraction of the
region inside the main separatrix. At large drive amplitudes, (d) and (h), the
excited state is engulfed by the chaotic layer. This is illustrated by the red
points in the chaotic layer.

The full nonlinear dynamics of the classical transmon is governed by a
dimensionless Hamiltonian that is formally identical to that of a driven
classical pendulum Koch  _et al._ (2007)

$\begin{split}\tilde{H}(\tilde{t})&=\frac{1}{2}\tilde{n}_{t}^{2}-\cos{\tilde{\varphi}_{t}}+\tilde{\varepsilon}_{t}\cos({\tilde{\omega}_{d}\tilde{t}})\tilde{n}_{t}\\\
&=\tilde{H}_{t}+\tilde{\varepsilon}_{t}\cos({\tilde{\omega}_{d}\tilde{t}})\tilde{n}_{t}.\end{split}$
(13)

Here, $\tilde{H}_{t}=H_{t}/E_{J}$ and $\tilde{t}=\omega_{p}t$ are rescaled
energy and time, with $\omega_{p}=\sqrt{8E_{J}E_{C}}$ the transmon plasma
frequency. Moreover, $\tilde{\varepsilon}_{t}=\varepsilon_{t}/\omega_{p}$, and
$\tilde{\omega}_{d}=\omega_{d}/\omega_{p}$ are the rescaled transmon drive
amplitude and drive frequency, respectively. The rescaled phase-space
coordinates are $\tilde{\varphi}_{t}=\varphi_{t}$ and $\tilde{n}_{t}=zn_{t}$
with Poisson bracket $\\{\tilde{\varphi}_{t},\tilde{n}_{t}\\}=z$, where
$z=\sqrt{8E_{C}/E_{J}}$ is the transmon impedance Cohen  _et al._ (2023). We
have omitted the gate charge $n_{g}$ from Eq. 13 since a static gate charge
does not affect classical dynamics. Moreover, as in Sec. III, the drive
frequency $\omega_{d}$ is a surrogate for the resonator frequency of the full
quantum model.

Given the periodic time dependence of Eq. 13, the solutions to Hamilton’s
equations are best studied using the Poincaré section obtained by plotting the
phase-space coordinates stroboscopically at multiples of the drive period
$\tilde{T}=2\pi/\tilde{\omega}_{d}$. In Fig. 10, we show the Poincaré sections
for a negative detuning (top panels) and for a positive detuning (bottom
panels) of the drive to the qubit frequency. In each case, we show the
Poincaré sections for increasing values of the drive amplitude
$\tilde{\varepsilon}_{t}$. These results are obtained from the same values of
the transmon parameters as above, here corresponding to
$\omega_{p}/2\pi=$6.526\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ and $z\sim
0.27$.

When $\tilde{\varepsilon}_{t}=0$, two distinct types of regular orbits fill
phase space, see panels (a) and (e). For $\tilde{H}_{t}<1$, the system
undergoes bounded periodic phase oscillations that are the classical analog of
the transmon bound states. In the Poincaré section, these bounded oscillations
correspond to the near circular orbits located in the center of phase space.
Below, we will refer to these orbits as the main regular region. On the other
hand, for $\tilde{H}_{t}>1$ the system undergoes full $\pm 2\pi$ phase
rotations, corresponding to phase slips in the quantum model Koch  _et al._
(2007). This type of motion appears in the Poincaré section as nearly
horizontal features that do not cross $\tilde{n}_{t}=0$. Separating these two
types of motion is the contour $\tilde{H}_{t}=1$ known as the separatrix
[olive line in Fig. 10(a,e)]. In the vicinity of the separatrix, small changes
to the initial conditions can lead to large changes in the system’s dynamics.

In addition to bounded and unbounded oscillations, for
$\tilde{\varepsilon}_{t}>0$ new types of motion are visible in the Poincaré
sections. The first type of motion common to all panels of Fig. 10 where
$\tilde{\varepsilon}_{t}>0$ is chaotic motion. Chaos emerges in the vicinity
of the separatrix for an arbitrarily weak drive amplitude and results from the
instability of the system with respect to small perturbation at energies close
to the separatrix. The size of the phase space region covered by chaotic
motion generally increases with the drive amplitude Zaslavsky (2005);
Wimberger (2014). Critically, this process reduces the area of phase space
supporting regular bounded oscillatory motion. This reduction of phase space
area available to the main regular region due to the emergence of chaotic
layers was identified as a key phenomenon responsible for the onset of
ionization in the classical model Cohen  _et al._ (2023).

The second type of motion seen in all panels of Fig. 10 for which
$\tilde{\varepsilon}_{t}>0$ is the presence of sets of islands surrounding the
main regular region, see, e.g., the features highlighted in plum red, olive,
and purple in panel (b) and the large croissant-shaped feature in panels
(f,g,h). These features are subharmonic periodic orbits corresponding to
$(n\\!:\\!m)$-nonlinear resonances occurring when
$n\tilde{\omega}(\tilde{I})=m\tilde{\omega}_{d}$ Zaslavsky (2005). In this
expression,
$\tilde{\omega}(\tilde{I})=\partial\tilde{H}_{t}/\partial\tilde{I}$, with
$\tilde{I}$ the rescaled action variable, is the energy-dependent frequency of
periodic motion of the pendulum which reduces to 1 at low energy. An
$(n\\!:\\!m)$ resonance can be identified in phase space by a chain of $n$
islands. We show below that in addition to the emergence of chaos the presence
of these resonances is a second key phenomenon responsible for the onset of
ionization within the classical model. Importantly, the phase-space area
occupied by a resonance depends on $n$, with increasing values of $n$ leading
to smaller resonances, i.e., weaker resonances. The area occupied by a
resonance and its distance from the center of the main regular region also
generally depends on $\tilde{\varepsilon}_{t}$. As a result, although there is
typically an abundance of $(n\\!:\\!m)$ resonances due to the nonlinearity of
$\tilde{\omega}(\tilde{I})$ Zaslavsky (2005), not all of those play an equally
important role in ionization.

###  IV.2 Bohr-Sommerfeld critical photon numbers

Figure 11:  Cut $(\tilde{\varphi}_{t}=0,\tilde{n}_{t})$ of phase space as a
function of $\tilde{\varepsilon}_{t}$ for (a) negative detuning and
($\tilde{\omega}_{d}=1.152$,
$\omega_{d}/2\pi=$7.515\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$,
$\Delta/2\pi=-$1.202\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$) (b) positve
detuning ($\tilde{\omega}_{d}=0.807$,
$\omega_{d}/2\pi=$5.267\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$,
$\Delta/2\pi=$0.998\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$). The Bohr-
Sommerfeld orbits are represented by pale colored lines and dark regions with
matching colors quantum fluctuations of $\pi\hbar_{\rm{eff}}$ on each side of
the Bohr-Sommerfeld orbits. The gray line centered at $\tilde{n}_{t}=0$ at
zero drive indicates the center of the main regular region, while the light
blue and light red lines closest to it are the Bohr-Sommerfeld orbits of the
qubit ground and excited states. In (b), the gray line at
$\tilde{n}_{t}\approx 1.5$ indicates the center of the $(1\\!:\\!1)$
resonance. White regions indicate points where the algorithm of Appendix G
fails to find a closed connected orbit due to chaotic or resonant motion. The
points where the Bohr-Sommerfeld orbits of the qubit states stop existing
(stars) yield the Bohr-Sommerfeld critical photon numbers
$\\{n_{\mathrm{crit},i_{t}}^{c}\\}$. As discussed in the text, weak resonances
are ignored. Panels (c) and (d) show the Poincaré sections using the same
color code corresponding to the drive amplitudes
$\tilde{\varepsilon}_{t}=0.167$ ($\bar{n}_{r}=20.6$) and
$\tilde{\varepsilon}_{t}=0.029$ ($\bar{n}_{r}=0.63$) indicated by the dashed
yellow lines in panels (a) and (b), respectively. For $E_{J}/E_{C}=110$ the
effective Planck constant is $2\pi\hbar_{\rm eff}=1.69445$. The horizontal
axis scans initial conditions $(\tilde{\varphi}_{t}=0,\tilde{n}_{t}^{0})$ with
$\tilde{n}_{t}^{0}\in[-2.0,2.0]$ by increments of $1/300$. The vertical axis
scans the drive amplitudes in increments of
$\delta\tilde{\varepsilon}_{t}=0.0015323$, corresponding to an increment of
$10\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$.

We now exploit the previous observations to predict the onset of ionization
within the classical model. Interestingly, the relevance of chaotic behavior
to the ionization of highly-excited Rydberg atoms was already pointed out
Delone  _et al._ (1983); Casati (1987); Wang and Chu (1989). Here, we follow
Ref. Cohen  _et al._ (2023) and use Bohr-Sommerfeld quantization to establish
a correspondence between the strobsocopically obtained orbits in the Poincaré
section and the Floquet modes of the transmon.

Since the rescaling of the dynamical coordinates in Eq. 13 does not preserve
volumes in phase space, the quantization must proceed using an effective
reduced Planck constant $\hbar_{\mathrm{eff}}=z=\sqrt{8E_{C}/E_{J}}$ Cohen
_et al._ (2023). The uncertainty principle then dictates that $2\pi\hbar_{\rm
eff}$ is the area of phase space occupied by a single quantum state. As such,
it imposes a limit to the resolution with which quantum states can probe the
continuous structure of the classical phase space. While the classical limit
is recovered by taking $\hbar_{\rm eff}\to 0$, we will see that classical
dynamics remain a good predictor of ionization for large but finite values of
$E_{J}/E_{C}$.

The quantized Bohr-Sommerfeld orbits for oscillatory motion are the classical
orbits which enclose an area equal to Messiah (2020)

$\mathcal{A}_{i_{t}}=2\pi\hbar_{\mathrm{eff}}\left(i_{t}+\frac{1}{2}\right),$
(14)

where $i_{t}\in\\{0,1,2,3,\dots\\}$ is the Bohr-Sommerfeld quantum number for
the transmon. The number of Bohr-Sommerfeld states contained in a given area
$\mathcal{A}$ of phase space is thus given by
$\lfloor\mathcal{A}/2\pi\hbar_{\rm{eff}}\rfloor$. As an illustration, the
orbits associated with the qubit ground and excited states are the blue and
red orbits in Fig. 10. They are the orbits of the main regular region which
enclose areas $\mathcal{A}_{0}=\pi\hbar_{\mathrm{eff}}$ and
$\mathcal{A}_{1}=3\hbar_{\mathrm{eff}}$, respectively. To account for quantum
fluctuations, in our analysis we extend the area of Bohr-Sommerfeld states to
include the area between $\mathcal{A}_{i_{t}}-\pi\hbar_{\rm{eff}}$ and
$\mathcal{A}_{i_{t}}+\pi\hbar_{\rm{eff}}$. Note that because we are quantizing
orbits obtained stroboscopically, the resulting Bohr-Sommerfeld orbits are the
classical analog of the Floquet modes. We provide numerical evidence for this
identification in Appendix F.

We expect the transmon drive to induce unwanted transitions for transmon state
$i_{t}$ whenever the Bohr-Sommerfeld orbit for that state cannot exist. Both
chaotic dynamics and the presence of a $(n\\!:\\!m)$ resonance can prevent the
existence of a given Bohr-Sommerfeld orbit. Indeed, the widening of the
chaotic layer with increasing $\tilde{\varepsilon}_{t}$ reduces the number of
Bohr-Sommerfeld orbits supported within the main regular region. In addition,
a Bohr-Sommerfeld orbit can collide with a resonance for some values of
$\tilde{\varepsilon}_{t}$. Moreover, as discussed further below, even a
resonance appearing in the chaotic layer away from the qubit ground or excited
orbits can facilitate ionization by reducing the regular area available to the
qubit states.

Because it is difficult to visualize from Poincaré sections how features in
phase space change with the drive amplitude, we instead plot in Fig. 11(a-b) a
single cut $(\tilde{\varphi}_{t}=0,\tilde{n}_{t})$ of phase space as a
function of $\tilde{\varepsilon}_{t}$. All three types of motions are easily
identified on this figure. Focusing for the moment on panel (a), regular
orbits correspond to the solid colored regions. The center of each Bohr-
Sommerfeld orbits are represented as pale colored lines and the horizontal
extent of the colored region about those lines accounts for quantum
fluctuations. Notably, the blue and red regions correspond to the qubit ground
and excited states, respectively. The chaotic layer appears as ‘noise’, and
the nonlinear resonances as white gaps in the regular regions with wider gaps
corresponding to stronger resonances.

At large enough drive amplitudes, all Bohr-Sommerfeld orbits are eventually
absorbed by the chaotic layer and the colored regions disappear in the noise.
This occurs earlier for orbits of large radius and last for the ground state
orbit. Moreover, merging of the regular region with the chaotic layer is
sometimes precipitated by resonances, e.g., see the large almost horizontal
features close to $\tilde{\varepsilon}_{t}=0.6$ and 0.9 in Fig. 11(a).
Consequently, the regions where regular orbits exist often have a jagged shape
reminiscent of a Christmas tree.

Another important observation is that the distance of resonances from the
center of the main regular region goes down with the drive amplitude. As a
result, resonances within the regular region collide with the Bohr-Sommerfeld
orbits at some values of $\tilde{\varepsilon}_{t}$. For a given transmon state
$i_{t}$, the values of
$\bar{n}_{r}=(\omega_{p}/2g)^{2}\tilde{\varepsilon}_{t}^{2}$ at which the
center of the Bohr-Sommerfeld orbit $i_{t}$ collides with a resonance or with
the chaotic layer yield a set of critical photon numbers
$\\{n_{\mathrm{crit},i_{t}}^{c}\\}$. These collisions are indicated by stars
for $0_{t}$ and $1_{t}$. Below, we assign an uncertainty on the value to those
critical photons number by accounting for quantum fluctuations, i.e., the
horizontal extent of the colored regions. Appendix G describes the systematic
procedure that was used to build the Christmas trees, as well as to extract
from them the values of $\\{n_{\mathrm{crit},i_{t}}^{c}\\}$ and of the lower
bound that can be expected due to quantum fluctuations.

We now discuss how the above framework enables us to understand and predict
ionization thresholds for the case of negative and positive detunings. In
doing so, we connect this classical understanding to the quantum phenomenology
described in the previous sections.

###  IV.3 Negative detuning: $\omega_{q}<\omega_{r}$

We first consider the case of negative qubit-drive detuning, see Fig. 10(a-d).
At any finite drive amplitude $\tilde{\varepsilon}_{t}$, a chaotic layer
emerges in the vicinity of the separatrix of $\tilde{H}_{t}$, see panel (b)
which is obtained at the small drive amplitude $\tilde{\varepsilon}_{t}=0.167$
corresponding to $\bar{n}_{r}=20.6$ photons. As a result, the area of the main
regular region shrinks with increasing drive amplitude at the expense of the
chaotic layer. Moreover, several nonlinear resonances can be identified in the
chaotic or regular regions, and we have highlighted in color orbits associated
with $(5\\!:\\!3)$, $(3\\!:\\!2)$, and $(7\\!:\\!5)$ resonances. At this drive
amplitude, neither the chaotic layer nor the nonlinear resonances disrupt the
orbits associated with the qubit logical states and we do not expect
ionization to occur.

At larger drive amplitudes, the orbit of the nonlinear resonances in the
regular region shrinks towards the center of the main regular region, see Fig.
11(a). At the drive amplitude corresponding to $\bar{n}_{r}=185.6$, the
$(7\\!:\\!5)$ resonance collides with the Bohr-Sommerfeld orbit associated
with the qubit excited state, see Fig. 10(c) where the collision is
illustrated by coloring the $(7\\!:\\!5)$ resonance in the color of the
excited state. At that point, the qubit excited state can be unstable.
Importantly, as discussed further below, this does not necessarily result in
ionization. Indeed, weak resonances occupying a small area in phase space can
cross the qubit states without damage. By further increasing the drive
amplitude, the orbit of the resonance shrinks and no longer coincides with
that of the qubit state. Given the abundance of nonlinear resonances, the
pattern of resonances crossing and then moving away from the orbits of the
qubit ground and excited states is a common occurrence. This is the reason for
the jagged shape of the Christmas tree of Fig. 11(a) where large resonance
lead to large windows of drive amplitude for which these orbits cannot exist.
This observation is in agreement with the experimental results of Ref. Sank
_et al._ (2016) which showed that ‘quiet spots’ with minimal measurement-
induced transitions could be found at average resonator photon number larger
than the value where transitions are observed.

At a still larger drive amplitudes, initializing the classical system at the
energy corresponding to the qubit excited states does not lead to a regular
orbit but rather fills the chaotic layer, see Fig. 10(d) obtained for
$\bar{n}_{r}=287.6$. At that drive amplitude, the excited state is thus
engulfed by the chaotic layer and is unstable. At yet larger drive amplitudes
(not shown), the ground state eventually suffers the same fate and becomes
unstable. It is interesting to note the presence in panel (d) of a pair of
$(3\\!:\\!2)$ resonances near the inner edge of the chaotic layer. Because to
each resonance is associated a new separatrix around which chaos emerges, the
presence of those $(3\\!:\\!2)$ resonances further reduces the phase space
available to the orbit associated with the qubit excited state, thereby
hastening its merger with the chaotic layer. Note that the shrinking of the
resonances towards the Bohr-Sommerfeld orbits associated with the qubit states
is the classical analog of the quantum mechanism for ionization discussed in
Sec. II.3 and Sec. III.1, where high-energy states are ac Stark-shifted to
lower energies until they hybridize with the qubit states.

Building on these observations, Fig. 12(a-b) shows the critical photon numbers
$\\{n_{\mathrm{crit},i_{t}}^{c}\\}$ as a function of detuning $\Delta$. The
solid red lines indicate a collision of the orbit associated with the qubit
(a) ground or (b) excited state with a resonance, while the hatched red lines
rather indicate a collision with the chaotic layer. At
$\Delta/2\pi=-$1.202\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ these lines
correspond to the stars in Fig. 11(a). While chaos is expected to always
induce ionization, this is not the case for all collisions with resonances.
This is because many resonances occur at low photon numbers and are thus too
weak to induce ionization. To account for this, the width of the solid red
lines is proportional to the phase-space area of the resonance colliding with
the qubit state orbit, which is in turn proportional to the resonance’s
strength. Using this information, the olive lines indicate the smallest
critical photon number for which the phase-space area of the corresponding
resonance is larger than $0.05\times 2\pi\hbar_{\rm eff}$ in panel (a) and
$0.10\times 2\pi\hbar_{\rm eff}$ in panel (b). These values are free
parameters and were chosen to best fit the classical model to the Floquet
analysis, see the discussion below. The occasionally jagged behavior of that
line in panel (a) occurs when the area of the relevant resonances becomes
comparable to the resolution with which we compute phase-space areas.

These Bohr-Sommerfeld critical photon numbers are compared to the period-
averaged transmon excitation $\langle\langle\hat{N}_{t}\rangle\rangle_{i_{t}}$
(colored background), a quantity which was shown to be a good indicator of
transmon ionization in Sec. III. We have averaged the results over $n_{g}$ to
remove the dependence on this parameter. In both panels, we observe abrupt
changes in $\langle\langle\hat{N}_{t}\rangle\rangle_{i_{t}}$ near the photon
number where classical resonances occur. The agreement is excellent when
discarding weak resonances which are not expected to have an effect on the
qubit, see the olive lines. That some resonances are not relevant to transmon
ionization is consistent with the discussion in Sec. III.1, where it was
observed that weak resonances between Floquet modes do not cause ionization
because they can be traversed diabatically. Other correspondences between the
chaotic behavior of the driven classical pendulum and Floquet quasienergies
are discussed in Ref. Cohen  _et al._ (2023).

Figure 12: Bohr-Sommerfeld critical photon numbers
$\\{n_{\mathrm{crit},i_{t}}^{c}\\}$ (red lines) and gate charge average of
$\langle\langle\hat{N}_{t}\rangle\rangle_{i_{t}}$ (colored background) as a
function of detuning $\Delta$ and photon number $\bar{n}_{r}$. The top panel
corresponds to the qubit ground state for (a) negative detunings and (c)
positive detunings, and the bottom panel for the qubit excited state for (b)
negative detunings and (d) positive detunings. The hatched red lines indicate
the critical photon numbers due to the chaotic layer while the solid red lines
indicate the critical photon numbers due to a resonance. The width of the
solid red lines is proportional to the phase-space area of the resonance while
that of the hatched red lines was chosen to correspond to an area of
$2\pi\hbar_{\rm eff}$. The shaded red areas give the lower-bound uncertainty
due to quantum fluctuations around the Bohr-Sommerfeld orbit. The black dashed
lines are the standard Jaynes-Cummings critical photon number $n_{{\rm
crit},i_{t}}$ defined in Appendix A for each state. The hatched white regions
indicate regions where the dispersive approximation does not hold and for
which critical photon numbers were not calculated. In (a) and (b), the olive
lines indicate the smallest critical photon number for which the phase-space
area of the corresponding resonance is larger than $0.05\times 2\pi\hbar_{\rm
eff}$ and $0.10\times 2\pi\hbar_{\rm eff}$, respectively. For both classical
and Floquet theory, the vertical axis is produced by varying
$\varepsilon_{t}/2\pi$ in increments of
$10\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$.

###  IV.4 Positive detuning: $\omega_{q}>\omega_{r}$

We now turn to the case of positive detuning, see Fig. 10(e-h). As at negative
detuning, a chaotic layer develops when the drive is made finite, see panel
(f) obtained for $\bar{n}_{r}=0.63$ photons. The most dramatic feature is,
however, the appearance of a large $(1\\!:\\!1)$ resonance inside the
separatrix of $\tilde{H}_{t}$. In Fig. 11(b), the domed-shaped feature
centered at $\tilde{n}_{t}=0$ at low amplitude correspond to the main regular
region while the resonance is the Christmas tree-like feature centered at
$\tilde{n}_{t}\sim 1.6$. This resonance results from the negative
anharmoniticy of the transmon which at positive detuning allows for
$\tilde{\omega}(\tilde{I})=\tilde{\omega}_{d}$ to be satisfied for some
$\tilde{I}$. This is the classical equivalent of the breakdown of the
dispersive approximation at zero photon number due to a 1-photon qubit-
resonator resonance discussed in Sec. II.4.

While this resonance typically does not collide with the Bohr-Sommerfeld orbit
associated with the qubit ground and excited states, its presence rapidly
restricts the phase space available to the main regular region. This is
exacerbated by the presence of a new separatrix surrounding the $(1\\!:\\!1)$
resonance and about which a chaotic layer develops, see panel (g) obtained at
$\bar{n}_{r}=4.88$ photons and the V-shaped region between the main regular
region and the resonance in Fig. 11(b). Increasing further the drive
amplitude, the qubit excited state is engulfed by the chaotic layer and
becomes unstable, see panel (h) obtained at $\bar{n}_{r}=17.02$ photons. This
is followed by the ground state also entering the chaotic layer at still
larger photon number.

The area of the $(1\\!:\\!1)$ resonance is large enough ($>2\pi\hbar_{\rm
eff}$) to support its own Bohr-Sommerfeld orbits, see Appendix F. In Fig.
11(b,d), these orbits share the same color code as the distinct orbits of the
main regular region because they enclose the same areas. These orbits are
located in phase space within the separatrix of $\tilde{H}_{t}$, and thus have
an energy that is below the top of the cosine potential well of the transmon.
This is in direct analogy with branch bunching about a state $i_{t}^{\star}$
located below the top of the well discussed in Sects. II.4 and III.2.
Moreover, additional $(n\\!:\\!m)$ resonances appear within the main regular
region, see Fig. 11(b). However, these are small and are usually not relevant
compared to the large $(1\\!:\\!1)$ resonance discussed above. Only at
relatively large positive detunings does the critical photon number associated
with the $(1\\!:\\!1)$ resonance become large enough to allow other
$(n\\!:\\!m)$ resonances to become relevant for ionization (not shown).

Fig. 12(c-d) shows the critical photon numbers
$\\{n_{\mathrm{crit},i_{t}}^{c}\\}$ as a function of positive detuning
$\Delta$ and average photon number $\bar{n}_{r}$, with panel (c) corresponding
to the ground state and panel (d) to the excited state. Because the small
$(n\\!:\\!m)$ resonances crossing the qubit state orbits are of no
consequences at this range of detuning, the critical photon number here
corresponds to the point where the qubit states are engulfed by the chaotic
layer (hatched red lines). For most detunings, the presence of the large
$(1\\!:\\!1)$ resonance results in a much smaller critical photon number than
at negative detuning. Also in contrast with the case of negative detuning, the
Bohr-Sommerfeld critical photon number monotonically increases with the
magnitude of the detuning since the $(1\\!:\\!1)$ dominates the ionization
mechanism.

The critical photon number is compared with the gate-charge average of
$\langle\langle\hat{N}_{t}\rangle\rangle_{i_{t}}$ (colored background).
Although the classical predictions tend to slightly overestimate the critical
photon number obtained from Floquet theory, we again observe a good agreement
between the two approaches over a large range of detunings and photon numbers.
Some features that are not captured are the oscillations in the critical
photon number at small detunings, see Fig. 13 for a clearer view of the
oscillating features. As discussed in Sec. II.4 and Sec. III.2, these arise
from quantum tunneling physics and can thus not be captured within a classical
model. In the broad strokes, however, the classical mode of the driven
nonlinear pendulum captures remarkably well the essence of transmon
ionization.

##  V Comparison to experimental results

We now compare the critical photon numbers extracted from the quantum
transmon-resonator, Floquet, and classical models to experimental data from
Ref. Khezri  _et al._ (2023). This experiment consists in preparing the qubit
in its ground or first excited states, filling the resonator with
$\bar{n}_{r}$ photons for a controlled amount of time, and then performing a
QND readout of the qubit state. The colored background of Fig. 13 shows the
conditional probability $P(i_{t}|i_{t})$ for the qubit to be measured in the
state it was prepared in as a function of detuning and average resonator
photon number. Panel (a) is for the ground state and panel (b) for the excited
state. A white background indicates that the interaction with the
$\bar{n}_{r}$ photons did not cause qubit leakage, while a dark background
indicates leakage. In both cases, sharp dips are observed as a function of
detuning. In that experiment, the readout resonator has a bare frequency
$\omega_{r}/2\pi=$4.75\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$, the coupling
constant is $g/2\pi=$91\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$, the charging
energy is $E_{C}/2\pi=$194\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$, and the
qubit frequency $\omega_{q}/2\pi$ can be tuned from
$5.5\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ to
$6.34\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$. In all cases, this corresponds
to the regime of positive detuning.

In Ref. Khezri  _et al._ (2023), a semiclassical model predicting the position
and amplitudes of these dips is presented. Here we compare the experimental
results to the critical photon number obtained from our fully quantum
transmon-resonator (pink line) and classically-driven transmon (red line)
models. The agreement with experimental results is excellent for both models
and qubit initial states. Notably, the theoretical values capture the full
detuning dependence of the critical photon number including the presence of
dips in $\bar{n}_{r}$. Notably, it not only captures the position of these
dips but also their height and width. Those dips correspond to the multiphoton
resonances which are expected to precipitate ionization at positive detuning,
see Sec. II.4. Considering the experiment’s slow time scale relative to gate
charge fluctuations, the red and pink shaded areas—marking the critical photon
numbers’ distribution from the 10th to 90th percentiles over 100 samples of
$n_{g}$—emphasize the necessity of including gate charge in these predictions.
Moreover, as expected from the discussion in the previous sections, ionization
typically occurs earlier for the qubit excited state than the ground state.
Because it predicts the full dependence of the ionization threshold with
detuning, including regions where $\bar{n}_{{\rm crit},0_{t}}$ and
$\bar{n}_{{\rm crit},1_{t}}$ are large in between resonances, the framework
presented here could be exploited to maximize readout performance.

Figure 13: Critical photon numbers of our three models compared to the
experimental results of Ref. Khezri  _et al._ (2023) for (a) ground state and
(b) excited state. For a fixed detuning between the tunable transmon and
resonator frequencies, the experiment consists of: preparing the transmon
qubit in $\ket{0_{t}}$ or $\ket{1_{t}}$, driving the resonator at the
conditional dressed frequency $\omega_{r,0_{t}}$ and $\omega_{r,1_{t}}$ and
populating it with $\bar{n}_{r}$ photons on average, waiting for the resonator
to decay back to the vacuum, and finally reading out the qubit to extract the
conditional probabilities $P(i_{t}|i_{t})$ of finding the outcome $i_{t}$
given that $\ket{i_{t}}$ was initially prepared. The heat maps correspond to
the experimental results. The red line shows our prediction of $\bar{n}_{{\rm
crit},i_{t}}^{q}$ using the branch analysis of the full transmon-resonator
system as explained in previous sections. The pink line shows $\bar{n}_{{\rm
crit},i_{t}}^{f}$ obtained from the Floquet branch analysis. Shaded areas
correspond to the distribution of $n_{{\rm crit},i_{t}}^{q}$ and $n_{{\rm
crit},i_{t}}^{f}$ between the 10th and 90th percentiles over 100 realization
of $n_{g}$. Finally, the white full line shows the prediction from the
classical dynamics, with a lower bound given by white dashed line. The long
vertical features in panel (b) at very low photon number and
$\Delta/2\pi\gtrsim$1.25\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$ are
attributed to the presence of two-level defects Khezri  _et al._ (2023).

The prediction from the classical model (white full lines) follows the general
trend of critical photon numbers as a function of detuning. The white dashed
lines are a lower bound to these predictions accounting for quantum
fluctuations. Because the classical model does not capture the multiphoton
resonances responsible for accelerating ionization, the classical predictions
do not capture the dips in critical photon number versus detuning. Yet, it is
remarkable that as long as the full cosine potential is accounted for,
obtaining excellent estimates of the photon number threshold for ionization
matching experimental results can be obtained simply by solving Hamilton’s
equations of motion.

##  VI Conclusion and outlook

We have introduced a comprehensive framework providing a unified physical
picture of the origin of transmon ionization. Our study is based on three
complementary levels of description of the qubit-resonator system used for
dispersive qubit readout in circuit QED. A first approach treats the qubit and
the resonator as coupled quantum systems. A second approach simplifies the
situation by replacing the driven resonator by a classical drive acting
directly on the transmon. A third approach considers an entirely classical
problem. All three approaches have in common that they rely on the full cosine
potential of the tranmon’s Josephson junction. The first two descriptions
identify multiphoton resonances as the mechanism responsible for ionization.
The positions of these resonances can be easily located, and estimates of
critical photon numbers can be efficiently obtained numerically with both
methods leading to predictions which are in agreement with each other and with
recent experimental results. Those approaches also highlight the strong
dependence of the critical photon numbers on gate charge, even deep in the
transmon regime. The third level of description provides a simple
interpretation based on the motion of a classical driven pendulum and its
chaotic behavior. Although this method misses some features captured by the
quantum models, the critical photon numbers it predicts are in qualitative
agreement with the quantum predictions and with experimental results. It is
remarkable that classical dynamics of the driven transmon is sufficient to
capture the physics of ionization.

Although we have used throughout this work fixed sets of parameters, we have
verified that the conclusions we have reached remain general. While, e.g., the
specific multiphoton resonances precipitating ionization may shift with
changes in parameters, the dependence of the critical photon numbers
$\bar{n}_{{\rm crit},i_{t}}^{q}$ and $\bar{n}_{{\rm crit},i_{t}}^{f}$ on
detuning preserves the same features.

We hope that the analyses presented in this work will facilitate the
experimental exploration of parameter regions wherein transmon ionization
occurs at large photon numbers. For example, this could be done by exploiting
regions of detunings located between multiphoton resonances where the critical
photon number can be large. Furthermore, despite the reduced dispersive shift,
the use of large qubit-resonator detunings as a readout strategy seems
promising and has not yet been extensively investigated. The increase in
critical photon number with increasing absolute detuning $|\Delta|$ (see,
e.g., Fig. 7) indicates the potential of this approach for improving qubit
readout. Recent experimental findings lend support to this hypothesis
Kurilovich  _et al._ . Moreover, because the speed at which resonances are
traversed affects the probability of ionization occurring, optimizing the
readout pulse shape with an understanding of the position of resonance could
prove to be beneficial.

Beyond dispersive readout, many experiments in circuit QED contain the basic
ingredients that are at the core of ionization, namely strong drives and weak
anharmonicities. These include but are not limited to other types of qubits
Manucharyan  _et al._ (2009); Hyyppä  _et al._ (2022), readout schemes Didier
_et al._ (2015); Chapple  _et al._ , parametric gates Caldwell  _et al._
(2018); Niskanen  _et al._ (2007), and quantum state stabilization Puri  _et
al._ (2017); Frattini  _et al._ (2022); Chávez-Carlos  _et al._ (2023). The
tools introduced in this work can easily be applied to these scenarios. These
tools are simple to interpret, have an easy construction, are easy to
implement, and are computationally efficient. We thus anticipate the framework
introduced in this work to become an integral part of the theoretical analysis
and support of current and future circuit QED experiments.

## Acknowledgements

The authors are grateful to Ross Shillito for helpful discussions. This work
is supported by a collaboration between the US DOE and other Agencies. This
material is based upon work supported by the U.S. Department of Energy, Office
of Science, National Quantum Information Science Research Centers, Quantum
Systems Accelerator. Additional support is acknowledged from NSERC, the Canada
First Research Excellence Fund, the Ministère de l’Économie et de l’Innovation
du Québec, and Fonds de recherche du Québec - Nature et technologies.

##  Appendix A Dispersive critical photon number

In the main text, we compared the various critical photon numbers defined in
our work with the critical photon number $n_{{\rm crit},i_{t}}$ that defines
the dispersive limit. Here, we summarize how these dispersive critical photon
numbers are calculated.

The commonly used Jaynes-Cummings critical photon number for the qubit
transition is defined as $n_{\rm
crit}=|(\omega_{0_{t}1_{t}}-\omega_{r})/2g_{0_{t}1_{t}}|^{2}$ Blais  _et al._
(2004, 2021). This is the photon number for which the interaction term in Eq.
1 becomes resonant with the qubit transition. However, this definition is
insufficient to capture the breakdown of the dispersive regime in a multilevel
system like the transmon. It is instead necessary to define a critical photon
number for every transition D’Anjou and Burkard (2019). To do so, we simply
transpose the definition for the qubit transition to all other transitions,

$\displaystyle\begin{aligned} n_{{\rm
crit},k_{t}l_{t}}\equiv\left|\frac{\omega_{k_{t}l_{t}}-\omega_{r}}{2g_{k_{t}l_{t}}}\right|^{2}.\end{aligned}$
(15)

Note that this definition naturally includes counter-rotating transitions as
the case where $\omega_{k_{t}l_{t}}\sim-\omega_{r}$. We define the dispersive
critical photon number for an initial transmon state $|i_{t}\rangle$ as the
smallest of the $n_{{\rm crit},k_{t}l_{t}}$ that involves the initial state:

$\displaystyle n_{{\rm crit},i_{t}}=\min_{k_{t}l_{t}}\\{n_{{\rm
crit},k_{t}l_{t}}|k_{t}=i_{t}\text{ or }l_{t}=i_{t}\\}.$ (16)

Due to the dipole-like selection rules of the transmon regime, the limiting
transitions for the qubit states are always $0_{t}\leftrightarrow 1_{t}$ and
$1_{t}\leftrightarrow 2_{t}$. As a result, $n_{{\rm crit},0_{t}}$ and $n_{{\rm
crit},1_{t}}$ are essentially independent of the gate charge $n_{g}$. Equation
15 gives similar predictions as other definitions of the state-dependent
critical photon numbers for the transmon Blais  _et al._ (2021), but has the
advantage of being applicable to any multilevel system.

##  Appendix B Nonperturbative corrections to the branches at negative
detuning

In Sec. II.3 we argued that nonperturbative corrections to the approximate
eigenstates obtained via Schrieffer-Wolff perturbation theory
$|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle i_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle i_{t},n_{r}$\cr} }}}\rangle\equiv
e^{-\hat{S}}|i_{t},n_{r}\rangle\approx|\overline{i_{t},n_{r}}\rangle$ was
ultimately the cause of branch swapping. Here we discuss at length the
mechanism through which these nonperturbative corrections occur.

As an illustrative example, let us consider the swapping of $B_{1_{t}}$ and
$B_{7_{t}}$ which occurs at roughly 84 photons, see Fig. 2. If we were only to
consider the bare transmon energies, swapping between these branches would be
unexpected. Indeed, even if one considers the possibility of a multiphoton
process connecting such states, the detuning between them is still large
$(E_{7_{t}}-E_{1_{t}}-4\omega_{r})/2\pi\approx$1.98\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$\gg
g/2\pi$. Moreover, the interaction cannot directly add or remove 4 photons
from the resonator, and any such process that appears within perturbation
theory should be small. It is thus at this stage not evident why the swapping
should occur between these branches.

To address how the requisite term mediating this process arises, we must keep
more terms in the infinite-order expansion than in Eq. 6. To that end, we
recall that the standard approach Koch  _et al._ (2007); Schrieffer and Wolff
(1966) dictates that by choosing

$\displaystyle\hat{S}=-i(\hat{S}_{t}\hat{a}^{\dagger}+\hat{a}\hat{S}_{t})$
(17)

with

$\displaystyle\hat{S}_{t}=\sum_{i_{t},j_{t}=0}^{\infty}\frac{g_{i_{t}j_{t}}}{\omega_{i_{t}j_{t}}-\omega_{r}}|i_{t}\rangle\langle
j_{t}|$ (18)

we obtain

$\displaystyle
e^{\hat{S}}\hat{H}_{tr}e^{-\hat{S}}=\omega_{r}\hat{a}^{\dagger}\hat{a}+\hat{H}_{t}+\hat{a}^{\dagger}\hat{a}\sum_{i_{t},j_{t}=0}^{\infty}\chi_{i_{t}j_{t}}|i_{t}\rangle\langle
j_{t}|$
$\displaystyle+\left(\hat{a}^{\dagger}\hat{a}^{\dagger}\sum_{i_{t},j_{t}=0}^{\infty}\eta_{i_{t}j_{t}}|i_{t}\rangle\langle
j_{t}|+\mathrm{h.c.}\right)+\sum_{i_{t},j_{t}=0}^{\infty}\Lambda_{i_{t}j_{t}}|i_{t}\rangle\langle
j_{t}|$ $\displaystyle+O(g^{3}),$ (19)

where the matrix matrix elements are

$\displaystyle\chi_{i_{t}j_{t}}$
$\displaystyle=\sum_{k_{t}=0}^{\infty}g_{i_{t}k_{t}}g_{k_{t}j_{t}}\left(\frac{\omega_{i_{t}k_{t}}}{\omega_{r}^{2}-(\omega_{i_{t}k_{t}})^{2}}+\frac{\omega_{j_{t}k_{t}}}{\omega_{r}^{2}-(\omega_{j_{t}k_{t}})^{2}}\right),$
(20) $\displaystyle\eta_{i_{t}j_{t}}$
$\displaystyle=\frac{1}{2}\sum_{k_{t}=0}^{\infty}g_{i_{t}k_{t}}g_{k_{t}j_{t}}\left(\frac{1}{\omega_{r}-\omega_{k_{t}j_{t}}}-\frac{1}{\omega_{r}+\omega_{k_{t}i_{t}}}\right),$
(21) $\displaystyle\Lambda_{i_{t}j_{t}}$
$\displaystyle=-\frac{1}{2}\sum_{k_{t}=0}^{\infty}g_{i_{t}k_{t}}g_{k_{t}j_{t}}\left(\frac{1}{\omega_{r}-\omega_{k_{t}j_{t}}}+\frac{1}{\omega_{r}-\omega_{k_{t}i_{t}}}\right),$
(22)

and we have dropped all terms of order $g^{3}$ and higher. Note that we use a
different notation than that presented in Ref. Koch  _et al._ (2007). We note
that Xiao  _et al._ (2023) introduce a diagramatic approach to obtain high-
order effective Hamiltonian for a driven transmon in the absence of a
resonance.

By keeping the diagonal terms in Eq. 19, we recover Eq. 6 and thus correctly
obtain the dispersion to order $g^{2}$. A salient feature of Eq. 19, however,
is that to that same order in the coupling we have terms that create pairs of
photons in the resonator $\hat{a}^{\dagger}\hat{a}^{\dagger}$. These terms can
also cause transitions between any transmon states with an amplitude
$\eta_{i_{t}j_{t}}$. If we were to explicitly keep track of terms to higher-
order in $g$, we would similarly conclude that to fourth order there appears
an analogous term
$\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}^{\dagger}$.
Nevertheless, it would seem that this term is irrelevant, considering that it
is of fourth order in a nominally small parameter.

As we emphasize in the main text, however, in this new basis the partially-
diagonalized eigenstates are ac-Stark shifted. At large enough photon number,
the $\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}^{\dagger}\hat{a}^{\dagger}$
term that appears in the new basis can couple dressed states which are nearly
degenerate. If the condition

$\displaystyle\frac{|\langle\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 1_{t},n_{r}+4$\cr}
}}}|\hat{H}_{tr}|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 7_{t},n_{r}$\cr}
}}}\rangle|}{\left|\langle\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 7_{t},n_{r}$\cr}
}}}|\hat{H}_{tr}|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 7_{t},n_{r}$\cr}
}}}\rangle-\langle\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 1_{t},n_{r}+4$\cr}
}}}|\hat{H}_{tr}|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 1_{t},n_{r}+4$\cr} }}}\rangle\right|}\gg 1$ (71)

is satisfied for some $n_{r}$, the eigenstates are a highly-entangled
superposition of $|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 1_{t},n_{r}+4$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 1_{t},n_{r}+4$\cr} }}}\rangle$ and
$|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 7_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 7_{t},n_{r}$\cr} }}}\rangle$ and there is an
avoided crossing in the spectrum. Before this occurs, the branch analysis
classifies $|\overline{1_{t},n_{r}}\rangle\approx|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=2.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\displaystyle 1_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.5pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\textstyle 1_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.25pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptstyle 1_{t},n_{r}$\cr} }}{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule
height=0.4pt,width=1.0pt\kern 0.5pt}\hfill\kern-0.5pt\cr\kern
1.0pt\cr$\scriptscriptstyle 1_{t},n_{r}$\cr} }}}\rangle$ and
$|\overline{7_{t},n_{r}}\rangle\approx|\mathop{\mathchoice{\vbox{
\offinterlineskip\halign{#\cr\kern-0.5pt\xleaders\hbox{\kern 0.5pt\vrule