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$, but after the
crossing the characters of the states have been switched and we have
$|\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 7_{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 7_{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 7_{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 7_{t},n_{r}-4$\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
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$.

We can confirm this interpretation to be valid by plotting the overlaps of the
approximate eigenstates with the relevant numerically-obtained eigenstates as
in Fig. 14. We go to one order higher in perturbation theory than in Eq. 19 by
defining

$\displaystyle|\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}}e^{-\hat{T}}|i_{t},n_{r}\rangle,$ (80)

where $\hat{T}$ is proportional to $g^{2}$ chosen to eliminate the off-
diagonal terms to that order in Eq. 19. Its explicit form is simple to obtain
and is merely an exercise in commutator algebra and thus we do not reproduce
it here. At low photon number, we see that
$|\langle\overline{1_{t},n_{r}}|\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|^{2}\approx|\langle\overline{7_{t},n_{r}}|\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|^{2}\approx 1$,
indicating the validity of the dispersive approximation. Near 84 photons, we
see these overlaps experience a sharp drop, after which they become
effectively zero. Instead, the overlaps
$|\langle\overline{1_{t},n_{r}}|\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}-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 7_{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 7_{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 7_{t},n_{r}-4$\cr} }}}\rangle|^{2}$ and
$|\langle\overline{7_{t},n_{r}}|\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|^{2}$ rise and
become finite, indicating branch swapping. This makes it clear that near these
photon numbers, these two approximate eigenstates are strongly mixed and the
condition of Eq. 71 is satisfied. The fact that the overlaps do not remain
close to unity on each side of the branch swapping is a consequence of using a
finite-order expansion for the partially-dressed eigenstates. Indeed, such an
expansion becomes less accurate at large $n_{r}$ since the transmon-resonator
interaction strength scales with photon number. They could be made arbitrarily
close to unity before and after the crossing by going to even higher order in
perturbation theory. At the crossing itself, however, the observed branch
swapping is always nonperturbative in the coupling $g$.

We emphasize that we focused on $B_{1_{t}}$ and $B_{7_{t}}$ since the branch
analysis had already identified that these two had swapped. Without the branch
analysis, we would need to verify the resonance condition, Eq. 71, for every
pair of branches. Further, we would need to ensure that the perturbative
energies $\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 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}
}}}|\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 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$ are accurate. At
higher photon number, this means going to higher order in perturbation theory,
a needlessly tedious task given that the branch analysis already contains this
information.

Finally, let us stress that the order of the multiphoton process responsible
for the swapping of the ground or excited state branch with any other can be
diagnosed by considering the modular spectrum of the bare transmon states. For
instance, here we knew that the four-photon process was responsible for the
swapping of the seventh and first branch since
$(E_{7_{t}}-E_{1_{t}}-4\omega_{r})/2\pi=$1.98\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$$
was the smallest detuning between these transmon states and any other
multiphoton process. The order of these multiphoton resonances thus changes
with the detuning $\Delta$. For instance, going from small to large negative
detuning, the three peaks we see in the detuning sweep of Fig. 7(a) are a
consequence of the multiphoton process of order 4, 3, and 2 between the ground
state and branches $B_{7_{t}}$, $B_{6_{t}}$ and $B_{5_{t}}$, respectively.

Figure 14: Overlaps between the eigenstates of branches $B_{1_{t}}$ and
$B_{7_{t}}$ and the corresponding perturbative approximations
$|\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$ obtained via
Schrieffer-Wolff perturbation theory, Eq. 80. Before the branch swapping at
approximately 84 photons, the eigenstates in each branch are well-approximated
by their perturbative counterpart. Near the swapping, perturbation theory
fails and the eigenstates are in a roughly equal superposition of the
perturbative eigenstates. For larger photon number, the character of the
eigenstates has switched. Deviation of the overlaps from unity beyond the
swapping point is a consequence of the finite-order approximation for the
$|\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$, which fails at large
photon numbers.

##  Appendix C Effective drive on the transmon

In Sec. III, we obtain Eq. 10 by applying a displacement transformation
$\hat{D}(\alpha)=\exp(\alpha\hat{a}^{\dagger}-\alpha\hat{a})$ on the driven
qubit-resonator master equation Eq. 3. In that transformation, the amplitude
$\alpha$ is time dependent and takes the form
$\alpha(t)=\int_{0}^{t}d\tau\,\mathcal{E}_{d}(\tau)e^{-(i\omega_{r}+\kappa/2)(t-\tau)}$.
With that choice, the transformation eliminates the resonator drive from the
Hamiltonian and induces an effective transmon drive
$\mathcal{E}_{t}(t)\hat{n}_{t}=2g\text{Im}[\alpha(t)]\hat{n}_{t}$.

Taking the resonator drive to be
$\mathcal{E}_{d}(t)=\varepsilon_{d}\sin(\omega_{d}t)$ with constant amplitude
$\varepsilon_{d}$, we have that

$\begin{split}\alpha(t)=&\frac{\varepsilon_{d}}{2}\left[\frac{e^{-i\omega_{d}t}-e^{-(i\omega_{r}+\kappa/2)t}}{\omega_{-}-i\kappa/2}\right.\\\
&\qquad-\left.\frac{e^{i\omega_{d}t}-e^{-(i\omega_{r}+\kappa/2)t}}{\omega_{+}-i\kappa/2}\right],\end{split}$
(81)

where $\omega_{\pm}=\omega_{r}\pm\omega_{d}$. The amplitude of the transmon
drive therefore reads

$\begin{split}\mathcal{E}_{t}(t)=&g\varepsilon_{d}\left[\frac{\cos(\omega_{d}t-\phi_{-})-\cos(\omega_{r}t-\phi_{-})e^{-\kappa
t/2}}{\sqrt{\omega_{-}^{2}+(\kappa/2)^{2}}}\right.\\\
&\quad-\left.\frac{\cos(\omega_{d}t+\phi_{+})-\cos(\omega_{r}t-\phi_{+})e^{-\kappa
t/2}}{\sqrt{\omega_{+}^{2}+(\kappa/2)^{2}}}\right],\end{split}$ (82)

where $\tan\phi_{\pm}=-2\omega_{\pm}/\kappa$. Because $\omega_{+}\gg
g\varepsilon_{d}$, we can ignore the second line of the above expression.
Furthermore, assuming that $|\omega_{-}|\ll\kappa$, the drive amplitude takes
the simpler form

$\mathcal{E}_{t}(t)=\frac{2g\varepsilon_{d}}{\kappa}(1-e^{-\kappa
t/2})\cos(\omega_{d}t).$ (83)

This is the expression that we use in Sec. III.

Beyond this limiting regime, the effective drive on the transmon Eq. 82 can
still be well approximated by a fast periodic drive with an envelope that
grows slowly as $(1-e^{-\kappa t/2})$, although not monotonically due to a
finite difference $\omega_{r}-\omega_{d}$. In particular, for large
$|\omega_{r}-\omega_{d}|$ compared to $\kappa$, an avoided crossing in the
Floquet spectrum might be dynamically crossed multiple times.

An important remark is that, when $|\omega_{r}-\omega_{d}|\gtrsim\kappa$, a
better approximation to the full transmon-resonator dynamics is still obtained
by using Eq. 83 rather than Eq. 82 for just a driven transmon. The reason is
that the resonator responds at a qubit-state-dependent frequency that is
typically closer to the readout frequency $\omega_{d}$ than to $\omega_{r}$,
thus the monochromatic-tone approximation is more accurate.

##  Appendix D Floquet formalism

While for the numerical simulations we extract the Floquet Hamiltonian
directly from the evolution operator over one period, $\hat{U}(T)\equiv
e^{-i\hat{H}_{F}T}$, it is instructive to express the Floquet problem in
Shirley space, also known as the infinite replicas picture Shirley (1965).

The replicas picture is obtained by solving the time-dependent Schrödinger
equation in frequency space. Defining the Floquet quasienergies and modes of
the time-dependent periodic Floquet Hamiltonian via
$\hat{H}_{F}\ket{\phi_{i_{t}}}=\epsilon_{i_{t}}\ket{\phi_{i_{t}}}$, a
particular solution to the Schrödinger equation is the Floquet state
$\ket{\psi_{i_{t}}(t)}=e^{-i\epsilon_{i_{t}}t}\ket{\phi_{i_{t}}(t)}$. The
quasienergy is defined within the first Brillouin zone
$\epsilon_{i_{t}}\in[-\omega_{d}/2,\omega_{d}/2]$. The Floquet modes are
periodic in time, $\ket{\phi_{i_{t}}(t)}=\ket{\phi_{i_{t}}(t+T)}$, where
$\ket{\phi_{i_{t}}(t)}=\hat{U}(t)e^{it\hat{H}_{F}}\ket{\phi_{i_{t}}(0)}$ and
$\ket{\phi_{i_{t}}}\equiv\ket{\phi_{i_{t}}(0)}$.

In the case of the driven transmon, expanding the Hamiltonian $\hat{H}_{t}(t)$
of Eq. 10 and the Floquet modes in Fourier components as
$\hat{H}_{t}(t)=\hat{H}_{t}+(\varepsilon_{t}/2)(e^{i\omega_{d}t}+e^{-i\omega_{d}t})\hat{n}_{t}$
and $\ket{\phi_{i_{t}}(t)}=\sum_{m}e^{im\omega_{d}t}\ket{\phi_{i,m}}$ and
substituting them in Schröginger’s equation
$i\partial_{t}\ket{\psi_{t}(t)}=\hat{H}(t)\ket{\psi_{t}(t)}$, leads to the
eigenvalue problem

$\begin{split}&(\hat{H}_{t}+m\omega_{d})\ket{\phi_{i_{t},m}}+\frac{\varepsilon_{t}}{2}\hat{n}_{t}(\ket{\phi_{i_{t},m+1}}+\ket{\phi_{i_{t},m-1}})\\\
&=\epsilon_{i_{t}}\ket{\phi_{i_{t},m}}.\end{split}$ (84)

The infinite replicas eigenvalue problem written in Eq. 84 can be interpreted
more clearly if written in matrix form as

$\begin{split}&\begin{pmatrix}\ddots&\ddots&\ddots&\\\
\ddots&\hat{H}_{t}-\omega_{d}&\frac{\varepsilon_{t}}{2}\hat{n}_{t}&0\\\
\ddots&\frac{\varepsilon_{t}}{2}\hat{n}_{t}&\hat{H}_{t}&\frac{\varepsilon_{t}}{2}\hat{n}_{t}&\ddots\\\
&0&\frac{\varepsilon_{t}}{2}\hat{n}_{t}&\hat{H}_{t}+\omega_{d}&\frac{\varepsilon_{t}}{2}\hat{n}_{t}\\\
&&\ddots&\ddots&\ddots&\ddots\end{pmatrix}\begin{pmatrix}\vdots\\\
\ket{\phi_{i_{t},-1}}\\\ \ket{\phi_{i_{t},0}}\\\ \ket{\phi_{i_{t},1}}\\\
\vdots\end{pmatrix}\\\ &\qquad=\epsilon_{i_{t}}\begin{pmatrix}\vdots\\\
\ket{\phi_{i_{t},-1}}\\\ \ket{\phi_{i_{t},0}}\\\ \ket{\phi_{i_{t},1}}\\\
\vdots\end{pmatrix}.\end{split}$ (85)

The original time-dependent eigenvalue problem has thus been mapped to a
static eigenvalue problem involving infinite replicas of transmons ordered in
a one-dimensional lattice. The transmon at site $m$ has a static Hamiltonian
$\hat{H}_{t}+m\omega_{d}$ and is coupled to nearest neighbors $(m\pm 1)$ sites
via the operator $(\varepsilon_{t}/2)\hat{n}_{t}$. From that point of view,
this “tunneling” operator couples transmon states belonging to neighboring
replica sites. In the absence of a drive $\varepsilon_{t}=0$, the Floquet
mode’s component $m$ is just
$\ket{\phi_{i_{t},m}[\varepsilon_{t}=0]}=\ket{i_{t},m}$, i.e., the bare
transmon state $i_{t}$ at site $m$. Moreover, in the replicas picture, the
bare states $\ket{i_{t},m}$ and $\ket{j_{t},n}$ can hybridize at finite
$\varepsilon_{t}$ provided that these states can be connected via a tunneling
path along the one-dimensional lattice. Moving back to the original time-
dependent picture, the interpretation is that there can be a transition
between levels $i_{t}$ and $j_{t}$ provided $n-m$ drive photons at frequency
$\omega_{d}$ are absorbed ($n>m$) or emitted ($n<m$).

With a second-order perturbative expansion in $\varepsilon_{t}$, we can
compute a correction to the bare ($\varepsilon_{t}=0$) frequencies
$E_{i_{t}}+m\omega_{d}$. The ac Stark shift of level $i_{t}$ is

$\begin{split}\chi_{i_{t},m}=&\frac{\varepsilon_{t}^{2}}{4}\sum_{j_{t}}\left(\frac{|\bra{i_{t}}\hat{n}_{t}\ket{j_{t}}|^{2}}{E_{i_{t},m}-E_{j_{t},m+1}}+\frac{|\bra{i_{t}}\hat{n}_{t}\ket{j_{t}}|^{2}}{E_{i_{t},m}-E_{j_{t},m-1}}\right)\\\
=&\frac{\varepsilon_{t}^{2}}{2}\sum_{j_{t}}\frac{\omega_{i_{t}j_{t}}|\bra{i_{t}}\hat{n}_{t}\ket{j_{t}}|^{2}}{\omega_{d}^{2}-\omega_{i_{t}j_{t}}^{2}}\\\
=&2\bar{n}_{r}\sum_{j_{t}}\frac{\omega_{i_{t}j_{t}}|g_{i_{t}j_{t}}|^{2}}{\omega_{d}^{2}-\omega_{i_{t}j_{t}}^{2}},\end{split}$
(86)

where we have used the definitions $\omega_{i_{t}j_{t}}=E_{j_{t}}-E_{i_{t}}$,
$g_{i_{t}j_{t}}=g\bra{i_{t}}\hat{n}_{t}\ket{j_{t}}$, and
$\varepsilon_{t}=2g\sqrt{\bar{n}_{r}}$. Note that this $\chi_{i_{t},m}$ shift
is independent of the replica index $m$, as expected as all replicas are
equivalent. Importantly, this expression is exactly the frequency shift of the
transmon level $i_{t}$ with an average number of photons $\bar{n}_{r}$ in the
resonator and for $\omega_{d}=\omega_{r}$ that is expected from Eq. 7 obtained
with the full quantum model.

##  Appendix E Floquet-Landau-Zener transition probability

The Landau-Zener formula for the transition probability at an anticrossing
between two Floquet quasienergies gives accurate predictions, even for fast
pulses Drese and Holthaus (1999); Ikeda  _et al._ (2022). In our mapping of
the driven transmon-resonator Hamiltonian to a Floquet problem, we do not
directly change the detuning of the quasienergies linearly in time as would be
the case for the standard Landau-Zener argument. Instead, the relevant
parameters are the amplitude of the drive on the resonator $\varepsilon_{d}$
and the resonator decay rate $\kappa$, which in turn changes the effective
drive amplitude on the transmon following Eq. 83. The latter is the quantity
which controls the detuning between the quasienergies.

To deduce the effective speed $v$ at which resonances are traversed as a
function of the drive amplitude $\varepsilon_{t}(t)$, we follow Ref. Ikeda
_et al._ (2022) in approximating the difference in the quasienergies close to
the avoided crossing as

$\begin{split}\epsilon_{j_{t}}[\varepsilon_{t}(t)]-\epsilon_{i_{t}}[\varepsilon_{t}(t)]&\approx\sqrt{\Delta_{\rm
ac}^{2}+[v(t-t_{\rm ac})]^{2}}\\\ &\approx\Delta_{\rm
ac}+\frac{1}{2}\frac{[v(t-t_{\rm ac})]^{2}}{\Delta_{\rm ac}}.\end{split}$ (87)

In this expression, $t_{\rm ac}$ is the time at which the minimum gap
$\Delta_{\rm ac}$ in the avoided crossing is reached. The value of the drive
amplitude at that time is denoted $\varepsilon_{t,{\rm
ac}}\equiv\varepsilon_{t}(t_{\rm ac})$. Expanding the left-hand-side of Eq. 87
to second order (the first derivatives vanish at the avoided crossing), we
obtain

$\begin{split}&\epsilon_{j_{t}}[\varepsilon_{t}(t)]-\epsilon_{i_{t}}[\varepsilon_{t}(t)]\approx\epsilon_{j_{t}}[\varepsilon_{t,{\rm
ac}}]-\epsilon_{i_{t}}[\varepsilon_{t,{\rm ac}}]\\\
&+\bigg{|}\frac{d^{2}\epsilon_{j_{t}}[\varepsilon_{t}]}{d\varepsilon_{t}^{2}}\bigg{|}_{\varepsilon_{t,{\rm
ac}}}\left(\frac{d\varepsilon_{t}(t)}{dt}\bigg{|}_{t_{\rm
ac}}\right)^{2}(t-t_{\rm ac})^{2}.\end{split}$ (88)

Combining the two above expressions, we find that the speed $v$ is given by
the relation

$\frac{v^{2}}{2\Delta_{\rm
ac}}=\bigg{|}\frac{d^{2}\epsilon_{j_{t}}[\varepsilon_{t}]}{d\varepsilon_{t}^{2}}\bigg{|}_{\varepsilon_{t,{\rm
ac}}}\left(\frac{d\varepsilon_{t}(t)}{dt}\bigg{|}_{t_{\rm ac}}\right)^{2}.$
(89)

Note that we assumed here that the second derivatives of
$\epsilon_{i_{t}}[\varepsilon_{t}]$ and $\epsilon_{j_{t}}[\varepsilon_{t}]$
are equal at the avoided crossing.

Knowing the gap size $\Delta_{\rm ac}$ and the speed $v$, the Landau-Zener
diabatic transition probability can be computed as $P_{\rm
LZ}=\exp[-\pi\Delta_{\rm ac}^{2}/2v]$. For the parameters used in Fig. 8(c),
we obtain $P_{\rm LZ}=0.35$ for the diabatic crossing, and $1-P_{\rm LZ}=0.65$
for the adiabatic crossing. The quasienergy gap and the speed are numerically
obtained from the tracked quasienergies and the functional form of the drive
amplitude $\varepsilon_{t}(t)$.

##  Appendix F Comparison of the Bohr-Sommerfeld energies with Floquet theory

In this appendix, we show that the Bohr-Sommerfeld quantization procedure
gives a sensible approximation to quantum theory for standard transmon
parameters. This is done by comparing the results obtained from the Bohr-
Sommerfeld quantization approach to the Floquet analysis developed in Sec. III
for $n_{g}=0$. More precisely, we compare the time-averaged transmon energy
over the coordinates $\\{\varphi_{i_{t}}(\tau),n_{i_{t}}(\tau)\\}$ for the
Bohr-Sommerfeld orbits,

$\langle\langle
H_{t}\rangle\rangle_{i_{t}}=\lim_{t\rightarrow\infty}\frac{1}{t}\int_{0}^{t}d\tau\,H_{t}\bigl{(}\varphi_{i_{t}}(\tau),n_{i_{t}}(\tau)\bigr{)},$
(90)

to the time-averaged energy of the Floquet modes Cohen  _et al._ (2023),

$\langle\langle\hat{H}_{t}\rangle\rangle_{i_{t}}=\frac{1}{T}\int_{0}^{T}d\tau\,\bra{\phi_{i_{t}}(\tau)}\hat{H}_{t}\ket{\phi_{i_{t}}(\tau)}.$
(91)

In these plots, the gray symbols correspond to the Floquet modes energy given
by Eq. 91, the colored symbols to the time-averaged energies of the Bohr-
Sommerfeld orbits given by Eq. 90. For the latter, we use the same color code
as in Fig. 11. The scale of the vertical axis is chosen such that zero energy
corresponds to the bottom of the cosine well while $2E_{J}$ to the energy of
the top of the well, i.e. the energy of the separatrix of $H_{t}$. The
significance of the colored full and dashed lines will be explained below.

Figure 15:  Time-averaged energies $\langle\langle
H_{t}\rangle\rangle_{i_{t}}$ of the Bohr-Sommerfeld orbits (colored dots with
color code of Fig. 11) and time-averaged energies
$\langle\langle\hat{H}_{t}\rangle\rangle_{i_{t}}$ of the Floquet modes at
$n_{g}=0$ (gray dots) as a function of the rescaled drive amplitude
$\tilde{\varepsilon}_{t}$. The solid lines indicate the energies of the center
of the main regular region (pink) and the $(1\\!:\\!1)$ resonance (indigo).
The dashed lines indicate the energy of the boundaries of the regular regions,
namely the main regular region (pink), the $(1\\!:\\!1)$ resonance (indigo),
and the innermost and outermost orbits external to the $(1\\!:\\!1)$ resonance
at positive detuning (green). We show (a) the case of negative detuning
($\tilde{\omega}_{d}=1.152$) and (b) the case of positive detuning
($\tilde{\omega}_{d}=0.807$). The Bohr-Sommerfeld $\langle\langle
H_{t}\rangle\rangle_{i_{t}}$ are obtained by averaging the Hamiltonian over
$\mathcal{N}=1000$ periods using $40$ points per period. The Floquet
$\langle\langle\hat{H}_{t}\rangle\rangle_{i_{t}}$ are obtained using $31$
points per period. The drive amplitude is varied in increments of
$\delta\varepsilon_{t}/2\pi=\omega_{p}\delta\tilde{\varepsilon}_{t}/2\pi=$10\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$.
The Bohr-Sommerfeld orbits and the region boundaries were identified using the
procedure described in Appendix G.

Focusing first on the case of negative detuning, we first observe that the
average energy of the Floquet modes (gray symbols) located near the bottom of
the well change smoothly with drive amplitude at small drive amplitude. On the
other hand, modes located near the separatrix rapidly hybridize resulting in
erratic behavior with the drive amplitude. In Ref. Cohen  _et al._ (2023),
this was shown to be a signature of chaotic behavior. As the drive amplitude
is further increased, the chaotic layer widens until it absorbs the Floquet
modes associated with the qubit logical states. Note that the resonances
responsible for ionization can be seen as sharp features of the average energy
connecting Floquet modes at the bottom of the well with modes close to the top
of the well.

Before chaotic behavior emerges, the energy of the Bohr-Sommerfeld orbits
(colored symbols) match very well the average Floquet energies. It is
interesting to note that this is the case even in this situation where there
are few levels in the well. The classical method is thus quite robust.
Moreover, the discontinuities that are observed in the energies of the Bohr-
Sommerfeld orbits coincide with the ranges of $\tilde{\varepsilon}_{t}$ for
which the corresponding orbits cease to exist because of collisions with
$(n\\!:\\!m)$ resonances. The full and dashed pink lines denote the average
energies of the orbit at the center of the main regular region [corresponding
to the gray line in Fig. 11(a)] and the average energy of the outermost orbit
of the main regular region, respectively. These two lines delimit the energy
window into which regular orbits exist, i.e., where Bohr-Sommerfeld states are
supported. The correspondence between the pink dashed line and the region
where the average Floquet energy shows chaotic behavior is striking.

At positive detunings we again see an excellent agreement between the
classical and the Floquet predictions at small drive amplitudes, see Fig.
15(b). Notably, both approaches show the emergence of a set of states in the
vicinity of $i_{t}=4$ whose energy changes smoothly with drive amplitude. This
is the bunching layer caused by a $(1\\!:\\!1)$ resonance discussed in Sec.
II.4 and Sec. III.2. Chaotic behavior emerges not only close to the separatrix
of $\tilde{H}_{t}$ but also close to the bunching layer. This is the behavior
expected from Figs. 10 and 11 and results from the apperance of a new
separatix close to the resonance. Increasing further the drive amplitude, the
chaotic layer widens until all states inside the cosine potential become
chaotic.

Because of the presence of the large $(1\\!:\\!1)$ resonance, there are three
different regions that can support Bohr-Sommerfeld states, namely the main
regular region, the $(1\\!:\\!1)$ resonance itself, and the orbits surrounding
both the main regular region and the $(1\\!:\\!1)$ resonance, see Fig. 11(b).
In Fig. 15(b), the full and dashed pink lines correspond to the average energy
of the orbits at the center and inner edge of the main regular region,
respectively. Furthermore, the full and dashed indigo lines are the energy of
the center and inner edge of the $(1\\!:\\!1)$ resonance. The lower and upper
dashed green lines are the average energy of the innermost and outermost
orbits external to the $(1\\!:\\!1)$ resonance. The agreement between the
Bohr-Sommerfeld energies and Floquet modes energies in the three identified
zones is excellent. We have verified that the correspondence between the
bunching layer Floquet states and the $(1\\!:\\!1)$ resonance holds better as
$\hbar_{\rm{eff}}$ is made smaller (not shown). In the limit
$\hbar_{\rm{eff}}\rightarrow 0$ ($E_{J}/E_{C}\rightarrow\infty$), the lowest
energy Bohr-Sommerfeld state coincides with the energy of the center of the
main regular region and ionizes at around $\tilde{\varepsilon}_{t}\sim 0.3$
where the pink full and dashed lines meet in Fig. 15. Consequently, the Bohr-
Sommerfeld critical photon number associated to the $(1\\!:\\!1)$ resonance
has a saturating value.

Finally, we note that the center of the $(1\\!:\\!1)$ resonance is a point of
stability with higher time-averaged energy than the surrounding points. Thus,
the bunched states are excitations within an effective inverted potential well
generated by the drive. Such an inverted effective potential well bears a
superficial resemblance to the Kapitza pendulum, i.e. a rigid pendulum
stabilized ‘upside-down’ by periodically changing the height of the pivot
point. With the $(1\\!:\\!1)$ resonance of interested here, however, the
stabilized orbits oscillate around $\varphi=0$ and not around
$\varphi=\pm\pi$.

##  Appendix G Calculation of the Bohr-Sommerfeld critical photon numbers

In this appendix we describe the method used to compute the Bohr-Sommerfeld
critical photon numbers $\\{n_{\mathrm{crit},i_{t}}^{c}\\}$ introduced in Sec.
IV.2. The key is to determine whether the Bohr-Sommerfeld orbit associated
with a given qubit state $i_{t}$ exists. If the orbit is chaotic or overlaps
with a $(n\\!:\\!m)$ resonance, there exists no closed and fully connected
regular orbit with the required area, and the qubit state is then declared
unstable. The different photon numbers at which the Bohr-Sommerfeld orbit
associated to qubit state $i_{t}$ becomes nonexistent yield the Bohr-
Sommerfeld critical photon numbers $\\{n_{\mathrm{crit},i_{t}}^{c}\\}$.

For a given value of the rescaled drive frequency $\tilde{\omega}_{d}$ (i.e.,
a fixed value of the detuning $\Delta$) and of the rescaled drive amplitude
$\tilde{\varepsilon}_{t}$, we solve Hamilton’s equations to obtain the orbits
over $\mathcal{N}$ drive periods corresponding to all initial conditions along
the axis $(\tilde{\varphi}_{t}=0,\tilde{n}_{t}^{0})$, with
$\tilde{n}_{t}^{0}\in[-2.0,2.0]$ varied in increments of $1/300$. We then
attempt to compute the area of all orbits using Green’s theorem. More
precisely, we choose some phase-space field
$\mathbf{F}=\left\\{P(\tilde{\varphi}_{t},\tilde{n}_{t}),Q(\tilde{\varphi}_{t},\tilde{n}_{t})\right\\}$
such that $\nabla\times\mathbf{F}$ is a unit vector perpendicular to the
phase-space plane. With that choice, Green’s theorem yields

$\begin{split}\mathcal{A}&=\iint_{\mathcal{D}}(\nabla\times\mathbf{F})\cdot
d\mathbf{A}=\oint_{\mathcal{C}}\mathbf{F}\cdot\mathbf{dr}\\\
&\approx\sum_{\nu}\left[P(\tilde{\varphi}_{t}^{\nu},\tilde{n}_{t}^{\nu})\Delta\tilde{\varphi}_{t}^{\nu}+Q(\tilde{\varphi}_{t}^{\nu},\tilde{n}_{t}^{\nu})\Delta\tilde{n}_{t}^{\nu}\right],\end{split}$
(92)

where the sum is over all points $\nu$ in the orbit, and where
$\Delta\tilde{\varphi}_{t}^{\nu}=\tilde{\varphi}_{t}^{\nu}-\tilde{\varphi}_{t}^{\nu-1}$
and $\Delta\tilde{n}_{t}^{\nu}=\tilde{n}_{t}^{\nu}-\tilde{n}_{t}^{\nu-1}$ are
the coordinate differences between adjacent points of the orbit. A simple
choice for $\mathbf{F}$ is $P(\tilde{\varphi}_{t},\tilde{n}_{t})=0$ and
$Q(\tilde{\varphi}_{t},\tilde{n}_{t})=\tilde{\varphi}_{t}$, yielding

$\mathcal{A}\approx\left|\sum_{\nu}\tilde{\varphi}_{t}^{\nu}\Delta\tilde{n}_{t}^{\nu}\right|.$
(93)

The absolute value is present to obtain a positive area independently of the
integration direction along the orbit. In general, the points of an orbit in
the Poincaré section are not naturally ordered along the orbit because they
are obtained stroboscopically. We must therefore order the points of an orbit
to compute an area using Eq. 93. If most points of an orbit can be ordered on
a closed and connected near-continuous orbit, an area can be successfully
calculated and a closed orbit with that area can exist. In all other cases, a
closed connected orbit cannot exist.

In the attempt to order the points along an orbit, we start at the initial
point of the orbit. We then find the point in the Poincaré section that is
closest to the initial point in phase space. The nearest neighbor of this
second point is then selected, excluding the points that were already
selected. The procedure is repeated until it returns to the initial point or
until the nearest neighbor is farther than a threshold distance $d$ chosen as

$d=\begin{cases}\frac{2\pi
r_{c}m}{\mathcal{N}}&\text{if}\quad\frac{\mathcal{N}d_{\rm{min}}}{2\pi
m}<r_{c}<\frac{\mathcal{N}d_{\rm{max}}}{2\pi m},\\\
d_{\rm{min}}&\text{if}\quad r_{c}\leq\frac{\mathcal{N}d_{\rm{min}}}{2\pi
m},\\\ d_{\rm{max}}&\text{if}\quad
r_{c}\geq\frac{\mathcal{N}d_{\rm{max}}}{2\pi m}\text{ or no main
region.}\end{cases}$ (94)

Here, $m$ is a scalar and $r_{c}$ is the distance between the initial
condition $(0,\tilde{n}_{t}^{0})$ and the center of the main regular region
(determined as described below). The first condition in Eq. 94 sets the
distance threshold in proportion to the expected average distance between
orbit points. Indeed, the quantity $2\pi r_{c}/\mathcal{N}$ is the average arc
length between orbit points distributed uniformly on a circle of radius
$r_{c}$. Thus, $m$ is the number of average arc lengths tolerated between
nearest neighbors. The second and third conditions in Eq. 94 set a minimum and
maximum distance between nearest neighbors, respectively. The minimum distance
$d_{\rm{min}}$ ensures that the threshold $d$ remains finite very close to the
center of the main regular region, although such small orbits are not relevant
for the parameter regime explored in this work provided that the area
corresponding to a circle of radius $\mathcal{N}d_{\mathrm{min}}/2\pi m$ is
smaller than the smallest area a Bohr-Sommerfeld orbit can enclose, namely
$\mathcal{A}_{0}$. The maximal distance $d_{\rm{max}}$ ensures that the
threshold remains finite far from the main regular region or when the main
regular region stops existing. This mitigates the number of chaotic orbits for
which the algorithm occasionally finds a closed connected orbit (these appear
as “noise” in the white regions of Fig. 11). The algorithm records a
successful area calculation if it returns to the initial point and if the
number of ordered points is larger than $\mathcal{N}/2$ and otherwise records
a failure. The latter condition prevents identifying a single island of a
disconnected subharmonic resonant orbit as a closed connected orbit. For the
above procedure, we used $\mathcal{N}=1000$, $2\pi m/\mathcal{N}\approx 0.65$,
$d_{\rm{min}}=0.10$ and $d_{\rm{max}}=0.80$, which gave us fast and reliable
results.

Setting the threshold $d$ using Eq. 94 requires knowledge of the position of
the center of the main regular region at all drive amplitudes. At
$\tilde{\varepsilon}_{t}=0$, the center sits at $(0,0)$ in phase space. If the
drive amplitude is slightly increased, the center slightly moves away from its
former position such that an orbit with initial condition $(0,0)$ now encloses
the new center. The position of the new center is estimated as the center of
mass of that orbit. The method is iterated until the main regular region stops
existing. This is detected by the orbit wandering too far away from the
previous center position. To track the center, we increment
$\tilde{\varepsilon}_{t}$ by $0.01$, use orbits with $\mathcal{N}=1500$
periods, and set the wandering off threshold distance to $1.0$. In Fig. 11,
the center of the main regular region is shown as gray lines starting at
$(0,0)$ at zero drive amplitude in panels (a-b) while it appears as gray dots
in panels (c-d).

The above procedure allows us to compute the phase-space area of orbits as a
function of $\tilde{\omega}_{d}$ and $\tilde{\varepsilon}_{t}$. We first vary
$\tilde{\omega}_{d}$ [$\Delta$] between $0.5$
[$3.02\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$] and $1.5$
[$-3.49\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$] in increments of $0.01$, but
exclude very small detunings $\tilde{\omega}_{d}\in(0.90,1.03)$
[$\Delta\in(-0.47,0.41)\,\textrm{GHz}$] where the dispersive approximation is
not expected to hold. For each value of $\tilde{\omega}_{d}$, we then vary
$\tilde{\varepsilon}_{t}$ in increments of
$\delta\tilde{\varepsilon}_{t}=\delta\varepsilon_{t}/\omega_{p}=0.0015323$,
corresponding to
$\delta\varepsilon_{t}/2\pi=$10\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$$. The
results are shown in Fig. 11 for a fixed value of a (a) negative detuning and
(b) positive detuning. The colored regions indicate the initial conditions
where an area was successfully calculated. In particular, the Bohr-Sommerfeld
orbits are highlighted as pale lines. The dark colored regions represent the
extent of quantum fluctuations around a given Bohr-Sommerfeld orbit. More
precisely, they contain all orbits with an area within $\pi\hbar_{\rm eff}$ of
that orbit at fixed $\tilde{\varepsilon}_{t}$. The white regions indicate the
initial conditions where the algorithm failed to calculate an area. The Bohr-
Sommerfeld critical photon numbers $\\{n_{\mathrm{crit},i_{t}}^{c}\\}$ are the
points when the Bohr-Sommerfeld line for state $i_{t}$ collides with the white
regions as the drive amplitude is increased. These points are represented by
stars in Fig. 11. A lower bound to the $\\{n_{\mathrm{crit},i_{t}}^{c}\\}$
from quantum fluctuations is obtained from the points where the boundaries of
the dark colored regions collide with the white regions as the drive amplitude
is increased. These points are not shown in Fig. 11.

Finally, we quantify the strength of the resonances at the critical photon
numbers by the area they occupy in phase space, see Fig. 12 of the main text.
To obtain the area of a resonance, we start at a critical point in Fig. 11
(star). We then vary $\tilde{n}_{t}^{0}$ across the resonance at fixed
$\tilde{\varepsilon}_{t}$. The area of the resonance is then obtained as the
difference of the areas of the two orbits that jointly enclose the resonance.
We repeat this procedure for every relevant resonance.

## References

  * Blais  _et al._ (2021) Alexandre Blais, Arne L. Grimsmo, S. M. Girvin, and Andreas Wallraff, “Circuit quantum electrodynamics,” Rev. Mod. Phys. 93, 025005 (2021).
  * Koch  _et al._ (2007) Jens Koch, Terri M. Yu, Jay Gambetta, A. A. Houck, D. I. Schuster, J. Majer, Alexandre Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, “Charge-insensitive qubit design derived from the Cooper pair box,” Phys. Rev. A 76, 042319 (2007).
  * Wallraff  _et al._ (2005) A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, “Approaching unit visibility for control of a superconducting qubit with dispersive readout,” Phys. Rev. Lett. 95, 060501 (2005).
  * Blais  _et al._ (2004) Alexandre Blais, Ren-Shou Huang, Andreas Wallraff, S. M. Girvin, and R. J. Schoelkopf, “Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation,” Phys. Rev. A 69, 062320 (2004).
  * Walter  _et al._ (2017) T Walter, P Kurpiers, S Gasparinetti, P Magnard, A Potočnik, Y Salathé, M Pechal, M Mondal, M Oppliger, C Eichler, and A Wallraff, “Rapid high-fidelity single-shot dispersive readout of superconducting qubits,” Phys. Rev. Appl. 7, 054020 (2017).
  * Sank  _et al._ (2016) Daniel Sank, Zijun Chen, Mostafa Khezri, J. Kelly, R. Barends, B. Campbell, Y. Chen, B. Chiaro, A. Dunsworth, A. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, T. White, J. Wenner, Alexander N. Korotkov, and John M. Martinis, “Measurement-induced state transitions in a superconducting qubit: Beyond the rotating wave approximation,” Phys. Rev. Lett. 117, 190503 (2016).
  * Khezri  _et al._ (2023) Mostafa Khezri, Alex Opremcak, Zijun Chen, Kevin C. Miao, Matt McEwen, Andreas Bengtsson, Theodore White, Ofer Naaman, Daniel Sank, Alexander N. Korotkov, Yu Chen, and Vadim Smelyanskiy, “Measurement-induced state transitions in a superconducting qubit: Within the rotating-wave approximation,” Phys. Rev. Appl. 20, 054008 (2023).
  * Minev  _et al._ (2019) Z. K. Minev, S. O. Mundhada, S. Shankar, P. Reinhold, R. Gutiérrez-Jáuregui, R. J. Schoelkopf, M. Mirrahimi, H. J. Carmichael, and M. H. Devoret, “To catch and reverse a quantum jump mid-flight,” Nature 570, 200–204 (2019).
  * Swiadek  _et al._ (2023) François Swiadek, Ross Shillito, Paul Magnard, Ants Remm, Christoph Hellings, Nathan Lacroix, Quentin Ficheux, Dante Colao Zanuz, Graham J. Norris, Alexandre Blais, Sebastian Krinner, and Andreas Wallraff, “Enhancing dispersive readout of superconducting qubits through dynamic control of the dispersive shift: Experiment and theory,” (2023), arXiv:2307.07765 [quant-ph] .
  * Johnson  _et al._ (2012) J. E. Johnson, C. Macklin, D. H. Slichter, R. Vijay, E. B. Weingarten, John Clarke, and I. Siddiqi, “Heralded state preparation in a superconducting qubit,” Phys. Rev. Lett. 109, 050506 (2012).
  * Touzard  _et al._ (2019) S. Touzard, A. Kou, N. E. Frattini, V. V. Sivak, S. Puri, A. Grimm, L. Frunzio, S. Shankar, and M. H. Devoret, “Gated conditional displacement readout of superconducting qubits,” Phys. Rev. Lett. 122, 080502 (2019).
  * Sunada  _et al._ (2022) Y. Sunada, S. Kono, J. Ilves, S. Tamate, T. Sugiyama, Y. Tabuchi, and Y. Nakamura, “Fast readout and reset of a superconducting qubit coupled to a resonator with an intrinsic Purcell filter,” Phys. Rev. Appl. 17, 044016 (2022).
  * Lledó  _et al._ (2023) Cristóbal Lledó, Rémy Dassonneville, Adrien Moulinas, Joachim Cohen, Ross Shillito, Audrey Bienfait, Benjamin Huard, and Alexandre Blais, “Cloaking a qubit in a cavity,” Nat. Commun. 14, 6313 (2023).
  * Chen  _et al._ (2023) Liangyu Chen, Hang-Xi Li, Yong Lu, Christopher W. Warren, Christian J. Križan, Sandoko Kosen, Marcus Rommel, Shahnawaz Ahmed, Amr Osman, Janka Biznárová, Anita Fadavi Roudsari, Benjamin Lienhard, Marco Caputo, Kestutis Grigoras, Leif Grönberg, Joonas Govenius, Anton Frisk Kockum, Per Delsing, Jonas Bylander, and Giovanna Tancredi, “Transmon qubit readout fidelity at the threshold for quantum error correction without a quantum-limited amplifier,” npj Quantum Inf. 9, 26 (2023).
  * Kjaergaard  _et al._ (2020) Morten Kjaergaard, Mollie E. Schwartz, Jochen Braumüller, Philip Krantz, Joel I.-J. Wang, Simon Gustavsson, and William D. Oliver, “Superconducting qubits: Current state of play,” Annu. Rev. Condens. Matter Phys. 11, 369–395 (2020).
  * Boissonneault  _et al._ (2008) Maxime Boissonneault, J. M. Gambetta, and Alexandre Blais, “Nonlinear dispersive regime of cavity QED: The dressed dephasing model,” Phys. Rev. A 77, 060305 (2008).
  * Boissonneault  _et al._ (2009) Maxime Boissonneault, J. M. Gambetta, and Alexandre Blais, “Dispersive regime of circuit QED: Photon-dependent qubit dephasing and relaxation rates,” Phys. Rev. A 79, 013819 (2009).
  * Petrescu  _et al._ (2020) Alexandru Petrescu, Moein Malekakhlagh, and Hakan E. Türeci, “Lifetime renormalization of driven weakly anharmonic superconducting qubits. II. The readout problem,” Phys. Rev. B 101, 134510 (2020).
  * Hanai  _et al._ (2021) Ryo Hanai, Alexander McDonald, and Aashish Clerk, “Intrinsic mechanisms for drive-dependent Purcell decay in superconducting quantum circuits,” Phys. Rev. Res. 3, 043228 (2021).
  * Shillito  _et al._ (2022) Ross Shillito, Alexandru Petrescu, Joachim Cohen, Jackson Beall, Markus Hauru, Martin Ganahl, Adam G.M. Lewis, Guifre Vidal, and Alexandre Blais, “Dynamics of transmon ionization,” Phys. Rev. Appl. 18, 034031 (2022).
  * Cohen  _et al._ (2023) Joachim Cohen, Alexandru Petrescu, Ross Shillito, and Alexandre Blais, “Reminiscence of classical chaos in driven transmons,” PRX Quantum 4, 020312 (2023).
  * Xiao  _et al._ (2023) Xu Xiao, Jayameenakshi Venkatraman, Rodrigo G. Cortiñas, Shoumik Chowdhury, and Michel H. Devoret, “A diagrammatic method to compute the effective Hamiltonian of driven nonlinear oscillators,” (2023), arXiv:2304.13656 [quant-ph] .
  * Boissonneault  _et al._ (2010) Maxime Boissonneault, J. M. Gambetta, and Alexandre Blais, “Improved superconducting qubit readout by qubit-induced nonlinearities,” Phys. Rev. Lett. 105, 100504 (2010).
  * Shillito (2023) R. Shillito,  _Simulation and Optimization of Superconducting Qubit Control and Readout_ , Ph.D. thesis, Université de Sherbrooke (2023).
  * Carmichael (2007) Howard J Carmichael,  _Statistical Methods in Quantum Optics 2_ (Springer Berlin Heidelberg, 2007).
  * Note (1) Emitting 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.
  * Muñoz Arias  _et al._ (2023) Manuel H. Muñoz Arias, Cristóbal Lledó, and Alexandre Blais, “Qubit readout enabled by qubit cloaking,” Phys. Rev. Appl. 20, 054013 (2023).
  * Breuer and Holthaus (1989) H.P. Breuer and M. Holthaus, “Quantum phases and Landau-Zener transitions in oscillating fields,” Phys. Lett. A 140, 507–512 (1989).
  * Grifoni and Hänggi (1998) Milena Grifoni and Peter Hänggi, “Driven quantum tunneling,” Phys. Rep. 304, 229 – 354 (1998).
  * Hone  _et al._ (1997) Daniel W. Hone, Roland Ketzmerick, and Walter Kohn, “Time-dependent Floquet theory and absence of an adiabatic limit,” Phys. Rev. A 56, 4045–4054 (1997).
  * Drese and Holthaus (1999) K. Drese and M. Holthaus, “Floquet theory for short laser pulses,” Eur. Phys. J. D 5, 119–134 (1999).
  * Oliver  _et al._ (2005) William D. Oliver, Yang Yu, Janice C. Lee, Karl K. Berggren, Leonid S. Levitov, and Terry P. Orlando, “Mach-Zehnder interferometry in a strongly driven superconducting qubit,” Science 310, 1653–1657 (2005).
  * Shevchenko  _et al._ (2010) S.N. Shevchenko, S. Ashhab, and Franco Nori, “Landau–Zener–Stückelberg interferometry,” Phys. Rep. 492, 1–30 (2010).
  * Ikeda  _et al._ (2022) Tatsuhiko N. Ikeda, Satoshi Tanaka, and Yosuke Kayanuma, “Floquet-Landau-Zener interferometry: Usefulness of the Floquet theory in pulse-laser-driven systems,” Phys. Rev. Res. 4, 033075 (2022).
  * Kramer and MacKinnon (1993) B Kramer and A MacKinnon, “Localization: theory and experiment,” Rep. Prog. Phys. 56, 1469 (1993).
  * Zaslavsky (2005) G.M. Zaslavsky,  _Hamiltonian Chaos and Fractional Dynamics_ (OUP Oxford, 2005).
  * Wimberger (2014) Sandro Wimberger,  _Nonlinear dynamics and quantum chaos: An introduction_ (Springer, Berlin, 2014) p. 206.
  * Delone  _et al._ (1983) N B Delone, B P Kraĭnov, and D L Shepelyanskiĭ, “Highly-excited atoms in the electromagnetic field,” Sov. phys., Usp. 26, 551 (1983).
  * Casati (1987) G Casati, “Relevance of classical chaos in quantum mechanics: The hydrogen atom in a monochromatic field,” Phys. Rep. 154, 77–123 (1987).
  * Wang and Chu (1989) Kwanghsi Wang and Shih-I Chu, “Dynamics of multiphoton excitation and quantum diffusion in Rydberg atoms,” Phys. Rev. A 39, 1800–1808 (1989).
  * Messiah (2020) Albert Messiah,  _Quantum mechanics_ (Dover Publications, Inc., Garden City, New York, 2020).
  * (42) P. Kurilovich, T. Connolly, C. Boettcher, S. Diamond, H. Nho, A. Ding, D. Weiss, V. Kurilovich, V. Fatemi, L. Glazman, and M. Devoret, In preparation.
  * Manucharyan  _et al._ (2009) Vladimir E. Manucharyan, Jens Koch, Leonid I. Glazman, and Michel H. Devoret, “Fluxonium: Single Cooper-pair circuit free of charge offsets,” Science 326, 113–116 (2009).
  * Hyyppä  _et al._ (2022) Eric Hyyppä, Suman Kundu, Chun Fai Chan, András Gunyhó, Juho Hotari, David Janzso, Kristinn Juliusson, Olavi Kiuru, Janne Kotilahti, Alessandro Landra, _et al._ , “Unimon qubit,” Nat. Commun. 13, 6895 (2022).
  * Didier  _et al._ (2015) Nicolas Didier, Jérôme Bourassa, and Alexandre Blais, “Fast quantum nondemolition readout by parametric modulation of longitudinal qubit-oscillator interaction,” Phys. Rev. Lett. 115, 203601 (2015).
  * (46) Alex Chapple, Alexander McDonald, Manuel Muños-Arias, Christian de Correc, Michelle Lachapelle, and Alexandre Blais, In preparation.
  * Caldwell  _et al._ (2018) S. A. Caldwell, N. Didier, C. A. Ryan, E. A. Sete, A. Hudson, P. Karalekas, R. Manenti, M. P. da Silva, R. Sinclair, E. Acala, N. Alidoust, J. Angeles, A. Bestwick, M. Block, B. Bloom, A. Bradley, C. Bui, L. Capelluto, R. Chilcott, J. Cordova, G. Crossman, M. Curtis, S. Deshpande, T. El Bouayadi, D. Girshovich, S. Hong, K. Kuang, M. Lenihan, T. Manning, A. Marchenkov, J. Marshall, R. Maydra, Y. Mohan, W. O’Brien, C. Osborn, J. Otterbach, A. Papageorge, J.-P. Paquette, M. Pelstring, A. Polloreno, G. Prawiroatmodjo, V. Rawat, M. Reagor, R. Renzas, N. Rubin, D. Russell, M. Rust, D. Scarabelli, M. Scheer, M. Selvanayagam, R. Smith, A. Staley, M. Suska, N. Tezak, D. C. Thompson, T.-W. To, M. Vahidpour, N. Vodrahalli, T. Whyland, K. Yadav, W. Zeng, and C. Rigetti, “Parametrically activated entangling gates using transmon qubits,” Phys. Rev. Appl. 10, 034050 (2018).
  * Niskanen  _et al._ (2007) A O Niskanen, K Harrabi, F Yoshihara, Y Nakamura, S Lloyd, and J S Tsai, “Quantum coherent tunable coupling of superconducting qubits,” Science 316, 723–726 (2007).
  * Puri  _et al._ (2017) Shruti Puri, Samuel Boutin, and Alexandre Blais, “Engineering the quantum states of light in a Kerr-nonlinear resonator by two-photon driving,” npj Quantum Inf. 3, 18 (2017).
  * Frattini  _et al._ (2022) Nicholas E Frattini, Rodrigo G Cortiñas, Jayameenakshi Venkatraman, Xu Xiao, Qile Su, Chan U Lei, Benjamin J Chapman, Vidul R Joshi, S M Girvin, Robert J Schoelkopf, Shruti Puri, and Michel H Devoret, “The squeezed kerr oscillator: spectral kissing and phase-flip robustness,” (2022), arXiv:2209.03934 [quant-ph] .
  * Chávez-Carlos  _et al._ (2023) Jorge Chávez-Carlos, Rodrigo G Cortiñas, Miguel A Prado Reynoso, Ignacio García-Mata, Victor S Batista, Francisco Pérez-Bernal, Diego A Wisniacki, and Lea F Santos, “Driving superconducting qubits into chaos,” (2023), arXiv:2310.17698 [quant-ph] .
  * D’Anjou and Burkard (2019) Benjamin D’Anjou and Guido Burkard, “Optimal dispersive readout of a spin qubit with a microwave resonator,” Phys. Rev. B 100, 245427 (2019).
  * Schrieffer and Wolff (1966) J. R. Schrieffer and P. A. Wolff, “Relation between the anderson and kondo hamiltonians,” Phys. Rev. 149, 491–492 (1966).
  * Shirley (1965) Jon H. Shirley, “Solution of the Schrödinger equation with a Hamiltonian periodic in time,” Phys. Rev. 138, B979–B987 (1965).