# Reminiscence of classical chaos in driven transmons

Joachim Cohen, Alexandru Petrescu, Ross Shillito, Alexandre Blais

► To cite this version:

Joachim Cohen, Alexandru Petrescu, Ross Shillito, Alexandre Blais. Reminiscence of classical chaos in driven transmons. PRX Quantum, 2023, 4 (2), pp.020312. 10.1103/PRXQuantum.4.020312. hal-03758792

HAL Id: hal-03758792

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

Submitted on 12 Jul 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Reminiscence of Classical Chaos in Driven Transmons

Joachim Cohen,¹,*, Alexandru Petrescu,¹,²,† Ross Shillito,¹ and Alexandre Blais¹,³

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

²LPENS, Département de physique, Ecole normale supérieure, Centre Automatique et Systèmes (CAS), MINES ParisTech, Université PSL, Sorbonne Université, CNRS, Inria, Paris 75005, France

³Canadian Institute for Advanced Research, Toronto, Ontario M5G1M1, Canada

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

(Received 9 November 2022; accepted 21 March 2023; published 21 April 2023)

Transmon qubits are ubiquitously used in superconducting quantum information processor architectures. Strong drives are required to realize fast, high-fidelity, gates and measurements, including parametrically activated processes. Here, we show that even off-resonant drives, in regimes routinely used in experiments, can cause strong modifications to the structure of the transmon spectrum rendering a large part of it chaotic. Accounting for the full nonlinear dynamics of the transmon in a Floquet-Markov formalism, we find that these chaotic states, often neglected through the hypothesis that the anharmonicity is weak, strongly impact the lifetime of the transmon's computational states. In particular, we observe that chaos-assisted quantum phase slips greatly enhance band dispersions. In the presence of a measurement resonator, we find that approaching chaotic behavior correlates with strong transmon-resonator hybridization, and an average resonator response centered on the bare resonator frequency. These results lead to a photon-number threshold characterizing the appearance of chaos-induced quantum demolition effects during strong-drive operations, such as dispersive qubit readout. The phenomena described here are expected to be present in all circuits based on low-impedance Josephson junctions.

DOI: 10.1103/PRXQuantum.4.020312

## I. INTRODUCTION

Low-impedance Josephson-junction circuits, where the Josephson energy dominates over the charging energy, are fundamental building blocks of superconducting quantum processors. Although the most widely used superconducting qubit is the transmon [1], capacitively shunted Josephson junctions appear in other species of qubits, such as the heavy fluxonium [2,3], the $0 - \pi$ qubit [4] and the capacitively shunted flux qubit [5,6]. Josephson junctions can also be used as simple nonlinear elements for parametrically activated multiwave mixing [7–10], or as linear inductive elements in a Josephson junction array to realize superinductances [11,12].

To meet the requirements of quantum information processing with fast and high-fidelity operations, strong driving fields that are off-resonant from the qubit are often

used in parametrically activated coupling [7,13], multi-qubit gates [14,15], and dispersive readout [16]. However, off-resonant drives of even moderate amplitude are often observed to cause spurious qubit transitions [17–19]. This is the case of the dispersive readout whose quantum non-demolition (QND) character is observed only at very small drive amplitudes, corresponding to a few photons ($\bar{n} \sim 2$) populating the measurement resonator [17]. Models based on perturbative expansion in the qubit-resonator coupling [20] or in qubit anharmonicity and drive amplitude [21,22] have been explored to understand the origin of these unwanted transitions. At large photon numbers ($\bar{n} \gtrsim 100$), qubit-resonator resonances [23] and structural instabilities [13,24–26] have been shown to result in spurious transitions. Recently, a numerical study of the full time dynamics of the transmon has shown that qubit-resonator resonances can lead to leakage of the transmon population to states lying above the Josephson-junction cosine potential, something that has been referred to as ionization [27]. In that study, ionization was shown to coincide with the loss of QNDness experimentally observed at low photon numbers in Ref. [17].

Here, we show that the often-neglected highly excited states of the transmon can play an important role for current experimental parameters and drive amplitudes, even in the absence of ionization [13,23,28]. We find that a

<footnote>*cohenj38@gmail.com</footnote>

<footnote>†alexandru.petrescu@minesparis.psl.eu</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
subset of the highly excited states in the spectrum of a driven transmon, the so-called chaotic layer, significantly and unexpectedly dresses the charge dispersion of the low-energy transmon spectrum, with detrimental effects on the qubit dephasing time. This phenomenon can be understood as chaos-assisted [29] quantum phase slips [30]. The increased dependence on the offset charge even for the low-energy states suggests that in the presence of strong drives, models, which rely on a perturbed harmonic oscillator such as the Kerr nonlinear oscillator [31], cannot give an accurate description of the system, in particular, because the selection rules derived from these models are no longer applicable [23]. In addition, we show that the chaotic layer makes steady-state populations deviate significantly from the Boltzmann distribution [26,32,33].

We draw upon the Floquet theory [32–34] of nonlinear oscillators [35] to distinguish between the chaotic and regular states [36,37] of a driven transmon. We introduce a rescaling of the transmon Hamiltonian from which an effective Planck constant emerges, $\hbar_{\text{eff}} = \sqrt{8E_C/E_J}$ ($\hbar = 1$). In the transmon regime, where $\hbar_{\text{eff}}$ is small, the chaotic dynamics of the classical driven transmon becomes more resolved in the quantum spectrum as the number of chaotic states increases. In particular, we show that the spectrum of a single driven transmon is correlated [38]. In transmon systems, this type of analysis has been performed in the many-body regime [39]. Our study reveals that, within the range of current experimental parameters, the size of the classical chaotic domain is strongly sensitive to the drive frequency, something which can lead to instabilities in the quantum dynamics even at low drive power. Our results suggest ways to avoid experimental realizations of transmons from being plagued by these instabilities.

Moreover, by simulating the full transmon plus resonator circuit QED model, we show that chaos develops along two directions, that of increasing resonator Fock state number and that of increasing drive power. In addition to validating the study of the single driven transmon, we show that the transmon and the resonator strongly hybridize in the chaotic phase. As a result, in the context of the dispersive readout where chaotic effects are present, spurious qubit transitions become possible. In particular, we introduce a critical photon number around which chaos-induced non-QND effects are expected. We also predict that spurious effects below this threshold should be exponentially reduced with $\sqrt{E_J/8E_C}$.

The remainder of this paper is organized as follows. In Sec. II, we introduce a rescaled version of the Hamiltonian along with the key parameters of the dynamics, discuss chaos in the classical driven transmon, and briefly discuss the impact of classical chaos on the quantum system. Section III tackles the spectral properties of the system and the dependence of the instability on drive frequency. In Sec. IV, we consider the coupling of the transmon to a bath, and the impact of the chaotic layer on the coherence

properties of the low-energy sector is analyzed. Section V focuses on the interplay of chaos and the validity of the dispersive approximation in the full-circuit QED setup. Section VI discusses spurious non-QND effects originating from the interaction with chaotic states.

## II. PERIODICALLY DRIVEN TRANSMON

The Hamiltonian of a transmon in a typical circuit QED setup takes the form [1,40]

$$ \mathbf{H}(t) = 4E_C(\mathbf{n} - n_g)^2 - E_J \cos(\phi) + \mathbf{n}\mathbf{F}, \quad (1) $$

where $\mathbf{n}$ and $\phi$ are, respectively, the charge and phase operators, $E_C$ and $E_J$ are the charging and Josephson energies, and $n_g$ an offset charge. The phase $\phi$ is compact and takes its values in the range $(-\pi, \pi]$. The operator $\mathbf{F}$ represents a classical driving field on the resonator or the coupling to a measurement resonator in a circuit QED setup [40]. As such, $\mathbf{F} = F_c(t) + \mathbf{F}_q$ can be expressed as the sum of a classical part and of a quantum part. The quantum part represents the displaced quadrature of the readout resonator, i.e., $\mathbf{F}_q = ig(\mathbf{a}^\dagger - \mathbf{a})$, where $g$ is the light-matter coupling. On the other hand, the classical part is assumed to take the form $F_c(t) = \varepsilon_d \cos(\omega_d t)$ and represents either a direct capacitive drive on the transmon or the classical amplitude of the resonator field. In this work, we are concerned with the limit of strong drives, where the classical part dominates over quantum fluctuations, i.e., $\varepsilon_d > g$. Expressing the drive amplitude $\varepsilon_d = 2g\sqrt{n}$ in terms of an equivalent number of resonator photons $\tilde{n}$, this strong drive limit corresponds to $\tilde{n} > 1$. In typical circuit QED (cQED) setup, one has $g/2\pi \sim 250$ MHz and $\tilde{n} \ge 2$, resulting in an effective drive $\varepsilon_d/2\pi \ge 700$ MHz.

In what follows, we neglect quantum fluctuations $\mathbf{F}_q$ to study the periodically driven transmon Hamiltonian

$$ \mathbf{H}(t) = 4E_C(\mathbf{n} - n_g)^2 - E_J \cos(\phi) + \varepsilon_d \cos(\omega_d t)\mathbf{n}. \quad (2) $$

We return to a full circuit QED model accounting for the presence of the resonator in Sec. V.

### A. Classical model

We first consider the classical limit of the driven transmon Hamiltonian where we replace the conjugate operators $\{\phi, \mathbf{n}\}$ by the phase-space coordinates $\{\phi, n\}$. In doing so, it is useful to rescale energy and time using the relations $\tilde{H} = H/E_J$ and $\tilde{t} = \omega_p t$, where $\omega_p = \sqrt{8E_J E_C}/\hbar$ is the plasma frequency of the transmon. Under this transformation, which preserves Hamilton’s equations, the classical Hamiltonian takes the form

$$ \tilde{H}(\tilde{t}) = \frac{(\tilde{n} - \tilde{n}_g)^2}{2} - \cos \tilde{\phi} + \tilde{\varepsilon}_d \cos(\tilde{\omega}_d \tilde{t}) \tilde{n}. \quad (3) $$

This corresponds to the Hamiltonian of a driven charged classical pendulum with dimensionless momentum $\tilde{n} = zn$
and position $\tilde{\phi} = \phi$ [1]. Here, $z = \sqrt{8E_C/E_J}$ is the characteristic impedance of the transmon. In this rescaled form, the three relevant parameters of the classical driven transmon are the rescaled drive amplitude $\tilde{\varepsilon}_d = \varepsilon_d/\omega_p$, the rescaled drive frequency $\tilde{\omega}_d = \omega_d/\omega_p$ and the rescaled offset charge $\tilde{n}_g = zn_g$.

In the absence of a drive, two different types of motion of the system can be distinguished. For $\tilde{H} < 1$ (i.e., $H < E_J$), the system undergoes small and bounded phase oscillations. On the other hand, for $\tilde{H} > 1$, the system experiences unbounded full $\pm 2\pi$ rotations of the phase. While manipulations of the transmon qubit are designed such as to only lead to small phase oscillations, the transmon can be promoted to states above $E_J$ by strong drives [13,23,26,27]. In the quantum case, the resulting full rotations correspond to quantum phase slips [1,30]. At the boundary of these two types of motion, defined by the trajectory of energy $\tilde{H} = 1$ also known as the separatrix, small perturbations can have a large impact, for example, causing bounded oscillations to turn into unbounded rotations. This structural instability results in an irregular, chaotic motion of the pendulum at finite drive amplitude in the vicinity of the separatrix [35].

A useful representation of the system dynamics, both regular and chaotic, is provided by the stroboscopic Poincaré sections obtained by plotting the value of the phase-space coordinates $\{\tilde{\phi}(\tilde{t}), \tilde{n}(\tilde{t})\}$ at every period $\tilde{T} = \omega_p/\omega_d$ of the drive for some initial condition [41]. Figure 1(a) first shows this in the absence of a drive. There, the two expected types of motions are clearly visible: the bounded oscillations leading to the closed orbits and the unbounded rotations to the nearly horizontal patterns. In the presence of the drive, see Fig. 1(b), the Poincaré sections break up into regular and chaotic regions. The regular regions consist of weakly perturbed Kolmogorov-Arnold-Moser tori, reminiscent of the motion of the unperturbed system, while the chaotic region develops around the separatrix [35]. The small tori located within the regular unbounded trajectories in Fig. 1(b) are due to two resonances where the drive frequency $\pm \tilde{\omega}_d$ matches the energies of the trajectories, which pass in the vicinity of $(\tilde{\phi}, \tilde{n}) = (0, \pm 2.5)$.

## B. Quantum model

To compare the quantum dynamics to the classical one, we quantize the rescaled Hamiltonian of Eq. (3). Importantly, because the rescaling does not preserve the phase-space volume, it leads to a renormalization of the Planck constant upon quantization with $\hbar_{\text{eff}} = \hbar\omega_p/E_J = z$. As a result, the commutation relation of the rescaled operators is $[\tilde{\phi}, \tilde{n}] = i\hbar_{\text{eff}}$. Consequently, in addition to the three parameters that determine the classical dynamics enumerated above, $\tilde{\varepsilon}_d$, $\tilde{\omega}_d$, and $\tilde{n}_g$, the driven quantum dynamics

<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Poincaré section of the undriven system and (b) of the system driven at $\tilde{\varepsilon}_d = 0.5$, $\tilde{\omega}_d = 1.34$ and with initial time $\tilde{T}/8$. In the driven case, a chaotic layer develops around the separatrix. Additional resonances within the unbounded states appear taking the form of tori located at the border of the separatrix. (c),(d) Mean energy per cycle as a function of drive amplitude $\tilde{\varepsilon}_d$ for $\tilde{\omega}_d = 1.34$, $n_g = 0.25$, and for two different values of $\hbar_{\text{eff}}$: (c) $\hbar_{\text{eff}}^{-1} = 7.91$ ($E_J/E_C = 500$) and (d) $\hbar_{\text{eff}}^{-1} = 3$ ($E_J/E_C = 72$). In analogy to the classical case, strong hybridization of the states develops about the separatrix located at $\langle\langle H \rangle\rangle/E_J = 2$. The insets marked by a red triangle show the Husimi functions of the Floquet modes of the first excited state at $\tilde{\varepsilon}_d = 0.5$. Its wave function remains globally regular and is reminiscent of the first Fock state. The insets marked by a red dot correspond to a state located close to the separatrix. The corresponding Husimi function is irregular and delocalized over the corresponding classical chaotic domain in (b). The black squares represent the phase-space area $\hbar_{\text{eff}}$ occupied by one state.</figcaption></figure>

of the transmon is characterized by a fourth parameter, $\hbar_{\text{eff}}$, characterizing quantum fluctuations.

The solutions to the time-dependent Schrödinger equation associated to the Hamiltonian in Eq. (2) are the time-dependent Floquet states $|\psi_k(t)\rangle = \exp(-i\varepsilon_k t)|\phi_k(t)\rangle$, characterized by the rescaled quasienergies $\tilde{\varepsilon}_k$ and the time-periodic Floquet modes $|\phi_k(t)\rangle$ [34]. The quasienergies of the Floquet modes are defined up to integer multiples of the drive frequencies, and hence are
not indicative of the amount of energy stored in the system. However, the mean energy per cycle for mode $|\phi_k(t)\rangle$ can be defined as [42]

$$ \langle\langle H \rangle\rangle = \frac{1}{T} \int_{0}^{T} dt \langle\phi_k(t)| \mathbf{H}(t) |\phi_k(t)\rangle, \quad (4) $$

where T is the period of the drive. This mean energy per cycle is plotted in Fig. 1(c) as a function of the drive amplitude $\tilde{\epsilon}_d$ and different Floquet states k for a very weakly non-linear transmon qubit with $\hbar_{\text{eff}}^{-1} = 7.91$ ($E_J/E_C = 500$). As the drive amplitude increases, the perturbation hybridized the states around the energy $\langle\langle H \rangle\rangle/E_J = 1$, corresponding to the energy of the separatrix in the classical system, and does so in an increasingly large bandwidth around that energy. Remarkably the main features subsist in the more experimentally relevant case of $\hbar_{\text{eff}}^{-1} = 3$ ($E_J/E_C = 72$), although the hybridization is visually less pronounced because the level separations are larger, see Fig. 1(d).

The diffusion of classical trajectories through the chaotic domain translates to delocalized Husimi functions of the Floquet modes for the driven quantum system [33,43]. This can be intuitively understood by the fact that the Floquet modes are eigenstates of the propagator over one period of the drive, which is the quantum analog of the stroboscopic Poincaré map defined above [44]. In the insets of Figs. 1(c) and 1(d), we plot the Husimi functions at time $T/8$ of the Floquet modes indicated by the red markers on the spectra at $\tilde{\epsilon}_d = 0.5$. Because the phase is defined on $(-\pi, \pi]$, we use the definition of Ref. [45] for a coherent state on a circle. The insets indicated by the red triangle correspond to the Husimi functions of the Floquet modes of the first excited state. Because it is outside of the region where the states are strongly mixed (i.e., outside of the chaotic layer in the classical system), its wave function remains globally regular and is reminiscent of the first Fock state. On the other hand, the insets marked by a red dot correspond to a state located close to the separatrix. The corresponding Husimi functions are irregular and delocalized over the region corresponding to the chaotic layer in the classical case, see Fig. 1(b). The black squares in the upper right corner of the insets indicate the phase-space area $\hbar_{\text{eff}}$ occupied by one state. A smaller $\hbar_{\text{eff}}$ results in a smaller amplitude of quantum fluctuations, and therefore in more resolved features in the Husimi functions.

### III. QUANTUM SIGNATURES OF CHAOS IN THE DRIVEN SPECTRUM

#### A. Level-spacing statistics

The chaotic nature of the states near the separatrix can be confirmed through the correlated nature of the Floquet spectrum, which manifests itself in the distribution of

<figure><img src="image_4.png" /><figcaption>FIG. 2. Level-spacing statistics for $\hbar_{\text{eff}}^{-1} = 3$ and $\tilde{\omega}_d = 1.34$ for the undriven (green) and driven (red) spectrum. The set of states comprises those of mean energy satisfying $0.6 < \langle\langle H \rangle\rangle/E_J < 1.5$. The statistics is generated using the Floquet spectra corresponding to 200 values of $n_g$ uniformly distributed over the interval $[0, 0.5]$. The driven spectrum distribution follows the Wigner-Dyson distribution (dashed line), while the distribution of energies for the undriven system follows the Poisson distribution (dotted line). The inset shows the integrated distribution $I(\Delta) = \int_0^\Delta dsP(s)$.</figcaption></figure>

level spacings $P(\Delta)$. Indeed, the strong level hybridization observed near the separatrix results in level repulsion and, in turn, to a Floquet spectrum with strong correlations in the distribution of the spacing between levels. In that situation, $P(\Delta)$ is expected to follow the Wigner-Dyson distribution (dashed line in Fig. 2) [38]. In contrast, in the regular regime levels are uncorrelated and their distribution is Poissonian, as is characteristic of random uncorrelated events (dotted line in Fig. 2) [38].

To display this distribution, we consider a set of N states with sorted quasienergies $\epsilon_1 \le \dots \le \epsilon_N$ lying within the first Brillouin zone, i.e., $|\epsilon_n| \le \omega_d/2$ for all n. The spacing between adjacent quasienergy levels is defined by $\Delta_n = (\epsilon_{n+1} - \epsilon_n)/\bar{\Delta}$ for $n = 1, \dots, N-1$, where $\bar{\Delta} = \omega_d/N$ is the mean level spacing. We also define $\Delta_N = (\epsilon_0 - \epsilon_N + \omega_d)/\bar{\Delta}$ at the boundary of the Brillouin zone. In the limit of large N, the distribution $P(\Delta)$ is expected to follow the Wigner-Dyson statistics for correlated spectrum [46], whereas it is expected to follow a Poisson distribution for an uncorrelated spectrum.

In analogy with the study of classical chaos, we define the chaotic domain using the mean energy per cycle $\langle\langle H \rangle\rangle$. As can first be seen in Fig. 1(c) for $\hbar_{\text{eff}}^{-1} = 7.91$, the mean energies of the chaotic states are concentrated around $E_J$, and are separated from the regular states by a gap. The relevant energy bandwidth of the chaotic zone depends on the drive amplitude and, for $\tilde{\epsilon}_d = 0.5$, we take the N states whose energies satisfy $0.6 < \langle\langle H \rangle\rangle/E_J < 1.5$. However, for the value of $\hbar_{\text{eff}}^{-1} = 3$ corresponding to a typical transmon qubit, the spectrum displays only a handful of chaotic states in this energy bandwidth. With about $N \sim 7$ levels, see Fig. 1(d), this is far from enough to generate
meaningful level-spacing statistics. To circumvent this problem, we use the offset charge $n_g$ uniformly distributed in the interval $[0, 0.5]$ to generate statistics. This is possible because, even in the transmon regime, the chaotic states are not confined in the cosine potential well making them very sensitive to the offset charge. We return in Sec. IV to some of the observable implications of this $n_g$ dependence of the chaotic states.

We plot in Fig. 2 the cumulative histograms of $P(\Delta)$ obtained from adding the distributions corresponding to 200 values of $n_g$ uniformly distributed in the interval $[0, 0.5]$ for $\tilde{\epsilon}_d = 0$ (green bars) and $\tilde{\epsilon}_d = 0.5$ (red bars). In the absence of a drive, the system is regular and $P(\Delta)$ follows the Poisson distribution expected for an uncorrelated spectrum. This also suggests that the undriven spectrum is well randomized by the variation of $n_g$. In contrast, at $\tilde{\epsilon}_d = 0.5$, $P(\Delta)$ approaches the Wigner-Dyson distribution. The number of states $N$ in the selected energy bandwidth has increased in the chaotic case. This is due to the hybridization of states with mean energy close to the separatrix region $\langle\langle H \rangle\rangle \sim E_J$, as mentioned above. More generally, the extent of the chaotic layer is seen to increase with the drive amplitude, something which is particularly clear in Fig. 1(c). The inset of Fig. 2 shows the integrated distribution $I(\Delta) = \int_0^\Delta ds P(s)$ for which the statistical variations are reduced because of the integration, thereby allowing for a clearer distinction between uncorrelated and correlated spectra [38].

## B. Drive-frequency dependence

In addition to spreading over an increasingly large energy bandwidth with increasing drive amplitude, the size of the chaotic region also strongly depends on the drive frequency. This is illustrated in Figs. 3(a)–3(d), which shows the Poincaré sections for $\tilde{\omega}_d = 0.75, 1.25, 2.25$, and $4.25$. The mean energy spectra of the corresponding quantum systems are shown in Figs. 3(e)–3(h) as a function of $\tilde{\epsilon}_d/\tilde{\omega}_d$ for $\hbar_{eff}^{-1} = 2.45$ and $n_g = 0.25$. To compare the effect of different drive frequencies, the ac-Stark shift computed from the quasienergies of the ground and first excited states is fixed to 100 MHz in (a)–(d) (see Sec. A). The corresponding drive amplitudes are represented by a vertical black dashed line in Figs. 3(e)–3(h). In (a)–(d), the width of the chaotic layer is observed to be maximal for intermediate drive frequencies $\tilde{\omega}_d \sim 1-2$. This observation is in qualitative agreement with the approximate expression for the width of the chaotic layer around the separatrix

$$W_c/E_J \approx \tilde{\epsilon}_d \tilde{\omega}_d \text{sech}\left(\frac{\pi \tilde{\omega}_d}{2}\right), \quad (5)$$

a result, which is valid for $\tilde{\omega}_d > 1$ [41,47]. Following this expression, the chaotic layer is expected to have a maximal width for $1 < \tilde{\omega}_d < 2$, and to exponentially decrease in width with increasing $\tilde{\omega}_d$ for $\tilde{\omega}_d > 2$.

The approximate expression for the width, however, does not account for resonances that result in the additional tori observed in Figs. 3(a)–3(d). In particular, in the range $0.6 \lesssim \tilde{\omega}_d \lesssim 3$, which includes the regime of operation of current experiments, we find that resonances play an essential role in drive-induced instabilities. The origin of these resonances can be qualitatively understood by expressing the classical Hamiltonian of Eq. (3) in the equivalent form

$$\tilde{H}(\tilde{t}) = \frac{(\tilde{n} - \tilde{n}_g)^2}{2} - \cos\left[\tilde{\phi} + \frac{\tilde{\epsilon}_d}{\tilde{\omega}_d} \sin(\tilde{\omega}_d \tilde{t})\right]. \quad (6)$$

Using the Jacobi-Anger expansion to first order in $\tilde{\epsilon}_d/\tilde{\omega}_d$, this can be approximated as

$$\tilde{H}(\tilde{t}) \approx \frac{(\tilde{n} - \tilde{n}_g)^2}{2} - J_0\left(\frac{\tilde{\epsilon}_d}{\tilde{\omega}_d}\right) \cos\tilde{\phi} + J_1\left(\frac{\tilde{\epsilon}_d}{\tilde{\omega}_d}\right) \left[\cos(\tilde{\phi} - \tilde{\omega}_d\tilde{t}) - \cos(\tilde{\phi} + \tilde{\omega}_d\tilde{t})\right], \quad (7)$$

where $J_k(z)$ is the $k$th Bessel functions of the first kind. The first line of Eq. (7), $\tilde{H}_{\tilde{\epsilon}_d} = (\tilde{n} - \tilde{n}_g)^2/2 - J_0(\tilde{\epsilon}_d/\tilde{\omega}_d) \cos\tilde{\phi}$, describes an undriven pendulum with a potential energy reduced by the factor $J_0(\tilde{\epsilon}_d/\tilde{\omega}_d)$. The second line is of smaller amplitudes and describes the first harmonic of the time-dependent perturbation. Higher-order harmonics of the drive are neglected because their amplitude is suppressed by the factor $J_k(\tilde{\epsilon}_d/\tilde{\omega}_d)$ and because the corresponding resonances occur at higher frequency, $\tilde{\omega} > 2\tilde{\omega}_d$. The effects of the perturbation can be understood by inserting in the second line of Eq. (7) the solution $(\tilde{\phi}(t), \tilde{n}(t))$ of the system under the time-independent Hamiltonian $\tilde{H}_{\tilde{\epsilon}_d}$ with initial conditions $(\tilde{\phi}_0, \tilde{n}_0)$. Depending on the drive frequency, resonances can occur either within the bounded states or the unbounded states.

The trajectories representing the bounded states of the pendulum can be generated in the Poincaré sections with the initial conditions $(\tilde{\phi}_0 < \pi, \tilde{n}_0 = 0)$. For small oscillations, the pendulum behaves as a slightly anharmonic oscillator, and its oscillation frequency decreases as the oscillation amplitude $\tilde{\phi}_0$ increases. For $|\tilde{\phi}_0| < 0.8\pi$, the nonlinearity does not play an important role and the pendulum frequency varies smoothly. This corresponds to a frequency range $0.65 \lesssim \tilde{\omega} < 1$ for the pendulum oscillation. In this case, under $\tilde{H}_{\tilde{\epsilon}_d}$ the trajectory that passes through $(\tilde{\phi}_0, 0)$ takes the standard form $\tilde{\phi}(t) \approx \tilde{\phi}_0 \sin(\tilde{\omega}t)$ where we have neglected the higher harmonics of the motion. Inserting this expression for $\tilde{\phi}(t)$ in the second line of Eq. (7) leads to the slowly rotating terms $J_1(\tilde{\epsilon}_d/\tilde{\omega}_d)J_{2k+1}(\phi_0) \cos[(\tilde{\omega}_d - (2k+1)\tilde{\omega})\tilde{t}]$.

The case $k=0$ corresponds to a 1:1 resonance, i.e., $\tilde{\omega}_d = \tilde{\omega} \in [0.65, 1]$, which strongly impacts the low-energy states of the system. At $\tilde{\omega}_d = 0.75$ this resonance appears as a second set of tori located close to the central
<figure><img src="image_5.png" /><figcaption>FIG. 3. (a)–(d) Poincaré sections at time $T/8$ for the drive frequencies (a) $\tilde{\omega}_d = 0.75$, (b) $\tilde{\omega}_d = 1.25$, (c) $\tilde{\omega}_d = 2.25$, and (d) $\tilde{\omega}_d = 4.25$. The drive amplitudes are chosen such that the frequency of the $0-1$ transition in the quantum system is ac-Stark shifted by 100 MHz for all $\tilde{\omega}_d$. (e), (f) Mean energy spectra as a function of $\tilde{\epsilon}_d/\tilde{\omega}_d$ at the same respective drive frequencies, with $\hbar_{\text{eff}}^{-1} = 2.45$ and $n_g = 0.25$. In each case, the drive amplitude yielding an ac-Stark shift of 100 MHz is indicated by the vertical black dashed line—the above Poincaré sections have been computed for these amplitudes. For the parameters of (a), (c), (d), a 1:1 resonance occurs between the drive and the system. In (a), this resonance results in a second set of tori within the bounded states of the system (red region). The overlap of the small oscillation and this resonance results in a large chaotic layer. In the quantum system (e), this leads to strong level hybridization. In (c), (d), the drive comes in resonance with the unbounded states, resulting in two additional out-of-phase sets of tori rather than one (red regions). From (c), (d), these tori move away from the center as the drive frequency increases. In (e), (f), 3:1 and 5:1 resonances emerge (green region). In (g), the proximity of the 1:1 (red) and 3:1 (green) resonances to the regular island also causes large instabilities, both in the classical and quantum systems. In (d), the resonances are far and the width of the chaotic layer is smaller, as expected from Eq. (5). In (b), the resonance is absent.</figcaption></figure>

tori, see the red-colored region in Fig. 3(a). Large instability results from the overlap between these two tori, with a chaotic layer arising at their separatrices [35]. The resulting increased width of the chaotic layer is not captured by the approximate expression of Eq. (5). The case $k=1$ leads to a 3:1 resonance for $\tilde{\omega}_d = 3\tilde{\omega} \in [2, 3]$. For $\tilde{\omega}_d = 2.25$ this resonance results in three sets of tori surrounding the central regular island, see the three green-colored regions in Fig. 3(c). Although of smaller amplitudes than the 1:1 resonance, the overlap of this resonance with the main set of tori is likely to produce unstable motion at the boundary of these. Note also the appearance of a weak 5:1 resonance at $\tilde{\omega}_d = 4.25$ resulting in five small tori within the regular island, see the green-colored region in Fig. 3(d).

In contrast, no direct resonance occurs for the case $1.1 < \tilde{\omega}_d < 1.5$, which is common for the dispersive readout, see

Fig. 3(b) for $\tilde{\omega}_d = 1.25$. Nevertheless, as expected from Eq. (5), the width of the chaotic layer is large at this drive frequency.

In the quantum case, the resulting large chaotic layer translates into a large hybridization of the states even at low drive amplitude, see the mean energy spectra Figs. 3(e)–3(g). This hybridization can be further observed from the modification of the rate matrices (see Sec. A). This state hybridization can lead to loss of the QND character of the dispersive readout. As a concrete example, the value of $\hbar_{\text{eff}}^{-1}$ and $\tilde{\omega}_d = 0.75$ of Fig. 3(f) is chosen to match the experimental parameters of Ref. [28] where the dispersive readout fidelity was observed to degrade for measurement photon numbers $\bar{n} > 2.5$, something which was attributed to measurement-induced mixing of unknown origin. Using the light-matter coupling of
$g/2\pi = 208$ MHz reported in Ref. [28], $\bar{n} = 2.5$ can be converted to an effective drive $\tilde{\varepsilon}_d = 0.105$ on the transmon. As can be observed in Fig. 3(e), the first excited state is “absorbed” in the chaotic layer for $\varepsilon_d \gtrsim 0.1$. Since chaotic states often lead to strong hybridization between the transmon and the readout resonator (see Sec. V), this hints at chaos-induced state mixing and non-quantum-demolition effects in the readout beyond that drive amplitude. Although the above mean-field analysis can qualitatively predict unstable behavior of the transmon, the agreement cannot be expected to be quantitative at low resonator photon numbers since vacuum fluctuations of the resonator field and qubit-resonator parametric processes occurring in the presence of a readout tone can play an important role. We address some of these mechanisms in Sec. V.

The drive can also come in resonance with unbounded states of the pendulum with energy $\tilde{H}_{\tilde{\varepsilon}_d} > -1 + 2J_0(\tilde{\varepsilon}_d/\tilde{\omega}_d)$ ~ 1. The trajectories at those energies take the approximate form $\tilde{\phi}(\tilde{t}) = \pm \tilde{\omega} \tilde{t} + F(\tilde{t})$, where $F$ is a function of period $2\pi/\tilde{\omega}$ [35]. Using this expression in Hamilton's equation of motion for the phase, we find that the oscillation frequency satisfies $\tilde{\omega} = |\langle\tilde{n} - \tilde{n}_g\rangle| \gtrsim 1.5$, where $|\langle\tilde{n} - \tilde{n}_g\rangle|$ is the averaged momentum of the trajectory over one period. This lower bound can vary with the drive amplitude through the effect of the reduction factor $J_0(\tilde{\varepsilon}_d/\tilde{\omega}_d)$. To leading order, substituting this expression for $\tilde{\phi}(\tilde{t})$ in the second line of Eq. (7) results in the term $A_0 J_1(\tilde{\varepsilon}_d/\tilde{\omega}_d) \cos[\pm(\tilde{\omega} - \tilde{\omega}_d)\tilde{t}]$, where $A_0$ is the zeroth Fourier component of $\cos[F(\tilde{t})]$. Hence, for $\tilde{\omega}_d = \tilde{\omega} \gtrsim 1.5$, this resonance appears in the form of two tori moving clockwise and anticlockwise, corresponding to $\pm \tilde{\omega}_d$, as depicted by the red-colored regions in Fig. 3(c) for $\tilde{\omega}_d = 2.25$ and Fig. 3(d) for $\tilde{\omega}_d = 4.25$. At $\tilde{\omega}_d = 2.25$, the proximity of the resonance with the regular island results in a large chaotic domain, which translates in strong state hybridization in the quantum case, see Fig. 3(g). Because the average momentum associated to this resonance is $|\langle\tilde{n} - \tilde{n}_g\rangle| = \tilde{\omega}_d$, the pair of tori move further away from the center of phase space with increasing drive frequency. For this reason, at $\tilde{\omega}_d = 4.25$, this resonance is far from the separatrix of the undriven system and does not affect the chaotic layer. As a result, for $\tilde{\omega}_d = 4.25$ and at larger drive frequencies, the width of the chaotic layer is exponentially suppressed as expected from Eq. (5). In this situation, the absence of instability yields a mean energy spectrum with little state hybridization, see Fig. 3(h). This fact is further illustrated by the regular structure of the charge matrix elements in Sec. A.

To summarize the above, the frequency ranges $0.65 < \tilde{\omega}_d < 1$ and $1.5 \le \tilde{\omega}_d \le 3$ lead to strong instabilities, in particular for the former. The parameter regime $1.1 < \tilde{\omega}_d < 1.5$ avoids resonances but has a large chaotic layer around the separatrix. For $\tilde{\omega}_d \ge 3.5$, resonances do not affect the bounded states (or very weakly through

resonances of order greater than 5), and the width of the chaotic layer is suppressed. As shown in the next sections, the presence of the chaotic layer impacts the coherence times of the transmon qubit, and it should therefore be minimized when operating the transmon with strong drive, e.g., in dispersive qubit readout. The frequency ranges $1.1 \le \tilde{\omega}_d \le 1.5$ and $\tilde{\omega}_d \ge 3.5$ seem to be more benign.

## IV. IMPACT ON COHERENCE PROPERTIES

Having established that the driven capacitively shunted Josephson junction exhibits signatures of chaos even at the large $\hbar_{eff}$ corresponding to the transmon regime, we now turn to the impacts of this observation on coherence properties of the transmon in the presence of a bosonic bath. In this situation, the total Hamiltonian now takes the form

$$ \mathbf{H}(t) = 4E_C(\mathbf{n} - n_g)^2 - E_J \cos(\phi) + \varepsilon_d \cos(\omega_d t)\mathbf{n} + i\mathbf{n} \sum_k g_k(\mathbf{b}_k^\dagger - \mathbf{b}_k) + \sum_k \omega_k \mathbf{b}_k^\dagger \mathbf{b}_k, \quad (8) $$

where $\mathbf{b}_k$ and $\mathbf{b}_k^\dagger$ are the annihilation and creation operators of the bath modes. Because of hybridization with states in the chaotic layer, which have a strong charge dispersion, it is important to keep the gate charge $n_g$ in Eq. (8) even when interested in the coherence properties of the low-lying eigenstates of the system.

### A. Rate matrix and steady-state population

Within the Floquet-Markov description of the system-bath coupling, the bath-induced transition rate from Floquet state j to i is given by [34]

$$ \Gamma_{ij} = \sum_k |n_{ijk}|^2 [\Theta(\Delta_{ijk}) + n_B(|\Delta_{ijk}|)] J(|\Delta_{ijk}|), \quad (9) $$

where $\Theta(x)$ denotes the Heaviside function, $n_B(x)$ is the thermal occupation number of the bath, and $\Delta_{ijk} = \varepsilon_j - \varepsilon_i - k\omega_d$. The spectral function $J(x)$ is assumed to that of an Ohmic bath, i.e., $J(x) \propto x \exp(-|x|/\omega_c)$, where $\omega_c$ is a high-frequency cutoff. In the expression for $\Gamma_{ij}$, k can be interpreted as the number of drive photons participating positively or negatively to the transition. We also introduce the charge operator matrix elements

$$ n_{ijk} = \frac{\omega_d}{2\pi} \int_0^{2\pi/\omega_d} dt \langle \phi_i(t) | \mathbf{n} | \phi_j(t) \rangle \exp(i\Delta_{ijk}t), \quad (10) $$

where $|\phi_j(t)\rangle$ are the Floquet modes of the driven system.

The rate matrices of the driven transmon at $\tilde{\varepsilon}_d = 0.4$ and a temperature of $T = 10$ mK are shown for $n_g = 0.5$ in Fig. 4(a) and for $n_g = 0.25$ in Fig. 4(b). As a comparison, the insets show the rate matrices in the undriven case. For a given transition $\Gamma_{ij}$, the blue (red) squares sum
<figure><img src="image_6.png" /><figcaption>FIG. 4. Top panels: rate matrices in units of the coupling strength to the bath, at $\tilde{\epsilon}_d = 0.4$, $\hbar_{eff}^{-1} = 3$, $\tilde{\omega}_d = 1.34$, and $T = 10$ mK for (a) $n_g = 0.5$ and (b) $n_g = 0.25$. The states are sorted by their mean energy $\langle\langle H \rangle\rangle$. The blue (red) squares correspond to the sum of rates involving an even (odd) number of drive photons $k$. The insets show the rate matrices of the undriven systems. At finite drive, the instability develops around the separatrix leading to an irregular block in the rate matrix also allowing for upward excitation rates. At the symmetric point $n_g = 0.5$, the red and blue squares do not overlap, while at $n_g = 0.25$, hybridization of parity sectors can be observed through the presence of purple squares even in the low-energy sector. The purple squares in (a) are due to numerical errors (see text). Bottom panels: steady-state population as derived from the above rate matrices for (c) $n_g = 0.5$ and (d) $n_g = 0.25$. In the driven system, a plateau forms over the states that are part of the chaotic block of the rate matrices, instead of the regular exponential distribution in the undriven case, as shown in the insets.</figcaption></figure>

the contributions from even (odd) values of $k$, with purple squares indicating contributions from both even and odd values. At zero drive, the upper triangular sector of the rate matrices (corresponding to upward transitions) contains negligible but nonzero elements due to the finite temperature (not visible in the insets). At finite drive, the instability develops around the states located on the separatrix (typically around the eighth excited state), forming an irregular block in the rate matrix. In particular, states within the chaotic layer are all coupled to one another through the charge operator. In addition, because of the drive photons, upward transitions are now apparent. The appearance of an irregular block in the rate matrix directly relates to the repulsive statistics of the quasienergies. In fact, chaotic systems can be accurately described by random matrices in the limit of a large number of chaotic states [38]. In Sec. E, we leverage this property to estimate the average charge matrix element between two chaotic states when all the low-energy states are chaotic.

In the transmon regime, the low-energy sector is almost independent of the offset charge and, if one neglects its influence, the Hamiltonian becomes effectively symmetric under the parity transformation $\phi \to -\phi$ and $n \to -n$ (see Sec. B). Because it neglects the gate charge, this symmetry is implicit in the Kerr nonlinear oscillator model of the transmon. Although this symmetry is exact only at $n_g = 0$ and $n_g = 0.5$, at zero drive it results in a suppression of the matrix elements of the charge operator $n_{i,i+2}$ in the low-energy sectors for all values of $n_g$, and forbids the transition $i \to i+2$, see insets of Figs. 4(a) and 4(b). The fact that this is only an approximate symmetry at $n_g = 0.25$ is apparent for the states in the separatrix region and above.

In the presence of a drive term $\epsilon_d \cos(\omega_d t)n$, the inversion symmetry only holds together with the time translation $t \to t + \pi/\omega_d$ [34]. This symmetry of the driven Hamiltonian defines even and odd parity sectors among the time-dependent Floquet states (see Sec. B). Because the charge operator is antisymmetric, under this generalized parity symmetry, transitions through the charge operator between two states of the same (opposite) parity can only involve an odd (even) number of drive photons $k$. As can be seen in Fig. 4(a) for $n_g = 0.5$, blue and red squares do not mix (i.e., there are no purple squares) indicating that the symmetry is respected for that gate charge. The purple squares appearing for states 18 and 19 only result from numerical precision errors. The situation is very different at $n_g = 0.25$ where purple squares appear in the chaotic layer but also for transitions involving the low-energy states. This results in transitions that are otherwise forbidden at the symmetric points $n_g = 0$ and $n_g = 0.5$. As discussed in further details in the next section, the breaking of this effective symmetry in the low-energy sector is a consequence of a strong increase of the band dispersion in the presence of drive. Interestingly, transitions forbidden by the apparent inversion symmetry of the transmon were experimentally observed under strong drives, but remained unexplained [23].

The rate matrix can also be used to compute the system's steady-state density matrix $\rho_{ss}$. Under the assumption of weak system-bath coupling, the steady state is diagonal in the Floquet basis

$$ \rho_{ss}(t) = \sum_j p_j |\phi_j(t)\rangle \langle\phi_j(t)|, \quad (11) $$

where the populations $p_j$ satisfy the rate equations $p_j = \sum_i \Gamma_{ji} p_i - \sum_i \Gamma_{ij} p_j$. The insets of Figs. 4(c) and 4(d) show these steady-state populations for the undriven systems for $\hbar_{eff}^{-1} = 3$ which is typical of the transmon regime. As expected, the populations follow the thermal distribution with populations quickly dropping below $10^{-10}$. In contrast, in the driven systems (here with $\tilde{\epsilon}_d = 0.4$), the steady-state populations form a plateau corresponding to
the chaotic layer. This behavior is typical of chaotic systems [32,33]. Increasing the drive amplitude, the chaotic layer and therefore the plateau grow until all the low-energy states are part of the plateau, including the ground state. This results in a dramatic decrease of the purity of the transmon's steady state, an observation, which is in agreement with numerical [26] and experimental results [13].

This discussion sheds light on the distinction between the time dynamics of the transmon ionization numerically observed in Ref. [27], and the ionization in the steady state [13,26]. In the former, one captures the dynamics of the dispersive readout of a transmon qubit. As the cavity rings up on a timescale $\kappa^{-1}$, where $\kappa$ is the photon loss rate, the effective field amplitude $\epsilon_d$ on the transmon increases. Starting in the ground or first excited states, the system follows the corresponding first or second line in the mean energy spectrum of Fig. 1(c). Ionization occurs when one of those lines crosses a large resonance with a chaotic state, i.e., when the mean energy suddenly increases in Fig. 1(c), something which happens at $\tilde{\epsilon}_d \approx 1.1$ for the ground state and at $\tilde{\epsilon}_d \approx 0.75$ for the first excited state. Once a chaotic state is populated, it decays either through the transmon-bath coupling or through the transmon-resonator coupling due to its numerous possible transitions at various frequencies (see Sec. V). On the other hand, ionization in the steady state is only a function of the rate matrix, which in turn depends on the matrix elements of the charge operator in the Floquet basis. Large matrix elements between the ground state and the chaotic states are likely to lead to ionization through bath-induced transitions. Although a resonance leads to large matrix elements, the contrary is not necessarily true.

## B. Chaos-assisted tunneling: effects on $T_2$

A defining feature of the transmon is its charge dispersion, which is exponentially suppressed with $8\hbar^{-1}_{\text{eff}} = \sqrt{8E_J/E_C}$, making it almost insensitive to charge noise [1]. Here, we show that this exponential suppression does not necessarily hold in the presence of a periodic drive.

In the previous section, we have seen that the Hamiltonian of the driven transmon can be approximately divided into regular blocks of states: phaselike states at the bottom of the cosine potential well; chargelike states above that potential, on which the drive acts pertubatively; and one chaotic block with a strongly correlated spectrum. From perturbation theory performed on the regular block, one would expect a drive to lead to a slow hybridization amongst regular low-lying states. In this situation, the ground and first excited states would therefore weakly inherit an offset-charge sensitivity from weak dressing with higher energy states, leading to an overall small increase of the band dispersion with the drive amplitude. In practice, we find that the presence of the chaotic layer results in the energy dispersion of these states to be significantly

<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) Band dispersion of the first excited state for the undriven (dashed lines) and driven (solid lines) transmon with $\tilde{\epsilon}_d = 0.4$, for $\hbar_{\text{eff}}^{-1} = \sqrt{E_J/8E_C}$ in the range [1, 7]. The experimentally relevant value of $\hbar_{\text{eff}}^{-1} = 3$ is highlighted in red. The exponential suppression of the dispersion of the undriven case is strongly modified by the presence of the drive, showing both an overall increase and sharp features. The latter are due to small resonances with chaotic states that have a large charge dispersion. (b) Fourier coefficients of the energy $\epsilon_1(n_g)$ of the first excited state at finite drive amplitude $\tilde{\epsilon}_d = 0.4$ as a function of the Fourier index $n$. These Fourier components are understood as phase-slip rates of order $n$. The spikes or “abrupt” resonances in (a) are responsible for the slowly decaying tail in (b) that translates to long-range tunneling between the Josephson potential wells. The inset shows the tunneling rates of the undriven systems with the expected exponential suppression with the distance between wells.</figcaption></figure>

modified when the system is driven. This is illustrated by the solid lines in Fig. 5(a) for the first excited state of the transmon with values of $\hbar_{\text{eff}}$ in the range 1 to 7 as labeled by the different colors. As a comparison, the dashed lines correspond to the undriven transmon and for which the exponential suppression of the charge dispersion is clearly observed. In contrast, in the driven transmon the energy bands are disrupted by peaks of multiple orders of magnitude in addition to being, on average, substantially larger than in the undriven case. The sharpness of these peaks is indicative of a weak resonance between the first excited state and strongly $n_g$-dependent states located in the chaotic layer.

These results are obtained by identifying the first excited state $|\phi_1^0\rangle$ of the undriven system, then tracking the corresponding Floquet mode $|\phi_1(t)\rangle$ as a function of the drive amplitude $\tilde{\epsilon}_d$. More precisely, $|\phi_1(t)\rangle$ is the Floquet state maximizing the overlap with $|\phi_1^0\rangle$ at time 0. As the maximum overlap decreases with drive amplitude
due to state hybridization, this method is likely to fail at high drive amplitude, typically for $\tilde{\epsilon}_d \ge 0.4$. To circumvent this problem, at every increase of 0.1 in the drive amplitude, we update the reference state to which the overlap is computed to the previously identified state, hence progressively accounting for state hybridization. The same method is used to track the ground state.

This phenomenon is closely related to chaos-assisted tunneling (CAT) [29]. In CAT, tunneling between two sets of disjoint regular states is facilitated by their coupling to delocalized states in the chaotic layer. Moreover, because of the participation of states within the chaotic layer, the tunneling rates are expected to vary widely with the control parameters [29]. For the transmon, chaos can assist tunneling between different wells of the cosine potential. Because the phase of the transmon is compact, tunneling between wells distant by $\delta\phi = 2\pi n$ translates to $n 2\pi$ swings of the transmon phase or, equivalently, to quantum phase slips of order $n$. The transmon states acquire their $n_g$ dependence through these full phase rotations [1]. In the undriven transmon, the rate of these events decreases exponentially with $n$ for large $\tilde{h}_{\text{eff}}^{-1}$.

To evaluate the phase-slip rate in the presence of a drive, we compute the Fourier transform of the energy bands of the system [48]. Indeed, the Fourier components $t_n = \int_{-0.5}^{0.5} dn_g \epsilon_1(n_g) e^{i2\pi n_g n}$ of the energy $\epsilon_1(n_g)$ of the first excited state corresponds to the rate of phase slips of order $n$ when the system is in the first excited state [49]. The components $t_n$ are plotted as a function of the index $n$ in Fig. 5(b). The inset shows the exponential suppression of $t_n$ expected for the undriven case. In the driven case (main panel), $t_n$ is no longer exponentially suppressed with $n$. Instead, it shows a long tail, indicating that the sharp features seen in Fig. 5(a) results from long-range hopping between the wells. Note that the tunneling rates $t_n$ remain globally suppressed with $\tilde{h}_{\text{eff}}^{-1}$, see Fig. 5(b). This is explained by an exponential suppression of the matrix elements $n_{ijk}$ between the regular states and the chaotic states with $\tilde{h}_{\text{eff}}^{-1}$. We discuss this point in Sec. VI.

In cold atoms trapped in a driven optical lattice, the signatures of CAT have been observed by measuring coherent oscillations between states localized in distinct wells [50,51]. Because the phase coordinate of the transmon is compact, tunneling oscillation between wells cannot be measured. However, phase slips (due to CAT or not) lead to a phase accumulation, which depends on the gate charge $n_g$ due to the Aharonov-Casher effect [12]. Because the gate charge is a fluctuating function of time, this results in enhanced dephasing of the transmon. The increased phase-slip rate due to CAT can thus, in principle, be witnessed through sharp variations of the pure dephasing rate $\gamma_\phi$ associated to the driven transmons logical states $|\phi_0(t)\rangle$ and $|\phi_1(t)\rangle$, as a function of $n_g$ and of the drive amplitude $\tilde{\epsilon}_d$.

Within the Floquet-Markov theory, the dephasing rates takes the form [52]

$$ \gamma_{\phi} = A_e |2g_{0,\phi}| \sqrt{|\log \omega_{IR} t_m|} + \sum_{k \neq 0} 2S(k\omega_d) |g_{k,\phi}|^2. \quad (12) $$

The first term represents $1/f$ charge noise, while the second term comes from the possible conversion of a photon loss to a dephasing event due to the hybridization of the logical undriven states $|0\rangle$ and $|1\rangle$. In the first term, $A_e$ is the amplitude of charge noise, $g_{k,\phi} = n_{11k} - n_{00k}$ where $n_{ijk}$ is a matrix element of the charge operator defined in Eq. (10), $\omega_{IR}$ is the infrared cut-off, and $t_m$ is the characteristic measurement time. Typical values are $\sqrt{|\log \omega_{IR} t_m|} \sim 4$ and $A_e = 10^{-4} e$ [53]. In the second term, $S(\omega)$ is the spectral function of the relevant bath, here assumed to be Ohmic. Its amplitude is chosen such that the lifetime of the undriven transmon is 50 µs.

The dephasing rate obtained from Eq. (12) and the numerically obtained Floquet modes $|\phi_0(t)\rangle$ and $|\phi_1(t)\rangle$ corresponding to the logical states in the presence of a drive is shown in Fig. 6(b) as a function of $\tilde{\epsilon}_d$ for 50 values of $n_g$ uniformly spaced in the range [0, 0.5] (red lines). All curves show a slow quadratic increase of $\gamma_\phi$ with the drive amplitude due to dielectric losses through the hybridization of $|0\rangle$ and $|1\rangle$ by the drive. This quadratic increase in the rate with the amplitude of the drive is a signature of a perturbative effect. More interestingly, the dephasing rate also displays sharp peaks whose position strongly depends on $n_g$ and $\tilde{\epsilon}_d$, as is expected for CAT. To better understand the origin of these structures, the mean energy of the different Floquet mode is plotted in Fig. 6(a) for $n_g = 0.13$ as a representative example. The value of $\gamma_\phi$ obtained for the same gate charge is highlighted in Fig. 6(b) (blue line). By comparing the two plots (see the vertical dashed lines), it becomes clear that the sharp increases in $\gamma_\phi$ correspond to resonances between regular states (ground or first excited state) and chaotic states. As mentioned above, the resulting hybridization of the computational states with states that have a strong charge dispersion leads to a sharp increase in the dephasing rate. The black line is an average over all realizations of the gate charges. Depending on the time scale of the charge fluctuations [54,55] and of time needed to measure $\gamma_\phi$, this average may be more representative of potential experimental observations.

In summary, in this section we have first shown that chaotic states form an irregular block in the rate matrix of a driven transmon coupled to a bath. The growth of the chaotic layer with drive results in the breaking of the effective parity symmetry of the low-energy states. The existence of the chaotic states also strongly impacts the steady-state population. Secondly, we find that the band dispersion of the ground and first excited states of the transmon is dramatically modified under drive due to chaos-assisted quantum phase slips. Due to $1/f$ charge
<figure><img src="image_8.png" /><figcaption>FIG. 6. (a) Mean energies of Floquet states of the driven transmon as a function of the drive amplitude, with $\hbar_{\text{eff}}^{-1} = 3$, $\tilde{\omega}_d = 1.34$, and $n_g = 0.13$. While the ground and first excited states of the transmon are tracked as we vary the drive amplitude (blue lines), the rest of the states are not tracked (black dots). (b) Dephasing rate of the transmon as a function of drive amplitude due to $1/f$ charge noise and dielectric losses (see text) for 50 values of $n_g$ distributed uniformly between 0 and 0.5 (red lines) and the same parameters used in (a). The blue line corresponds to $n_g = 0.13$, as in (a). Small resonances between the computational states and the chaotic states that can be observed in (a) result in sharp peaks in the dephasing rate, due to $1/f$ charge noise. The average dephasing rate over the offset charge (black line) is dominated by $1/f$ charge noise.</figcaption></figure>

noise, these processes impact the coherence time of the transmon.

## V. CIRCUIT QED: TRANSMON COUPLED TO A RESONATOR

We have so far treated the drive seen by the qubit as a purely classical field. To account for vacuum fluctuations and the richer structure of the energy levels in the presence of a cavity, we now consider a circuit QED setup where the transmon is capacitively coupled to a cavity. In this situation, the Hamiltonian takes the form [40]

$$ \mathbf{H}(t) = 4E_C(\mathbf{n} - n_g)^2 - E_J \cos(\phi) + \varepsilon_d \cos(\omega_d t)\mathbf{n} \\ + \omega_a \mathbf{a}^\dagger \mathbf{a} - i g \mathbf{n}(\mathbf{a} - \mathbf{a}^\dagger) \\ - (\mathbf{a} - \mathbf{a}^\dagger) \sum_k g_k (\mathbf{b}_k - \mathbf{b}_k^\dagger) + \sum_k \omega_k \mathbf{b}_k^\dagger \mathbf{b}_k. \quad (13) $$

The first line of this expression is the Hamiltonian of the driven transmon qubit, as in Eq. (2). The second line contains the free cavity Hamiltonian defined by a frequency $\omega_a$, together with the transmon-cavity capacitive coupling $g$. We take $g/2\pi = 250$ MHz and $\omega_a/2\pi = 8$ GHz, with the drive detuned from the cavity at $\omega_d/2\pi = 7.5$ GHz, while keeping the same parameters for the transmon qubit as in the previous sections. Finally, the last line represents the capacitive coupling of the cavity to a bosonic bath. The drive term on the qubit can be either seen as directly acting on the qubit in the laboratory frame, or as resulting from a drive on the resonator. In that latter situation, the Hamiltonian of Eq. (13) is to be understood as expressed in a displaced frame where the drive on the cavity has been removed, i.e., $\varepsilon_d = 2g\sqrt{\tilde{n}}$ where $\tilde{n}$ is the cavity steady-state photon number in the laboratory frame. In this section, we find the Floquet spectrum of Eq. (13) in the absence of coupling to the bath, and then calculate transition rates in linear response theory [26]. Details on the Floquet simulations are provided in Sec. C. In the following, the gate charge is set to $n_g = 0.25$. We have numerically checked that the phenomena observed here are qualitatively the same for other values of the offset charge.

## A. Structure of the Floquet spectrum

To characterize the structure of the Floquet spectrum, for each Floquet mode we compute expectation values of a pair of operators that are good quantum numbers in the undriven and decoupled Hamiltonian: the transmon excitation number $N_t = \sum_{i,n} i |in\rangle\langle in|$ and the resonator excitation number $N_r = \sum_{i,n} n |in\rangle\langle in|$, where $|in\rangle$ is a bare state of the joint transmon-cavity system. In Fig. 7 we show these quantities on a two-dimensional grid in the $((\langle\langle N_t \rangle\rangle), (\langle\langle N_r \rangle\rangle))$ plane, the (a)–(d) representing different drive amplitudes $\tilde{\varepsilon}_d$. As in Eq. (4), the double angle brackets represent the time-averaged expectation value of the operator in a given Floquet mode over a period of the drive. In the absence of coupling and drive, this grid is expected to be rectangular. In the coupled case, for high enough resonator photon number $\langle\langle N_r \rangle\rangle$ or for large enough drive power, some Floquet modes deviate strongly from the regular rectangular grid. These strong deviations are associated with state dressing, and a significant hybridization of the two subsystems. More precisely, horizontal deviations from the rectangular grid can be interpreted as hybridization of transmon states, while vertical deviations originate from hybridization of Fock states. The hybridization of Fock states usually occurs through entanglement of the two systems, implying also a hybridization of the transmon states.

The degree of hybridization between the two subsystems is measured by the purity of the transmon reduced density matrix. Encoding this purity into the color of the grid points in Fig. 7, deviations from the rectangular grid appear
<figure><img src="image_9.png" /><figcaption>FIG. 7. Grid of Floquet modes in the $(\langle N_r \rangle_t, \langle N_r \rangle_t)$ plane. The purity of the transmon reduced density matrix is encoded in the symbol color, for four different drive amplitudes: $\tilde{\epsilon}_d = 0(a)$, 0.25(b), 0.5(c), and 0.95(d). In a zero-drive decoupled system, the points should form a square grid. In the coupled system, as $\langle N_r \rangle$ and $\tilde{\epsilon}_d$ increase, level hybridization leads to grid-point displacement, and purity drops. With increasing drive, for $\langle N_r \rangle_t \approx 0$, there is a purity drop in Floquet modes corresponding to the chaotic layer. Steady state of Floquet-Markov master equation is represented by a red cross. The thickness of each arrow between a pair of Floquet modes is proportional to the corresponding rate $\Gamma_{ij}$ at zero temperature, apart from an offset to render the smallest rates visible. An arrow is plotted if $\Gamma_{ij} > 10^{-3}\kappa$ for (a),(b), or $10^{-2}\kappa$ for (c),(d). Rates of the undriven system are dominated by single-photon relaxation (vertical lines), and Purcell decay (left-to-right arrows for $\langle N_t \rangle_t \lesssim 6$) and rate flows tend to increase purity (a). Rates become increasingly nonlocal as the drive strength is increased. (b)–(d) Steady-state populations of each Floquet mode are encoded in the radii of circles around each grid point. For sufficiently low drive, the steady state is dominated by the vacuum state (a)–(c), unlike steady states above the ionization threshold (d). See Fig. 14 for the same data at higher drive power.</figcaption></figure>

to correlate with drops in purity, i.e., an increase of entanglement entropy between the transmon and the resonator. In particular, purity drops are drastic for states whose transmon excitation number corresponds to the chaotic layer identified in the previous sections for similar drive powers.

We can gain further understanding of the purity drop by first diagonalizing the time-dependent Hamiltonian of the driven transmon [first line of Eq. (13)], and then expressing the transmon-resonator coupling $-ig\mathbf{n}(\mathbf{a} - \mathbf{a}^\dagger)$ in the joint basis of the transmon Floquet states $\{|\tilde{i}\rangle\}$ and cavity Fock states. In a frame rotating at the drive frequency for

the resonator and at the quasienergies for the transmon, the Hamiltonian reduces to

$$ \mathbf{H}(t) = \Delta \mathbf{a}^\dagger \mathbf{a} - ig \sum_{i,j,k} n_{ijk} |\tilde{i}\rangle \langle \tilde{j} | (e^{i\Delta_{ij,k+1}t} \mathbf{a} - e^{i\Delta_{ij,k-1}t} \mathbf{a}^\dagger), \quad (14) $$

where $\Delta = \omega_a - \omega_d$. The charge-matrix elements $n_{ijk}$ and the energy differences $\Delta_{ijk}$ are defined in Sec. IV, and the Floquet modes of the transmon are evaluated at $t = 0$. Large hybridization occurs between the states $|\tilde{i}, n\rangle$ and $|\tilde{j}, n+1\rangle$ if the coupling strength is of the order of the transition frequency, i.e., $2g\sqrt{n+1}n_{ijk} \sim \Delta_{ij,k-1} - \Delta$. Note that, for small detuning $\Delta$, the transition frequencies $\Delta_{ij,k-1} - \Delta$ are small for $k = 0, 1$, and consequently hybridization of the states is mainly due to the matrix elements $n_{ij,k=0}$ and $n_{ij,k=1}$. For chaotic states, we derive in Sec. E an estimate of the matrix elements $n_{ijk}$ and the frequencies $\Delta_{ij,k-1}$. Using these estimates and the previous expression allows us to obtain a threshold value on the photon number $n$ for strong transmon-resonator hybridization (see details in Sec. F). For the parameters of Fig. 7, we find that the transmon states become strongly coupled even at the lowest Fock states. This strong hybridization can result in large decay rates for chaotic states, as explained in the paragraph below.

Transition rates, defined within linear response theory for the charge operator of the resonator $-i(\mathbf{a} - \mathbf{a}^\dagger)/\sqrt{2}$ in a similar way as in Eq. (9), are shown as gray arrows connecting the grid points in Fig. 7. For simplicity, we take in this section a flat spectral function for the bath coupled to the resonator, i.e., $J(\omega) \propto \Theta(\omega)$. At low drive power, rates are predominantly local, consisting of single-photon relaxation (vertical downward-pointing arrows connecting neighboring grid points), and qubit Purcell decay (horizontal left-pointing arrows also connected neighboring grid points). The existence of one dominant rate allows identification of states corresponding to definite transmon excitation number [27,56]. In contrast, in the chaotic layer where the purity is low, at finite drive amplitude, the rate matrix becomes nonlocal connecting nonadjacent points. Note that here, dissipation on the resonator is treated in linear response theory. In Sec. H, we verify that typical loss rates of readout resonators ($\kappa/2\pi \sim 10$ MHz) do not alter the spectrum of the undriven system. In particular, hybridization is not weakened by the linewidth of the Fock states considered in Fig. 7.

Using the transitions rates, we can now find the steady state of the Floquet-Markov master equation [26]. In order to avoid the saturation of the resonator Hilbert space (see Sec. C), we impose a cutoff by setting all rates involving states with $\langle\langle N_r \rangle\rangle \ge 15$ to vanish. We find that this causes only a small quantitative change in the steady-state density matrix. The transmon and resonator excitation numbers in
this steady state are represented on Fig. 7 as a red cross for each drive power. Moreover, the occupation $p_i$ of each Floquet mode in the steady state is encoded in the area of a circle centered at each grid point. For sufficiently low drive power, the steady state is dominated by the state with $\langle\langle N_{r,t} \rangle\rangle \approx 0$, i.e., the dressed-state closest to the vacuum state. Beyond a threshold drive amplitude, the steady state of the Floquet-Markov Lindblad master equation has significant weights on the chaotic states (see Fig. 14 of Appendix C, which shows similar results as Fig. 7 but for larger drive amplitudes). This threshold corresponds to transmon ionization where states above the cosine potential of the transmon become populated [13,26,27]. Notably, strong hybridization leads to nonzero resonator photon number in the steady state, that is $\langle\langle N_r \rangle\rangle \gtrsim 2$ at $\tilde{\epsilon}_d = 0.95$, see Fig. 13.

Comparing to the mean-field study of a driven transmon of the previous sections, we find that the qualitative features of the spectrum of the off-resonantly driven circuit QED Hamiltonian match those of the driven transmon qubit. To illustrate this, we plot in Fig. 8(a) the transmon excitation number $\langle\langle N_t \rangle\rangle$ for the full circuit QED model of Eq. (13) (colored dots) versus drive amplitude; to avoid overcrowding the plot, we retain only those Floquet modes for which $\langle\langle N_r \rangle\rangle \approx 0$ (in practice, we retain only modes with $\langle\langle N_r \rangle\rangle < 0.73$). We compare this excitation number to the corresponding observable in the Floquet spectrum of Eq. (2) with the same parameters, but corresponding to a transmon driven by a classical field (black crosses). As in Fig. 7, for the former, color encodes the purity of the transmon-reduced density matrix. While the agreement between the full circuit QED simulation and the mean field driven transmon is excellent for regular states, deviations appear within the chaotic layer where the purity drops. Differences within the chaotic layer between the two models are expected due to strong hybridization between the transmon and the resonator (not accounted in the mean field model), and due to fine sensitivity of chaotic spectra on external parameters [38]. Nevertheless, in the two models, the threshold values of $\tilde{\epsilon}_d$ for the onset of chaos agree.

## B. Frequency response of the resonator

Performing the analog of a numerical two-tone spectroscopy experiment, in Fig. 8(b) we plot the pulled frequency response of the cavity versus drive amplitude $\tilde{\epsilon}_d$. To obtain these results, we first identify the cavity vacuum states corresponding to each transmon occupation number, i.e., $\langle\langle N_r \rangle\rangle \approx 0$. Then, for each of these vacuum states, labeled $i$, we identify the dominant transition corresponding to the largest matrix element $|n_{ijk}|$; with $j, k$ thus identified, we plot a point whose $x$ coordinate is $\tilde{\epsilon}_d$ and whose $y$ coordinate is the difference between $\Delta_{ijk}$ and the bare cavity frequency mod $\omega_d$. The corresponding transition rate is encoded in the symbol size. As in Fig. 8(a), the symbol color encodes the purity of the corresponding

vacuum state. Only states with purity above 0.85 are plotted, for which the cavity vacuum can be well identified, since the states are close to a tensor product state, so the condition $\langle\langle N_r \rangle\rangle \approx 0$ is easily fulfilled.

Figure 8(b) shows that, as expected, low-energy regular states of the transmon correspond to a traceable resonator pull that is monotonically decreasing with drive strength [27,40,56]. In particular, there is a group of four levels (corresponding to pulls up to slightly below 30 MHz at zero drive), which exhibit ac-Stark shifts monotonically decreasing with drive strength $\tilde{\epsilon}_d$ towards the bare frequency of the cavity, for sufficiently low drive. These levels correspond to resonator states pulled by the lowest

<figure><img src="image_10.png" /><figcaption>FIG. 8. (a) Agreement of transmon population $\langle\langle N_t \rangle\rangle$ between the full circuit QED simulation (color encodes purity of transmon-reduced density matrix; only states satisfying $\langle\langle N_r \rangle\rangle \le 0.73$ are shown) of Eq. (13) and the driven transmon model of Eq. (2) (black crosses). (b) Pulled resonator frequency as obtained from two-tone spectroscopy, with colors representing purities of Floquet modes corresponding to the resonator vacuum (only purities above 0.85 are retained, thus excluding cavity frequency pulls due to chaotic states). Pulls from perturbation theory Eq. (15) are shown as black crosses. (c) Resonator frequency response in the steady state. Only states where the steady-state population times the square of the transition matrix element exceeded 0.001$\kappa$.</figcaption></figure>
regular transmon states of Fig. 8(a), and the corresponding curves terminate for drive strengths where the corresponding transmon level strongly hybridizes with the chaotic layer. The pull is vanishing for the charge like states of the transmon, as can be identified by an inspection of the transmon excitation numbers $\langle\langle N_t \rangle\rangle$ (not shown). On the other hand, the resonator frequency pull associated with chaotic states is not shown here because the two systems strongly hybridize when the transmon is in a chaotic state, rendering impossible the identification of the resonator vacuum. Nonetheless, we address this point below when computing the steady-state density matrix.

Although it fails in the chaotic layer, perturbation theory in the transmon-cavity coupling can be used to compute the cavity pull in the regular regime, see the black symbols in Fig. 8(b). To second order, using Eq. (14), the cavity pull $\chi_i$ associated to the transmon Floquet state $i$ is given by [57]

$$ \chi_i = \sum_{j,k} g^2 n_{ijk}^2 \left( \frac{1}{\omega_a + \Delta_{ijk}} - \frac{1}{\omega_a - \Delta_{ijk}} \right). \quad (15) $$

Here, $n_{ijk}$ and $\Delta_{ijk}$ can be numerically computed by diagonalizing the time-dependent Hamiltonian of the single driven transmon [see Eq. (14) above], thereby fully accounting for the effect of the drive on the transmon. In this way, we avoid the heavy computational cost of the Floquet operator on the full transmon-cavity system. Perturbation theory (black symbols) agrees well with the full simulation only for regular states, which can be identified in Fig. 8(b) by their large purity. In particular, the perturbation theory for the cavity pull of the two computational states of the transmon agree with the full numerical calculations over the whole range of drive amplitudes for which the states remain regular. This suggests using the Floquet basis as a starting point for perturbative treatments of multimode coupling used in two-qubit gate analysis [14,58] or qubit readout [22]. The idea was recently applied to describe the coupling between two resonators induced by a driven ancilla in Ref. [59]. Note that for a small enough coupling strength $g$, the formula, Eq. (15), is expected to be valid also for chaotic transmon states.

Additionally, we show in Fig. 8(c) the frequency response of the cavity in the steady state. This is distinct from the result in Fig. 8(b), which illustrates the cavity pull by each transmon state individually. Specifically, for each value of $\tilde{\epsilon}_d$, we plot points with $y$ coordinate the transition frequency $\Delta_{\alpha\beta k}$, for all possible transitions. The point size is proportional to $\sqrt{p_\alpha|n_{\alpha\beta k}|}$, where $p_\alpha$ is the population of the Floquet mode $\alpha$ in the steady state. That is, Fig. 8(c) weighs the result Fig. 8(b) by steady-state population. A cutoff is imposed on small, indiscernible, weights to reduce plot size. For $\tilde{\epsilon}_d \ge 0.95$ the average cavity response is centered on the bare cavity frequency, with a large variance. In Refs. [13,26], the cavity response at the bare frequency was attributed to the Josephson potential essentially becoming

transparent under strong drives, and the transmon being in a mixture of charge like states. We rather observe here that this response is caused by the transmon entering a low-purity state with significant weight on a large number of chaotic Floquet modes, which corresponds to the transmon becoming strongly hybridized with the cavity and giving vanishing pull on average.

To summarize the results of this section, we have shown that in a driven-circuit QED system, chaos develops along the directions of increasing drive amplitude and increasing resonator photons. For typical circuit QED parameters, the main features of the full-circuit QED spectrum can be explained using the classically driven transmon model, which validates the mean-field analysis of the previous sections. In particular, the onsets of state ionization remain the same. This allows us to interpret the recently observed cavity response at the bare frequency under strong drive [13,24,26] as resulting from chaotic behavior of the driven transmon.

## VI. NON-QND EFFECTS ORIGINATING FROM THE CHAOTIC STATES

The purity drops and the connecting arrows in Fig. 7 can be used to quantify the quantum demolition character of the light-matter interaction in Eq. (13). For example, through its coupling with the resonator, the transmon can relax through the resonator bath. This effect, known as the Purcell decay, results in the left-pointing arrows connecting the low-energy states of the transmon in Fig. 7(a). Moreover, approximating the transmon as a weakly nonlinear oscillator, recent works have shown that the small anharmonicity of the transmon can result in non-QNDness via correlated photon emission of the transmon-resonator systems [21,22,60]. However, the magnitude of these corrections, controlled by the relative transmon anharmonicity and the relative detuning of the transmon with the resonator, remains small compared to the Purcell effect at sufficiently large detuning as considered in this paper. In addition, these correlated effects occur primarily at the qubit frequency, which could be efficiently suppressed by a Purcell filter. Here, we show that the presence of chaotic states in the transmon spectrum leads to non-QND effects that cannot be captured by a perturbative treatment.

In the previous sections, we have shown that part of the transmon spectrum is rendered chaotic when the system is driven and/or coupled to a resonator. Through the coupling with the resonator, see Eq. (14), regular states can decay to chaotic states through multiphoton transitions. This is further illustrated by the large nonlocal arrows in Fig. 7 connecting distant dots in the $(\langle\langle N_t \rangle\rangle, \langle\langle N_r \rangle\rangle)$ plane. The inherently strong nonlinearity of these states cannot be accounted for by perturbative models taking the relative anharmonicity as a small parameter. Moreover, as shown in Sec. IV, chaotic states strongly depend on the
gate charge $n_g$, while this parameter can be gauged away in oscillator models for the transmon where the phase becomes noncompact. We stress that, as compared to perturbative treatments where correlated effects occur at a few dominant frequencies, the transition frequencies involved in this decay are numerous and highly depend on the system parameters.

Spurious bath-assisted decay from the regular to the chaotic states can therefore result in a source of non-QNDness during strong-drive operations such as dispersive readout. In this section, in order to characterize such an effect, we introduce a “chaos-assisted” critical photon number $n_{ca}$, based on the onset of large hybridization between low-lying regular states and higher-energy chaotic states. Due to finite matrix elements between these states, transitions to the chaotic layer arise below $n_{ca}$. We propose to exponentially suppress this effect by increasing the ratio $E_J/E_C$.

### A. Chaos-assisted critical photon number

Going back to our analysis of a driven transmon in the absence of quantum fluctuations from the resonator of Sec. II, an “ionization threshold” between regular and chaotic levels occurs at a critical drive amplitude, and can be graphically identified from the plots of the mean energy of the transmon, see Figs. 1(c) and 1(d). For each regular state, the critical drive corresponds to the drive amplitude $\tilde{\epsilon}_{d,ca}$ for which the level gets “absorbed” into the chaotic layer. As can be observed in Figs. 1(c) and 1(d), as the drive amplitude increases, the chaotic layer grows into a conelike shape whose tip corresponds to the separatrix energy, and thus the boundary between chaotic layer and regular states depends on drive amplitude. The closer a state is to the separatrix at zero drive, the smaller the critical drive amplitude necessary to have it completely hybridize with the chaotic layer. This critical drive amplitude then matches the so-called “ionization threshold” of the respective regular state [27].

Having associated the critical drive amplitude to the boundary between regular and chaotic states, we can observe this boundary in Figs. 1(c) and 1(d) for two values of $\hbar_{eff}$. From these figures, we see that the dependence of the mean energy $\langle\langle H/E_J \rangle\rangle$ at the boundary as a function of $\tilde{\epsilon}_d$ is roughly independent of $\hbar_{eff}$, as this is ultimately tied to the energies of orbits at the regular-to-chaotic boundaries in the classical phase space in Fig. 1(b). Remarkably, this implies that the mean energy $\langle\langle H/E_J \rangle\rangle$ at the boundary solely depends on the parameters of the classical driven pendulum, i.e., $\tilde{\epsilon}_d = \epsilon_d/\omega_p$ and $\tilde{\omega}_d = \omega_d/\omega_p$.

As emphasized in Sec. II, what does depend on $\hbar_{eff}$ is the number of regular states below the chaotic layer. This number controls which mean energies $\langle\langle H/E_J \rangle\rangle$ are accessible by the quantum states, and ultimately the critical drive amplitudes at which regular states get absorbed

into the chaotic layer. As discussed in Sec. III B, the size of the regular central island in the Poincaré sections, e.g., in Figs. 3(a)–3(d), depends on the drive frequency and decreases with drive amplitude. When quantizing the system, in the semiclassical limit, we expect the number of states contained in the low-energy regular subspace to be directly proportional to the inverse of the effective Planck constant $\hbar_{eff}^{-1}$ and to the area of the regular island [36,37]. Hence, the critical drive amplitude for, say, the excited state, is defined as the amplitude at which the area of the regular island is equal to $2\hbar_{eff}$, i.e., it contains two states, the ground and the excited. The critical drive amplitude therefore increases with $\hbar_{eff}^{-1}$ for any given state: the larger the ratio $E_J/E_C$, the larger the critical drive amplitude needed to absorb the level into the chaotic layer.

Similarly, for the case of an undriven transmon coupled to a resonator, one can define a critical photon number $n_{ca} = (\epsilon_{d,ca}/2g)^2$, where $\epsilon_{d,ca}$ is the critical drive amplitude of the driven transmon without coupling to a resonator. We justify this definition by noting that the two systems, i.e., the driven transmon and the undriven transmon coupled to a resonator, have similar spectra (see Sec. D). This hints towards the fact that one could replace the driving field by a quantum field without affecting the ionization threshold. For photon numbers $n > n_{ca}$, we expect the joint transmon-resonator states corresponding to the first excited state of the transmon to hybridize strongly with the states in the chaotic layer. The validity of $n_{ca}$ is numerically demonstrated in Sec. G for the parameters of Sec. V.

A few comments are in order regarding the definition of $n_{ca}$. First, the definition of $n_{ca}$ differs from the usual critical photon number $n_{crit}$ of circuit QED defined in Refs. [16,40]. In these works, $n_{crit}$ is used to ensure the validity of the dispersive approximation through a small hybridization of the transmon (or two-level system) with the resonator. In contrast, $n_{ca}$ corresponds to a threshold for chaos to affect the computational states, and can also be seen as the onset of strong nonlinear effects. Moreover, the definition of $n_{ca}$ results from the analysis of the classical analog of the driven transmon. As such, it does not account for the influence of the offset charge, since the latter does not affect the dynamics of the classical system. As illustrated in Sec. G, the offset charge strongly affects the onset of chaos in the high excited transmon states. However, the photon numbers at which the ground and first excited states “ionize” have a weaker dependence on the offset charge, and $n_{ca}$ provides a good estimate.

### B. Exponential suppression of spurious decay to the chaotic layer with $\hbar_{eff}^{-1}$

In this section, we argue that going to a deeper regime of the transmon, i.e., larger ratio $E_J/E_C$, is a favorable regime in the presence of drive. Indeed, another benefit of working in the large $E_J/E_C$ regime for the transmon is
that the magnitude of the spurious interactions between the low-energy states and the chaotic states is exponentially suppressed with $\hbar_{\text{eff}}^{-1} = \sqrt{E_J/8E_C}$. This suppression mechanism is related to that of charge dispersion with $\sqrt{E_J/E_C}$ [1]. Charge dispersion manifests itself through phase-slip events, which are in turn exponentially suppressed with $\sqrt{E_J/E_C}$ in an undriven transmon. In a driven transmon, chaos-assisted phase slips can occur through the interaction of the low-energy states with the chaotic states (see Sec. IV). Since phase slips involve the full nonlinearity of the cosine potential, they result in multiphoton transitions [61]. Nonetheless, as explained below, these chaos-assisted events are also exponentially suppressed with $\sqrt{E_J/8E_C}$.

In a quantum system, the regular island and the chaotic layer are not dynamically separated as is the case in the classical system. Spurious transitions from the computational states to the chaotic states might occur well below the critical drive amplitude. The regular states have a finite overlap with the chaotic layer in phase space, resulting in finite coupling between regular and chaotic states. Note that CAT, studied in Sec. IV, is an example of such possible regular-to-chaos transition. When a bath or a resonator is coupled to the driven transmon through the charge operator, this results in a finite decay from the regular states to the chaotic layer. Being concerned with high-fidelity operations such as qubit readout, one wished to minimize such effects.

First, to avoid these effects, one can maximize the size of the regular central island by driving at a suitable frequency, as discussed in Sec. III B. Second, one can maximize the number of states contained in the regular island. Previous

<figure><img src="image_11.png" /><figcaption>FIG. 9. Unitless coupling $N_{ij,k \in \{-1,0\}}$ of the ground state (a) and first excited state (b) with the Floquet states of index $l \ge j$ over the Brillouin zones $k = -1, 0$ [see Eq. (16)], for three values of $\hbar_{\text{eff}}^{-1} = 2, 3, 4$ corresponding to the ratios $E_J/E_C = 32, 72, 128$. The drive amplitude is set to $\tilde{\epsilon}_d = 0.1$ (triangles and solid lines), and $\tilde{\epsilon}_d = 0$ for comparison (filled circles and dashed lines). The values of the matrix elements are averaged over 50 values of the offset charge. The coupling to the chaotic layer is exponentially suppressed with $\hbar_{\text{eff}}^{-1}$.</figcaption></figure>

works [43,62] have shown that for a state of a given excitation number, the tunneling rate to the chaotic layer is exponentially suppressed with $\hbar_{\text{eff}}^{-1}$. Intuitively, this follows from an exponential suppression, as a function of $\hbar_{\text{eff}}^{-1}$, of the phase-space overlap of a regular state with the chaotic layer decreases exponentially. In other words, the coupling of a regular state to the chaotic layer is mediated by the regular states in between, resulting in an exponential suppression with the number of intermediate states, thus with $\hbar_{\text{eff}}^{-1}$. Note that at large values of $\hbar_{\text{eff}}^{-1}$, other processes such as resonance-assisted tunneling can alter this exponential law [62].

In practice, this translates to an exponential suppression of the matrix elements of the charge operator between regular and chaotic Floquet modes. As mentioned in Sec. V, the matrix elements $n_{ijk}$ with $k = -1, 0$ between Floquet modes determine the ability of the driven transmon to hybridize with the resonator and potentially decay into the resonator bath similarly to a Purcell effect. For a Floquet state $|i\rangle$ of the driven transmon, we define its unitless coupling to the subspace of “high-energy” Floquet states above a threshold index $j > i$ as

$$N_{ij,k \in \{-1,0\}} = \sqrt{\sum_{l \ge j; k = -1,0} |n_{ilk}|^2}, \quad (16)$$

where the matrix elements $n_{ijk}$ are defined in Eq. (10), and where the Floquet states are sorted through their mean energies defined in Eq. (4). By its definition in Eq. (16), the dimensionless coupling $N_{ij,k \in \{-1,0\}}$ is a monotonically decreasing function. Relevant information about the coupling to the chaotic layer is encompassed in how fast this decrease is as a function of the threshold state index $j$.

In Fig. 9(a) [respectively (b)], for three values of $\hbar_{\text{eff}}^{-1}$, and the same parameters as in Fig. 5, we plot the coupling $N_{ij,k \in \{-1,0\}}$ as a function of the index $j$ for the ground state $i=0$ (respectively, the first excited state $i=1$), at zero drive (dotted line) and $\tilde{\epsilon}_d = 0.1$ (solid line). As the matrix elements fluctuate as a function of the offset charge, they are averaged over 50 values of $n_g$. At zero drive, the matrix elements are exponentially suppressed with the state index $j$ for all three values of $\hbar_{\text{eff}}^{-1}$. At finite drive, after a rapid decay, all three curves stabilize on a plateau marked by comparatively slower decrease in both Figs. 9(a) and 9(b). This plateau is due to the near equality of the matrix elements $n_{ilk}$ when $l$ belongs to the chaotic layer. While the width of the plateau should grow with $\hbar_{\text{eff}}^{-1}$, as the chaotic states are more numerous, its value is exponentially reduced with $\hbar_{\text{eff}}^{-1}$. Due to its proximity to the chaotic states, the first excited state is more strongly coupled to the chaotic states than the ground state, which translates to larger values of the dimensionless coupling $N_{ij,k \in \{-1,0\}}$ for the excited state than for the ground state, both within and away from the plateaus.
Note that at this relatively small drive amplitude ($\tilde{\epsilon}_d = 0.1$) and for a light-matter coupling strength $g/2\pi = 250$ MHz, the transmon-resonator interaction term involves transitions with the chaotic states at a coupling strength of approximately 20 MHz in the case $\hbar_{\text{eff}}^{-1} = 3$ ($E_J/E_C = 72$). The associated frequencies $\Delta_{ij,k-1} - \Delta$ strongly depend on the drive amplitude and the offset charge. With the fluctuations of the offset charge and the ac-Stark shift leading to a sweep across the spectrum during readout [23,27], these transitions might result in non-QND effects during transmon readout.

To summarize this section, we have introduced a chaos-assisted critical photon number $n_{ca}$ to characterize the onset of strong non-QND effects. Below this threshold, spurious transitions to the chaos states remain possible, but transition rates are exponentially suppressed with $\hbar_{\text{eff}}^{-1}$.

## VII. CONCLUSION

In summary, we have shown that the presence of classical chaos has repercussions in the driven-dissipative quantum dynamics of transmon qubits. One consequence is that part of the spectrum of the driven transmon becomes strongly correlated, i.e., it exhibits strong level hybridization, which favors chaos-assisted quantum phase slips that significantly enhance charge dispersion and qubit dephasing. This could also affect the coherence times of other Josephson-based qubits, such as the fluxonium, which involve arrays of hundreds of Josephson junctions [11, 12,63]. However, as pointed out in Ref. [64], the induc-tively shunted transmon remains stable even at large drive powers.

With full-circuit QED simulations, we have also shown that states corresponding to the chaotic layer have lower purity, which indicates that perturbation theories such as the dispersive approximation [16] or black-box quantization [31] become inapplicable. In particular, the Kerr nonlinear oscillator model of the transmon does not capture the large enhancement of charge dispersion and the strong correlations within the spectrum. This situation is analogous to recent studies of arrays of transmons [39], which show that a dispersive theory cannot reproduce the spectrum in the chaotic phase. Moreover, we show that the response of the resonator at large drive amplitudes, commonly referred to as “bright-stating” the qubit, corresponds to entering a steady state massively populating chaotic states, which exhibit high hybridization with the resonator, with the resonator pull averaging to zero [13,65].

This important population of the chaotic states also means that, in order for numerical studies to capture chaos-induced effects on the low-energy states of a single driven transmon, it is necessary to use a Hilbert space size, which contains the entire chaotic layer. That is, in some instances it may be necessary to revisit the conventional wisdom that using a few transmon states is sufficient for accurate

simulations of the driven transmon. Moreover, based on the study of instabilities in the classical system, it may be possible to identify favorable frequency placement for which the width of the chaotic layer is minimal, results that can be used to find optimal parameters for operations with strong drives such as dispersive qubit readout. The identification of the chaotic layer as a function of the classical parameters $\omega_d/\omega_p$ and $\epsilon_d/\omega_p$ leads to the definition of a chaos-induced critical photon number in the quantum system, which increases as a function of $E_J/E_C$. In particular, spurious transitions during strong-drive operations are expected to be minimized in the regime of large $E_J/E_C$.

## ACKNOWLEDGMENTS

We thank Michiel Burgelman, Mark Dykman, Mazyar Mirrahimi, and Alain Sarlette for fruitful discussions. This work was supported by NSERC, the Canada First Research Excellence Fund and the U.S. Army Research Office Grant No. W911NF-18-1-0411. 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.

## APPENDIX A: SUPPRESSION OF CHAOTIC BEHAVIOR AT HIGH-FREQUENCY DRIVING

In Sec. III B, we study the dependence of the chaotic layer on the drive frequency. To have a fair comparison between the effects at different frequencies, we fix the absolute value of the ac-Stark shift of the 0-1 transition to 100 MHz, and correspondingly choose the drive amplitude. This is made possible by tracking the ground and first excited states as a function of drive amplitude to obtain the difference of the quasienergies, $\epsilon_1 - \epsilon_0$. The negative absolute value of the ac-Stark shift is plotted in Fig. 10 as

<figure><img src="image_12.png" /><figcaption>FIG. 10. ac-Stark shift of the 0-1 transition as a function of drive amplitude $\tilde{\epsilon}_d$, for various frequencies, and at $\hbar_{\text{eff}}^{-1} = 2.45$. The dashed line indicates an ac-Stark shift of -100 MHz and the corresponding amplitudes at which the Poincaré section in Fig. 3 are plotted. Tracking is lost beyond a certain drive amplitude, which is well above that needed to work at a fixed shift of -100 MHz.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 11. Rate matrices and steady-state populations for the same parameters as in Fig. 3, with the drive amplitude corresponding to an ac-Stark shift of -100 MHz. The rate matrices (b),(d) are more regular, resulting in a smaller plateau in the steady-state distribution.</figcaption></figure>

a function of $\tilde{\epsilon}_d$ for various frequencies. Disruptions of the lines indicate that tracking is no longer possible.

As a complementary study to Sec. III B, in Figs. 11(a)–11(d), we show the rate matrices at an ac-Stark shift of 100 MHz for the same drive frequencies. The size of the chaotic block follows that of the chaotic layer in the Poincaré sections of Figs. 3(a)–3(d). In particular, the rate matrix corresponding to Fig. 3(d) remains close to that of the undriven system [see inset of Fig. 4(b)], although the 0-1 transition frequency is shifted by 100 MHz.

## APPENDIX B: INVERSION SYMMETRY SECTOR OF THE DRIVEN SYSTEM

In this section, we provide further explanation on the inversion symmetry discussed in Sec. IV. First, let us consider the eigenvalue equation $\mathbf{H}|\psi\rangle = E|\psi\rangle$ with $\mathbf{H}$ the Hamiltonian Eq. (2) of the undriven transmon ($\epsilon_d = 0$). The eigenstates of $\mathbf{H}$ are eigenstates of the parity operator $\mathbf{P}$, defined as $\mathbf{P}: \mathbf{n} \rightarrow -\mathbf{n}, \phi \rightarrow -\phi$, if and only if $n_g \equiv 0, 0.5 \mod 1$. To prove this, we define the boost operator along charge $\mathbf{U} = \exp(-i n_g \phi)$. After transformation under this boost, the Hamiltonian $\mathbf{H}' = \mathbf{U}\mathbf{H}\mathbf{U}^\dagger$ does not depend on $n_g$ and is therefore symmetric under $\mathbf{P}$. Thus parity acting onto an eigenstate must yield the same eigenstate if the spectrum is nondegenerate. That is, $\mathbf{P}|\psi'\rangle = e^{i\theta}|\psi'\rangle$, or $\langle -\phi|\psi'\rangle = e^{i\theta}\langle\phi|\psi'\rangle$, for all $\phi$, and some phase $\theta$ independent of $\phi$. Additionally, in the boosted frame, the eigenstates obey the “twisted-periodic” boundary condition $\langle\phi + 2\pi|\psi'\rangle = e^{-i2\pi n_g}\langle\phi|\psi'\rangle$. Writing the two conditions above at $\phi = \pi$ gives $\langle -\pi|\psi'\rangle = e^{i\theta}\langle\pi|\psi'\rangle = e^{-i2\pi n_g}\langle\pi|\psi'\rangle$, which implies that $e^{i\theta} = e^{i2\pi n_g}$. The wave functions are eigenstates of parity if and only if $\theta = 0, \pi$, that is for $n_g = 0.0, 0.5 \mod 1$.

As explained in the main text, in the transmon regime, the low-energy sector is almost independent of the offset charge and the Hamiltonian becomes *effectively* parity symmetric. When parity is a symmetry of $\mathbf{H}$, matrix elements of the charge operator $\langle i|\mathbf{n}|j\rangle \neq 0$ if and only if $i$ and $j$ belong to different parity sectors. Therefore if $\langle i|\mathbf{n}|i\rangle \neq 0$ then parity symmetry is broken. By the Hellmann-Feynman theorem, $\langle i|\mathbf{n}|i\rangle \propto dE_i/dn_g$, so whenever a band has nonzero group velocity, parity symmetry is broken. This is why sweet spots (zero group velocity) occur at $n_g = 0, 0.5 \mod 1$. Exponentially small group velocity, as in the transmon limit $E_J \gg E_C$, translates to exponentially weak breaking of parity. Thus, in the inset of Fig. 4(b) computed for $n_g = 0.25$ and therefore nonzero group velocity, purple squares appear, but they are exponentially dim. By extension, in the low-energy manifold, transitions that are forbidden by parity selection rules at $n_g = 0, 0.5 \mod 1$ are exponentially suppressed when parity symmetry is broken, leading to an *effective* parity symmetry. Note that it is this *effective* symmetry in the transmon limit $E_J \gg E_C$ that allows us to ignore charge-dispersion effects and express the transmon $\mathbf{H}$ as a Kerr nonlinear oscillator.

At zero drive, this *effective* symmetry results in a suppression of the matrix elements of the charge operator $n_{i,i+2}$, and forbids the transition $i \rightarrow i+2$ in the low-energy sectors for both $n_g = 0.5$ and $n_g = 0.25$. This can be seen in the rate matrices of the undriven systems—insets of Figs. 4(a) and 4(b)—where the squares corresponding to these transitions remain empty. Importantly, this symmetry is not respected at $n_g = 0.25$ for states in the separatrix region where charge dispersion become non-negligible.

In the presence of a drive $\epsilon_d \cos(\omega_d t)\mathbf{n}$, the inversion symmetry holds if one additionally applies the operation
$t \rightarrow t + \pi/\omega_d$ [34]. This symmetry of the driven Hamiltonian defines even and odd parity sectors among the time-dependent Floquet states. Note that the parity of a Floquet state depends on the Brillouin zone it belongs to in the undriven case. For example, in the undriven system, the eigenenergy of the first excited state decomposes as $E_1 = \epsilon_1 + k\omega_d$, where the quasienergy satisfies $|\epsilon_1| \le \omega_d/2$, and $k=1$ for the parameters used in Fig. 4. Within the Floquet formalism, the Floquet state, corresponding to the eigenstate $\exp(-iE_1t)|1\rangle$ of the time-evolution operator, reads $|\psi_1(t)\rangle = \exp(-i\epsilon_1 t)|\phi_1(t)\rangle$, with the Floquet mode $|\phi_1(t)\rangle = \exp(-i\omega_d t)|1\rangle$. Since the Floquet mode $|\phi_1(t)\rangle$ is invariant under the transformation $t \rightarrow t + \pi/\omega_d$ and $n \rightarrow -n$, it belongs to the even sector. When increasing the drive, the state $|\phi_1(t)\rangle$ will remain in the even sector, and expands only on even states $e^{-i(2k+1)\omega_d t}|2n+1\rangle$, $e^{-i2k\omega_d t}|2n\rangle$, with $k \in \mathbb{Z}$ and $n \in \mathbb{N}$.

Under this symmetry, transitions through the charge operator between two states of the same (opposite) parity can only involve an odd (even) number of drive photons $k$, as the charge operator already changes the parity. This results in superposition of blue and red squares in at $n_g = 0.5$, see Fig. 4(a). In contrast, at $n_g = 0.25$, purple squares, indicative of red and blue contributions, light up not only in the chaotic block, but also for transitions involving the low-energy states.

Transitions that do not respect the apparent symmetry of the transmon were experimentally observed in Ref. [23] but remained unexplained. There, the situation is slightly more complicated as a harmonic mode is capacitively coupled to the transmon, making this effectively a three-mode problem. Nonetheless, the situation is similar, and the same conclusions can be drawn. Indeed, the inversion symmetry for the full Hamiltonian holds upon adding the transformation $n_r \rightarrow -n_r$, where $n_r = -i(a - a^\dagger)/\sqrt{2}$ is the charge operator of the resonator. Because the resonator is low $Q$, transitions between states are most likely to be induced by the bath coupled to the resonator through the resonator charge operator $n_r$. Since $n_r$ maps one parity sector to another similarly to $n$, the single driven transmon analysis remains valid for the full driven-circuit QED setup.

## APPENDIX C: FLOQUET SIMULATIONS FOR CIRCUIT QED HAMILTONIAN

In this Appendix, we provide more details for the numerical simulations presented in Sec. V. To capture the entirety of the chaotic layer even at strong drives $\tilde{\epsilon}_d \ge 1$, we use a local Hilbert-space size of 35 for the transmon. For the resonator, we use 20 states, which is pertinent in the regime of off-resonant drives (here, the drive frequency is 500 MHz below the bare cavity frequency at 8 GHz). We characterize truncation errors by plotting the error in the bosonic commutation relation in Fig. 12, and observe significant errors of the commutator only for states with $\langle\langle N_r \rangle\rangle \ge 12$.

<figure><img src="image_14.png" /><figcaption>FIG. 12. Commutator error, encoded in symbol color, for each Floquet mode, expressed as $|1 - \langle\langle [\mathbf{a}, \mathbf{a}^\dagger] \rangle\rangle|$.</figcaption></figure>

We also observe a “reflection” at the Hilbert-space boundary [66], as indicated by monotonically decreasing values of $\langle\langle N_r \rangle\rangle$ versus $\langle\langle N_t \rangle\rangle$ (upper regions of the four panels of Fig. 12).

Boundary effects in the Hilbert space become important at strong drives, where states with high $\langle\langle N_r \rangle\rangle$ become populated, see Fig. 13. In particular, averaging the commutator error in the steady state of the driven-dissipative evolution, we find that it correlates well with the number of occupied Floquet modes as determined from the entanglement entropy of the steady-state density matrix [64] [Figs. 13(b) and 13(c)], which occurs at about the same threshold power $\tilde{\epsilon}$ as ionization, defined as a significant increase in $\langle\langle N_t \rangle\rangle$, see Fig. 13(a).

When computing the steady state, a minimal bound for the required truncation of the resonator Hilbert space can be estimated based on energy considerations. As discussed in Sec. V, one can first diagonalize the driven transmon Hamiltonian, and then write the transmon-resonator coupling in the new basis (Floquet basis for the transmon and Fock basis for the resonator). At strong drives, we need to account for the full chaotic layer of the driven transmon and its corresponding energy span. Although the mean energies of these states concentrate around $\langle\langle H \rangle\rangle \sim 2E_J$, their Fourier decompositions extend over multiple Brillouin zones of width $\omega_d$. More precisely, the highest populated Brillouin zone is approximately given by the energy of the charge state that lies on the external boundary of the chaotic layer. For instance, on Fig. 1(c) the border of the regular and chaotic domain is located around
<figure><img src="image_15.png" /><figcaption>FIG. 13. (a) Steady-state population of the transmon and resonator as a function of drive strength. Both expectation values become nonmonotonic beyond a threshold at $\tilde{\epsilon}_d \simeq 1$. (b) Commutator error (see also Fig. 12) in the steady state. (c) Number of occupied Floquet modes in the steady state [64] as calculated from the entanglement entropy $N_{\text{occ}} = \exp(-\sum_i p_i \log p_i)$.</figcaption></figure>

$\langle\langle H \rangle\rangle_{\text{max}} = 3.5E_J$ at $\tilde{\epsilon}_d = 0.5$, and around $\langle\langle H \rangle\rangle_{\text{max}} = 5.5E_J$ at $\tilde{\epsilon}_d = 1$.

Through the transmon-resonator coupling, the energy stored in the Floquet modes can be converted into resonator photons. In fact, the chaotic behavior of the system with increasing Fock-state number (strong hybridization) shows the ability for the resonator to absorb energy. The maximum energy $\langle\langle H \rangle\rangle_{\text{max}}$ at $\tilde{\epsilon}_d = 0.5$ (respectively, $\tilde{\epsilon}_d = 1$), can be converted into $\langle\langle H \rangle\rangle_{\text{max}}/\omega_a \approx 7$ (respectively, 11) resonator photons. Interestingly, this is in agreement with the range of $\langle\langle N_r \rangle\rangle$ of populated states in Fig. 7 (highlighted with circles). Therefore, the dimension of the resonator Hilbert space has to be chosen larger than the estimated highest populated Fock state.

## APPENDIX D: CHAOS IN THE UNDRIVEN CQED SYSTEM

We consider an undriven cQED system composed of a transmon coupled to a resonator with Hamiltonian

$$ \mathbf{H} = 4E_C(\mathbf{n} - n_g)^2 - E_J \cos(\phi) + \omega_a \mathbf{a}^\dagger \mathbf{a} - i \mathbf{g} \mathbf{n}(\mathbf{a} - \mathbf{a}^\dagger). \quad (\text{D1}) $$

This system is analog to a single driven transmon where the drive is replaced by a quantum field. We show that correlations develop in the spectrum with increasing Fock-state number. By comparing the spectrum of this system

<figure><img src="image_16.png" /><figcaption>FIG. 14. Same as Fig. 7, but for higher drive strengths. Steady-state expectation values of resonator photon number and transmon excitation number deviate significantly from zero. The steady state has significant weight over a large number of chaotic states. An arrow is plotted if $\Gamma_{ij} > 10^{-3}\kappa$ for (a),(b), or $10^{-2}\kappa$ for (c),(d).</figcaption></figure>

with the Floquet spectrum of a driven transmon where an effective drive amplitude is varied, we find that the spectra agree at low photon number and for regular states, but show significant deviations for chaotic states.

We bring the corresponding eigenenergies $E_i$ to the first Brillouin zone delimited by the resonator frequency $\omega_a$, by defining $\epsilon_i = E_i[\omega_a] \in [-\omega_a/2, \omega_a/2]$. In Fig. 15 (black dots), for each eigenstate $|i\rangle$ satisfying $\langle N_r \rangle_i < 20$, we show its energy $\epsilon_i$ and its mean photon number $\langle N_r \rangle_i$, in the $(\epsilon, \langle N_r \rangle)$ plane. We note the appearance of seemingly continuous lines as a function of $\langle N_r \rangle$ (note that these are made of unlabeled dots), that can be associated to transmon states [27]. We refer to these lines as *branches*.

For eigenstates with a given mean photon number $\bar{n} = \langle N_r \rangle$, we can compare the spectrum of the undriven system with the quasispectrum of a single transmon driven (no resonator) at an amplitude $\epsilon_d = 2g\sqrt{\bar{n}}$ (red dots in Fig. 15). Note that large anticrossings appear in both spectra, which is the signature of a correlated spectrum. Although the two spectra qualitatively agree at low photon number, significant deviations appear at large photon number.

As recently observed numerically in Ref. [27], in the undriven system, the frequency pull on the resonator goes to zero as $\langle N_r \rangle$ increases. To see this, let us consider the effect of a probe on the system in a given state $|i\rangle$, with the probe frequency being close to that of the resonator, at
<figure><img src="image_17.png" /><figcaption>FIG. 15. Black dots: eigenenergies of the undriven cQED Hamiltonian in Eq. (D1). For each eigenstate, the energy $\epsilon_i$ and its mean photon number $\langle N_r \rangle_i$ are plotted in the $(\epsilon_i, \langle N_r \rangle_i)$ plane. The energies $\epsilon_i$ are the eigenenergies folded back into the first Brillouin zone defined by the resonator frequency. Only states satisfying $\langle N_i \rangle < 20$ are shown. The parameters are the same as in Fig. 7, with $\omega_a/2\pi = 8$ GHz and $g/2\pi = 0.250$ GHz. Red dots: quasienergies of the Floquet Hamiltonian Eq. (2) with $\omega_d/2\pi = 8$ GHz as a function of $\tilde{n} = (\epsilon_d/2g)^2$, for states satisfying $\langle N_i \rangle < 20$. The behavior of the two spectra is qualitatively similar, but significant deviations appear at large photon number.</figcaption></figure>

$\omega_a + \delta$, where $\delta$ is small. A weak enough probe can only cause transitions to states with neighboring mean photon number $\langle N_r \rangle_i \pm 1$. In addition, since the probe frequency is $\omega_a + \delta$, it can only connect $|i\rangle$ with states of energy $\epsilon_i \pm \delta$, as all the energies are plotted modulo the resonator frequency. Hence, the probe causes transitions to neighboring states in the $(\epsilon_i, \langle N_r \rangle_i)$ plane, which correspond to neighboring states on the same branch. The frequency pull is then $\epsilon_j - \epsilon_i$, where $|j\rangle$ is the state on the same branch with $\langle N_r \rangle_j \approx \langle N_r \rangle_i \pm 1$. Hence, the frequency pulls correspond to the slopes of the apparent branches, which decrease with increasing photon number $\langle N_r \rangle$.

To explain this, let us consider an eigenstate $|i\rangle$ on one of the energy branches, at $\langle N_r \rangle_i$. We pick the state $|j\rangle$ that has maximum overlap with $\mathbf{a}^\dagger |i\rangle / |\mathbf{a}^\dagger |i\rangle|$, so that $i \to j$ is the brightest transition upon driving the resonator when the system is in state $|i\rangle$. By using

$$ \langle j | \mathbf{H} \mathbf{a}^\dagger | i \rangle = \langle j | [\mathbf{a}^\dagger \mathbf{H} + [\mathbf{H}, \mathbf{a}^\dagger]] | i \rangle, \quad (\text{D2}) $$

we obtain the general relation

$$ (E_j - E_i - \omega_a) = g \frac{\langle j | \mathbf{n} | i \rangle}{\langle j | \mathbf{a}^\dagger | i \rangle}. \quad (\text{D3}) $$

While $|j| |\mathbf{n}|$ grows sublinearly with the square root of the photon number $\sqrt{\langle N_r \rangle_i}$ (see Sec. E), $|j| |\mathbf{a}^\dagger|\mathbf{i}|$ scales as $\sqrt{\langle N_r \rangle_i}$ by definition of $|j\rangle$. Hence, the frequency difference satisfies

$$ |\epsilon_j - \epsilon_i| = |E_j - E_i - \omega_a| \propto g/\sqrt{\langle N_r \rangle_i} \to 0. \quad (\text{D4}) $$

## APPENDIX E: STANDARD DEVIATION OF THE TRANSMON DIPOLE MOMENT IN THE CHAOTIC PHASE

We consider a single transmon driven at high enough power to render all of the phaselike states chaotic. The Hilbert space is then simply composed of a finite number $N$ of chaotic states and an infinity of regular charge states. We define $U$ the unitary that diagonalizes the Floquet propagator $U_F = \mathcal{T} \exp[-i \int_0^T \mathbf{H}(t) dt]$, where $\mathcal{T}$ is the time-ordering operator. Because the system is chaotic, the matrix elements $U_{ij}$ of the unitary $U$ projected on the chaotic subspace can be modeled as random variables, distributed with variance $1/N$ and zero mean [38]. The charge-matrix element $n_{ij}$ in the Floquet basis (at a given time) are related to the charge-matrix element $\tilde{n}_{ij}$ in the Floquet basis at zero drive via the unitary $U$, through $\mathbf{n} = U\tilde{\mathbf{n}}U^\dagger$. We can therefore write

$$ n_{ij} = \sum_{r,s} U_{ir} \tilde{n}_{rs} U_{js}^* \quad (\text{E1}) $$

The $n_{ij}$'s can thus be modeled as random variables. Taking the average over these random variables of the modulus squared of the above equality, we obtain

$$ \begin{aligned} \mathbb{E}(|n_{ij}|^2) &= \sum_{r,s,t,u} \mathbb{E}(U_{ir} U_{js}^* U_{it}^* U_{ju}) \tilde{n}_{rs} \tilde{n}_{tu}^* \\ &= \sum_{r,s} \mathbb{E}(|U_{ir}|^2 |U_{js}|^2) |\tilde{n}_{rs}|^2 \\ &\approx \sum_{r,s} \mathbb{E}(|U_{ir}|^2) \mathbb{E}(|U_{js}|^2) |\tilde{n}_{rs}|^2 \\ &= \sum_{r,s} |\tilde{n}_{rs}|^2 / N_{\text{ch}}^2. \end{aligned} \quad (\text{E2}) $$

The second line is obtained by observing that only squared matrix elements $|U_{ir}|^2$ survive when taking the average value. In going from the second to third line, we obtain an approximation by neglecting the correlations between $|U_{ir}|^2$ and $|U_{js}|^2$. Finite correlations exist, since the $U_{ir}$ obey orthonormality conditions. However, using the probability distributions derived in Ref. [38] for random matrices, we check numerically that these second-order correlations are small if the number of chaotic states is large enough, i.e., $N_{\text{ch}} > 8$, and vanish in the large $N_{\text{ch}}$ limit. The last line is obtained using the equality $\mathbb{E}(|U_{ir}|^2) = 1/N_{\text{ch}}$.

The term $\sum_{r,s} |\tilde{n}_{rs}|^2$ can be easily estimated. Indeed, one can write $\sum_{r,s} |\tilde{n}_{rs}|^2 = ||\mathbf{P}_{\text{ch}} \mathbf{n} \mathbf{P}_{\text{ch}}||^2$, where $\mathbf{P}_{\text{ch}}$ is the projector on the chaotic subspace. As the regular states are charge states, a basis of the chaotic subspace are the $N_{\text{ch}}$ charge states $|n_{\text{charge}}\rangle$ with $n_{\text{charge}} < N_{\text{ch}}/2$. We can therefore calculate $||\mathbf{P}_{\text{ch}} \mathbf{n} \mathbf{P}_{\text{ch}}||^2$ in the charge basis, giving
$$ \begin{aligned} ||\mathbf{P}_{\text{ch}}\mathbf{n}\mathbf{P}_{\text{ch}}||^2 &= \sum_{n_{\text{charge}}=-N_{\text{ch}}/2}^{N_{\text{ch}}/2} n_{\text{charge}}^2 \\ &= \frac{(N_{\text{ch}}/2)(N_{\text{ch}}/2+1)(N_{\text{ch}}+1)}{3}. \end{aligned} $$

This leads to

$$ \sqrt{\mathbb{E}(|n_{ij}|^2)} \approx \sqrt{\frac{N_{\text{ch}}}{12}}. $$

In other words, the standard deviation of the dipole moment for a given transition increases with the size of the chaotic layer. Loosely speaking, the action of the unitary $U$ is to “randomize” the matrix elements of the charge operator. As a consequence, the regular structure of the matrix elements of the undriven system, e.g., the harmonic oscillatorlike structure $n_{ij} \propto \delta_{j,i+1} + \delta_{j,i-1}$ of the low-energy states, is not preserved.

For a numerical illustration of the above, we define the mean dipole moment for a Floquet state $i$ by

$$ \langle n_i \rangle = \sqrt{\frac{1}{M} \sum_{j=0}^{M} |n_{ij}|^2}, \quad (\text{E3}) $$

where $M$ is an integer larger than the number of chaotic states $N_{\text{ch}}$. The matrix elements $n_{ij}$ are taken between Floquet modes at time $t=0$. We also define the mean dipole moment

$$ \langle n \rangle = \sqrt{\frac{1}{M} \sum_{i=0}^{M} \langle n_i \rangle^2}. $$

Using the same argument as above, we have $\langle n \rangle = ||\mathbf{P}_M\mathbf{n}\mathbf{P}_M|| \approx \sqrt{M/12}$, where $\mathbf{P}_M$ is the projector on the manifold spanned by the charge states $|m\rangle$ such that $|m| \ge M/2$.

In Fig. 16(a), the mean dipole moment per Floquet state $\langle n_i \rangle$ is plotted for the first 25 Floquet states (sorted by mean energy) as a function of the drive amplitude. For clarity, we highlight in red states that are starting close to the separatrix, in blue the low-energy states, and in black the charge-like states. In green is plotted the mean dipole moment $\langle n \rangle$. Its constant value for $\tilde{\epsilon}_d < 2.5$ is a good indication that the first $M$ states do not hybridize with higher charge like states for this range of drive amplitude. Note that $M=25$ gives $\langle n \rangle = \sqrt{M/12} \approx 1.44$, which agrees with the green curve. For $\tilde{\epsilon}_d > 2.5$, we note that $\langle n \rangle$ increases, indicating the increasing hybridization of higher-energy states with the 25 states represented here. This is also an indication that at $\tilde{\epsilon} = 2.5$, the first 25 states are hybridized. As the drive amplitude increases, the mean dipole moments per Floquet states $\langle n_i \rangle$ converge slowly towards $\langle n \rangle$.

However, in the definition Eq. (E3) of the mean dipole, the diagonal matrix elements $n_{ii}$ also contribute to the average, yielding large dipole moments for charge like states.

<figure><img src="image_18.png" /><figcaption>FIG. 16. (a) Mean dipole moment $\langle n_i \rangle$ associated to each Floquet mode $i$, as a function of the drive amplitude, for the same parameters as in Fig. 1(d), except $n_g = 0.12$. The states are sorted through their mean energies. For each Floquet state $i$, the average in Eq. (E3) is performed over the first 25 transmon states $j$. For clarity, the low-energy states ($0 \le i \le 3$) are colored in blue, the states initially close to the separatrix ($4 \le i \le 8$) appear in red, and the above initially charge like states are colored in black. The average dipole moment $\langle n \rangle$ is shown in green. Its constant value from $\tilde{\epsilon}_d = 0$ to $2.5$ indicates that higher-energy states do not yet hybridize with the represented states here. For $\tilde{\epsilon}_d > 2.5$, we see that $\langle n \rangle$ starts to increase. (b) Same as in (a) but using the definition (E4) for the dipole moment $\langle n'_i \rangle$, where the average is performed over states with $j \ne i$. The dipole moment is low for charge like states until they enter the chaotic layer. The average dipole moment $\langle n' \rangle$ increases with the drive amplitude, due to the increasing number of charge like states with initially high dipole moment $\langle n_i \rangle$ hybridizing with chaotic states. The average dipole moment $\langle n \rangle$ is plotted for comparison as a green dashed line. At $\tilde{\epsilon}_d \approx 2.5$, the states are fully hybridized and $\langle n' \rangle$ is close to $\langle n \rangle$. For $\tilde{\epsilon}_d > 2.5$, hybridization with higher-energy states start to decrease $\langle n' \rangle$, where the average is performed only on the first 25 states. Besides, the dipole moment of low-energy states $\langle n'_i \rangle$ (in blue) dramatically increase upon entering the chaotic layer.</figcaption></figure>

As $n_{ii}$ does not result in transitions to other states, we define the more meaningful dipole moment for the state $i$,

$$ \langle n'_{i} \rangle = \sqrt{\frac{1}{M-1} \sum_{j=0, j \ne i}^{M} |n_{ij}|^2}, \quad (\text{E4}) $$

and the mean dipole moment

$$ \langle n' \rangle = \sqrt{\frac{1}{M} \sum_{i=0}^{M} \langle n'_{i} \rangle^2}. $$
These quantities are plotted in Fig. 16(b), with the same color code used in (a). If the values of $\langle n'_i \rangle$ are qualitatively the same as in Fig. 16(a) for the blue and red states, the behavior is dramatically different for the charge-like states, for which $\langle n'_i \rangle$ is small at low drive power. This is a consequence of the fact that the charge operator is mostly diagonal for charge-like states. Crucially, we note that as the drive strength increases, the mean dipole moments $\langle n'_i \rangle$ converge toward $\langle n' \rangle$, which increases (green curve). This is due to the increasing participation of charge-like states in the chaotic layer. Note that at $\tilde{\epsilon}_d = 2.5$, all 25 states are hybridized, and $\langle n' \rangle$ has come close to $\langle n \rangle$ [green dashed line in Fig. 16(b)].

Finally, note that the mean dipole moment of the ground state, corresponding to the lowest blue line in Fig. 16(b), undergoes a fourfold increase upon entering the chaotic layer for $\tilde{\epsilon}_d > 1$.

## APPENDIX F: ONSET OF STRONG HYBRIDIZATION BETWEEN THE TRANSMON AND THE RESONATOR IN THE PRESENCE OF DRIVE

In this section, we study in detail the onset of the hybridization of the transmon and resonator when the (driven) transmon is in a chaotic state. We start with the Hamiltonian of the driven transmon coupled to a resonator, Eq. (13) in the main text, reproduced here in the absence of the coupling to the bath

$$ \mathbf{H}(t) = 4E_C(\mathbf{n} - n_g)^2 - E_J \cos(\phi) + \epsilon_d \cos(\omega_d t)\mathbf{n} + \omega_a \mathbf{a}^\dagger \mathbf{a} - i \mathbf{g} \mathbf{n} (\mathbf{a} - \mathbf{a}^\dagger). $$

Moving to the Floquet basis $\{|\tilde{i}\rangle\}$ for the driven transmon, and in the frame rotating at the drive frequency, this Hamiltonian reduces to

$$ \mathbf{H}(t) = \sum_i \epsilon_i |\tilde{i}\rangle \langle\tilde{i}| + \Delta \mathbf{a}^\dagger \mathbf{a} - ig \sum_{i,j,k} n_{ijk} |\tilde{i}\rangle \langle\tilde{j}| (e^{i(k-1)\omega_d t} \mathbf{a} - e^{i(k+1)\omega_d t} \mathbf{a}^\dagger), \quad (F1) $$

where $\Delta = \omega_a - \omega_d$, and $\epsilon_i$ are the quasienergies of the driven transmon. Assuming that $|\Delta| \ll \omega_d$ and as $|\epsilon_i| < \omega_d/2$, the interaction term (second line) can come close to resonance for $k = 0, 1$ (term in $\mathbf{a}$) and $k = -1, 0$ (term in $\mathbf{a}^\dagger$). In order to study the impact of chaotic transmon states on the resonator, for simplicity we assume that the drive strength is large enough for all the phase-like states to be ionized, as in Sec. E.

We consider the tensor product state $|\tilde{i}\rangle \otimes |n\rangle$, composed of a chaotic transmon Floquet mode $|\tilde{i}\rangle$ and the nth Fock state of the resonator. The term $-i \tilde{g} n_{ij, -1} a^\dagger$ couples this state to the state $|\tilde{j}\rangle \otimes |n+1\rangle$ with a coupling

<figure><img src="image_19.png" /><figcaption>FIG. 17. Spectrum of the uncoupled driven system corresponding to the Hamiltonian in Eq. (F1) with $g = 0$, $\Delta/2\pi = 200$ MHz, and the same parameters as in Fig. 16. The energies $\epsilon_i$ represent states of the form $|\tilde{i}\rangle \otimes |n\rangle$, where only the first 13 transmon states $|i\rangle$ (sorted through mean energies) are represented here, along with the $n = 0, 1$ Fock states of the resonator. The states $|\tilde{i}\rangle \otimes |0\rangle$ (black dots) and $|\tilde{i}\rangle \otimes |1\rangle$ (blue dots) are translated by $\Delta$. The dots form lines, which are interrupted when hybridization with higher-energy levels occurs (accompanied by a sudden increase of mean energy). As explained in the text, at finite $g$, the two subspaces are coupled through the term $-i g n_{ij, -1} a^\dagger + \text{h.c.}$ At $\tilde{\epsilon}_d = 0.05$, the green line represents the dominant coupling the state $|\tilde{i}\rangle \otimes |0\rangle$ (black dot) with the state $|\tilde{0}\rangle \otimes |1\rangle$ (blue dots). The thickness of line is proportional to the coupling strength. The couplings to other states $|\tilde{j}\rangle \otimes |1\rangle$ are too small to be drawn (1000 times smaller). A dominant coupling results in an energy shift whose sign depends on the detuning. At $\tilde{\epsilon}_d = 1.15$, almost all states are strongly hybridized. As an example, we randomly pick one state $|\tilde{i}\rangle \otimes |0\rangle$ and we plot in red the lines representing the coupling to states of the one-photon subspace (with a slight horizontal offset for clarity). The couplings are numerous and of similar amplitude as shown by the thickness of the lines, with positive and negative detunings. All of these transitions contribute positively and negatively to shifting the energy of the state $|\tilde{i}\rangle \otimes |0\rangle$.</figcaption></figure>

strength $g n_{ij, -1} \sqrt{n}$. Let us first give an intuitive picture of the physics by projecting the Hamiltonian in Eq. (F1) on the vacuum and first few Fock states of the resonator. The corresponding situation is depicted in Fig. 17, where we plot the quasispectra of the uncoupled system corresponding to the Hamiltonian in Eq. (F1) with $g = 0$, for the states $|\tilde{i}\rangle \otimes |0\rangle$ (black dots) and $|\tilde{i}\rangle \otimes |1\rangle$ (blue dots). Only the first 13 states $|\tilde{i}\rangle$ (sorted through mean energies) are represented, which are all part of the chaotic layer of the driven transmon at $\epsilon_d \ge 1$. The two quasispectra are detuned by $\Delta$, due to the term $\Delta \mathbf{a}^\dagger \mathbf{a}$ in Eq. (F1). Here, we choose $\Delta/2\pi = 200$ MHz for visual clarity.

At $\epsilon_d = 0.05$, the possible transitions, due to the term $-i g n_{ij, -1} a^\dagger + \text{h.c.}$, involving the state $|\tilde{i}\rangle \otimes |0\rangle$ and the states $|\tilde{j}\rangle \otimes |1\rangle$, are shown as green vertical lines, where the thickness of the lines is proportional to the coupling strength. These lines couple one black dot and several blue
dots. As expected from the harmonic oscillatorlike structure of the matrix elements, one transition is dominant, corresponding to the transition with the state $|0\rangle \otimes |1\rangle$. The coupling with the latter is the main contribution to the energy shift of the state $|1\rangle \otimes |0\rangle$. Other green lines, representing coupling to other states, do not appear here as their coupling strength is more than a thousand times smaller than that of the dominant transition.

At $\epsilon_d = 1.15$, almost all states represented are strongly hybridized. We pick one state of the form $|\tilde{i}\rangle \otimes |0\rangle$ (represented as a black dot), and similarly, the possible transitions and coupling strength are represented by red lines. To be able to visually distinguish between the numerous lines, we slightly offset these lines from each other on the x axis. However, they represent the couplings of one black dot (corresponding to $|\tilde{i}\rangle \otimes |0\rangle$) to multiple blue dots at the same drive amplitude. Contrary to the situation at low drive power, the transitions are numerous and of similar coupling strength, along with positive and negative detuning. The resulting energy shift on $|\tilde{i}\rangle \otimes |0\rangle$ is null on average (in the sense of random matrix model, see Sec. E).

Strong hybridization between the resonator and the driven transmon can occur through the terms $-\text{ign}_{ij,k=-1,0} a^\dagger + \text{h.c.}$, depending on the ratio between the coupling strength $g_{\text{eff}} = 2gn_{ij,k=-1,0} \sqrt{n+1}$ and the effective detuning $\delta_{\text{eff}}$ of the transition $i \leftrightarrow j$ [40]. Below, we propose a possible estimate of the minimal Fock-state number $n$ for which this happens [40]. Recall from Sec. E that for chaotic states, we have $\sqrt{\mathbb{E}(|n_{ij,k=0}|^2)} \approx \sqrt{N_{\text{ch}}/12}$, where $N_{\text{ch}}$ is the number of chaotic states (typically 10 at the ionization point of the ground state for the parameters in Sec. V). This yields an effective coupling strength $g_{\text{eff}} = 2g\sqrt{n+1}\sqrt{N_{\text{ch}}/12}$. Note that this is an upper bound, as the contributions from the matrix elements $n_{ij,k=-1}$ are usually smaller. Besides, the mean level spacing of the driven transmon is $\omega_d/N_{\text{ch}}$. Accounting for the detuning $\Delta$ of the resonator with respect to the drive, the effective detuning between the set of states $|\tilde{i}_1\rangle \otimes |n\rangle$ and the states $|\tilde{i}_2\rangle \otimes |n+1\rangle$, can be as low as $\delta_{\text{eff}} = \min(\omega_d/N_{\text{ch}}, |\omega_d/N_{\text{ch}} - |\Delta||)$. Note that due to the repulsive statistics, the energy levels do not bunch, and the resulting variance of the effective detuning is rather small. On average, one has $\delta_{\text{eff}} = \omega_d/N_{\text{ch}}$.

The condition $g_{\text{eff}} \sim \delta_{\text{eff}}$ reads $2g\sqrt{n+1}\sqrt{N_{\text{ch}}/12} \sim \omega_d/N_{\text{ch}}$, leading to the definition of a critical photon number,

$$n_{\text{crit}}^{\text{ch}} \sim \frac{3\omega_d^2}{g^2 N_{\text{ch}}^3},$$

above which strong transmon-resonator hybridization occurs. Note that this relation is valid for chaotic transmon states. This criterion is sensitive to the fluctuations of the energy levels and of the charge-matrix elements as a function of the system parameters, and therefore can only give a rough estimate. The strong hybridization between

<figure><img src="image_20.png" /><figcaption>FIG. 18. Same as Fig. 7, for $g/2\pi$ reduced by a factor $\sqrt{6}$.</figcaption></figure>

the systems occurs in the steady state if this condition is satisfied for the lowest Fock states as well.

As an example, for the parameters of Sec. V, at a drive amplitude such that $N_{\text{ch}} \sim 12$ ($N_{\text{ch}}$ can be estimated by how many grid points of the first row come together in Fig. 7), we find $g_{\text{eff}}/2\pi \sim 0.5$ GHz and $\delta_{\text{eff}}/2\pi \sim 0.6$ GHz, leading to $n_{\text{crit}}^{\text{ch}} \sim 0.5$. A large drop in purity is expected to occur even for the lowest Fock states [see Fig. 7(d)]. As a second example, we consider the parameters of Sec. V system with $g = 0.250/\sqrt{6}$ GHz instead of $g = 0.250$ GHz in Fig. 18. This reduction of the coupling $g$ makes the value of $n_{\text{crit}}^{\text{ch}}$ an order of magnitude higher.

## APPENDIX G: CHAOS-ASSISTED CRITICAL PHOTON NUMBER DERIVED FROM POINCARÉ SECTIONS: AGREEMENT WITH DRIVEN QUANTUM SYSTEM

In this section, we numerically check that the critical photon number $n_{\text{ca}}$ introduced in Sec. VI provides a good estimate for the onset of strong nonlinear effects, for the set of parameters used in Sec. V. We determine $n_{\text{ca}}$ as follows. First, from the Poincaré sections of the classical driven pendulum, we extract the area $A$ of the central regular island for various drive amplitudes. Using the Bohr-Sommerfeld quantization rule, $A = 2\pi\hbar_{\text{eff}}(N + 1/2)$, we convert the area to a number of quantum states $N$, here taking $\hbar_{\text{eff}}^{-1} = 3$ ($E_J/E_C = 72$). Converting the drive amplitude to an equivalent photon number $\tilde{n} = (\epsilon_d/2g)^2$ (with $g/2\pi = 0.25$ GHz as in Sec. V), we obtain the photon numbers $\tilde{n}_N$ at which the central island contains
<figure><img src="image_21.png" /><figcaption>FIG. 19. The mean energy spectra of the driven transmon is plotted for six values of the offset charge, as a function of the photon number $\bar{n} = (\varepsilon_d/2g)^2$. The parameters are the same as in Sec. V. The dashed lines indicate the photon numbers at which the area of the central island of the corresponding Poincaré section can contain N quantum states (see text). The red dashed line indicates $n_{ca}$, the photon number above which strong nonlinear effects start affecting the first excited state of the transmon.</figcaption></figure>

$N$ states. The critical photon number is defined as $n_{ca} = \bar{n}_{N=2}$. In the semiclassical limit, the first excited state is expected to ionize for $\bar{n}$ satisfying $\bar{n}_{N=1} \le \bar{n} \le \bar{n}_{N=2}$, and the ground state for $\bar{n}_{N=0} \le \bar{n} \le \bar{n}_{N=1}$.

In Fig. 19, we plot the mean energy spectra of the driven transmon for six values of the offset charge, as a function of the photon number. The dashed lines indicate the values of the photon numbers $n_N$. As expected, the photon numbers at which the ground and first excited states ionize, fall between $\bar{n}_{N=0}$ and $\bar{n}_{N=1}$, and $\bar{n}_{N=1}$ and $\bar{n}_{N=2}$, respectively.

## APPENDIX H: EFFECT OF THE RESONATOR DISSIPATION ON THE SPECTRUM

In Sec. V, we have not accounted for the effect of the resonator dissipation on the spectrum. In the limit of very large $\kappa$, the two systems decouple and chaos disappears. However, for a given Fock state $|n\rangle$ of the resonator, the effect of dissipation on the system (not on resonance, i.e., $\omega_d \ne \omega_a$) should be small as long as the decay rate of $n\kappa$ is small compared to the effective coupling strength $g\sqrt{n}$. This leads to the condition $\kappa \ll g\sqrt{n}$. For typical experimental parameters, $\kappa/2\pi = 10$ MHz and $g/2\pi = 250$ MHz, dissipation will affect Fock states with $n \gtrsim (g/\kappa)^2 = 625$. As a numerical check, we diagonalize the non-Hermitian undriven Hamiltonian

$$ \mathbf{H} = 4E_C(\mathbf{n} - n_g)^2 - E_J \cos(\phi) + \left(\omega_a - i\frac{\kappa}{2}\right) \mathbf{a}^\dagger\mathbf{a} - ig\mathbf{n}(\mathbf{a} - \mathbf{a}^\dagger). \quad (H1) $$

<figure><img src="image_22.png" /><figcaption>FIG. 20. Eigenenergies of the undriven system described by the Hamiltonian in Eq. (D1), with the same parameters as in Fig. 7(a). For each eigenstate, the energy $\varepsilon_i$ and its mean photon number $\langle N_r \rangle_i$ is plotted in the $(\varepsilon_i, \langle N_r \rangle_i)$ plane. The energies $\varepsilon_i$ are the eigenenergies modulo the resonator frequency. Only states satisfying $\langle N_t \rangle < 20$ are shown. Blue dots correspond to $\kappa/2\pi = 0$, and red dots to $\kappa/2\pi = 10$ MHz (see Sec. H). In both cases, the spectrum exhibits correlation at large $\langle N_r \rangle$. The spectrum is unchanged for this range of Fock states under weak dissipation.</figcaption></figure>

In Fig. 20, the real part of the energy $\varepsilon_i$ and the photon number $\langle N_r \rangle_i$ of each state are plotted in the $(\varepsilon_i, \langle N_r \rangle_i)$ plane for $\kappa = 0$ (blue dots) and $\kappa/2\pi = 10$ MHz (red dots). Qualitatively, the spectrum is the same in both cases. In particular, anticrossings, which are signatures of resonances of the coupled system, do not seem to be attenuated.

[1] 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).

[2] N. Earnest, S. Chakram, Y. Lu, N. Irons, R. K. Naik, N. Leung, L. Ocola, D. A. Czaplewski, B. Baker, Jay Lawrence, Jens Koch, and D. I. Schuster, Realization of a System with Metastable States of a Capacitively Shunted Fluxonium, Phys. Rev. Lett. 120, 150504 (2018).

[3] Yen-Hsiang Lin, Long B. Nguyen, Nicholas Grabon, Jonathan San Miguel, Natalia Pankratova, and Vladimir E. Manucharyan, Demonstration of Protection of a Superconducting Qubit from Energy Decay, Phys. Rev. Lett. 120, 150503 (2018).

[4] András Gyenis, Pranav S. Mundada, Agustin Di Paolo, Thomas M. Hazard, Xinyuan You, David I. Schuster, Jens Koch, Alexandre Blais, and Andrew A. Houck, Experimental Realization of a Protected Superconducting Circuit Derived from the $0-\pi$ Qubit, PRX Quantum 2, 010339 (2021).

[5] J. Q. You, Xuedong Hu, S. Ashhab, and Franco Nori, Low-decoherence flux qubit, Phys. Rev. B 75, 140515 (2007).
[6] Fei Yan, Simon Gustavsson, Archana Kamal, Jeffrey Birenbaum, Adam P Sears, David Hover, Ted J. Gudmundsen, Danna Rosenberg, Gabriel Samach, S Weber, Jonilyn L. Yoder, Terry P. Orlando, John Clarke, Andrew J. Kerman, and William D. Oliver, The flux qubit revisited to enhance coherence and reproducibility, Nat. Commun. 7, 12964 (2016).

[7] Mazyar Mirrahimi, Zaki Leghtas, Victor V. Albert, Steven Touzard, Robert J. Schoelkopf, Liang Jiang, and Michel H. Devoret, Dynamically protected cat-qubits: A new paradigm for universal quantum computation, New J. Phys. 16, 045014 (2014).

[8] Z. Leghtas, S. Touzard, I. M. Pop, A. Kou, B. Vlastakis, A. Petrenko, K. M. Sliwa, A. Narla, S. Shankar, M. J. Hatridge, M. Reagor, L. Frunzio, R. J. Schoelkopf, M. Mirrahimi, and M. H. Devoret, Confining the state of light to a quantum manifold by engineered two-photon loss, Science 347, 853 (2015).

[9] Raphaël Lescanne, Marius Villiers, Théau Peronnin, Alain Sarlette, Matthieu Delbecq, Benjamin Huard, Takis Kontos, Mazyar Mirrahimi, and Zaki Leghtas, Exponential suppression of bit-flips in a qubit encoded in an oscillator, Nat. Phys. 16, 509 (2020).

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

[11] Nicholas A. Masluk, Ioan M. Pop, Archana Kamal, Zlatko K. Minev, and Michel H. Devoret, Microwave Characterization of Josephson Junction Arrays: Implementing a Low Loss Superinductance, Phys. Rev. Lett. 109, 137002 (2012).

[12] Vladimir E. Manucharyan, Nicholas A. Masluk, Archana Kamal, Jens Koch, Leonid I. Glazman, and Michel H. Devoret, Evidence for coherent quantum phase slips across a Josephson junction array, Phys. Rev. B 85, 024521 (2012).

[13] Raphaël Lescanne, Lucas Verney, Quentin Ficheux, Michel H. Devoret, Benjamin Huard, Mazyar Mirrahimi, and Zaki Leghtas, Escape of a Driven Quantum Josephson Circuit into Unconfined States, Phys. Rev. Appl. 11, 014030 (2019).

[14] Sarah Sheldon, Easwar Magesan, Jerry M. Chow, and Jay M. Gambetta, Procedure for systematically tuning up cross-talk in the cross-resonance gate, Phys. Rev. A 93, 060302 (2016).

[15] Hanhee Paik, A. Mezzacapo, Martin Sandberg, D. T. McClure, B. Abdo, A. D. Córcoles, O. Dial, D. F. Bogorin, B. L. T. Plourde, M. Steffen, A. W. Cross, J. M. Gambetta, and Jerry M. Chow, Experimental Demonstration of a Resonator-Induced Phase Gate in a Multiqubit Circuit-QED System, Phys. Rev. Lett. 117, 250502 (2016).

[16] 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).

[17] 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).

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

[19] A. Grimm, N. E. Frattini, S. Puri, S. O. Mundhada, S. Touzard, M. Mirrahimi, S. M. Girvin, S. Shankar, and M. H. Devoret, Stabilization and operation of a Kerr-cat qubit, Nature 584, 205 (2020).

[20] 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).

[21] Moein Malekakhlagh, Alexandru Petrescu, and Hakan E. Türeci, Lifetime renormalization of weakly anharmonic superconducting qubits. I. Role of number nonconserving terms, Phys. Rev. B 101, 134509 (2020).

[22] 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).

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

[24] I. Pietikäinen, S. Danilin, K. S. Kumar, A. Vepsäläinen, D. S. Golubev, J. Tuorila, and G. S. Paraoanu, Observation of the Bloch-Siegert shift in a driven quantum-to-classical transition, Phys. Rev. B 96, 020501 (2017).

[25] I. Pietikäinen, J. Tuorila, D. S. Golubev, and G. S. Paraoanu, Photon blockade and the quantum-to-classical transition in the driven-dissipative Josephson pendulum coupled to a resonator, Phys. Rev. A 99, 063828 (2019).

[26] Lucas Verney, Raphaël Lescanne, Michel H. Devoret, Zaki Leghtas, and Mazyar Mirrahimi, Structural Instability of Driven Josephson Circuits Prevented by an Inductive Shunt, Phys. Rev. Appl. 11, 024003 (2019).

[27] 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).

[28] 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).

[29] Steven Tomsovic and Denis Ullmo, Chaos-assisted tunneling, Phys. Rev. E 50, 145 (1994).

[30] K. A. Matveev, A. I. Larkin, and L. I. Glazman, Persistent Current in Superconducting Nanorings, Phys. Rev. Lett. 89, 096802 (2002).

[31] Simon E. Nigg, Hanhee Paik, Brian Vlastakis, Gerhard Kirchmair, S. Shankar, Luigi Frunzio, M. H. Devoret, R. J. Schoelkopf, and S. M. Girvin, Black-Box Superconducting Circuit Quantization, Phys. Rev. Lett. 108, 240502 (2012).

[32] Heinz-Peter Breuer, Wolfgang Huber, and Francesco Petruccione, Quasistationary distributions of dissipative nonlinear quantum oscillators in strong periodic driving fields, Phys. Rev. E 61, 4883 (2000).
[33] Roland Ketzmerick and Waltraut Wustmann, Statistical mechanics of Floquet systems with regular and chaotic states, Phys. Rev. E 82, 021114 (2010).

[34] Milena Grifoni and Peter Hänggi, Driven quantum tunneling, Phys. Rep. 304, 229 (1998).

[35] Boris V. Chirikov, A universal instability of many-dimensional oscillator systems, Phys. Rep. 52, 263 (1979).

[36] I. C. Percival, Regular and irregular spectra, J. Phys. B: At. Mol. Phys. 6, L229 (1973).

[37] M. V. Berry, Regular and irregular semiclassical wavefunctions, J. Phys. A: Math. Gen. 10, 2083 (1977).

[38] Fritz Haake, *Quantum Signatures of Chaos* (Springer, Berlin, Germany, 2010).

[39] Christoph Berke, Evangelos Varvelis, Simon Trebst, Alexander Altland, and David P. DiVincenzo, Transmon platform for quantum computing challenged by chaotic fluctuations, Nat. Commun. 13, 1 (2022).

[40] Alexandre Blais, Arne L. Grimsmo, S. M. Girvin, and Andreas Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).

[41] Georgin Moiseevich Zaslavski, R. Z. Sagdeev, D. A. Usikov, and A. A. Chernikov, *Weak Chaos and Quasi-Regular Patterns*, edited by A. R. Translator Sagdeeva, Cambridge Nonlinear Science Series (Cambridge University Press, 1991).

[42] Roland Ketzmerick and Waltraut Wustmann, Statistical mechanics of Floquet systems with regular and chaotic states, Phys. Rev. E 82, 021114 (2010).

[43] A. Bäcker, R. Ketzmerick, S. Löck, and L. Schilling, Regular-to-Chaotic Tunneling Rates Using a Fictitious Integrable System, Phys. Rev. Lett. 100, 104101 (2008).

[44] R. Graham, Quantization of Poincaré Maps, Europhys. Lett. (EPL) 7, 671 (1988).

[45] José A. González and Mariano A. del Olmo, Coherent states on the circle, J. Phys. A: Math. Gen. 31, 8841 (1998).

[46] M. L. Mehta, *Random Matrices*, revised and enlarged 2nd ed., edited by M. L. Mehta (Academic Press, San Diego, 1991), p. 1.

[47] Nikolaus Bubner and Robert Graham, Quantum dynamics in a chaotic separatrix layer, Phys. Rev. A 43, 1783 (1991).

[48] Maxime Martinez, Olivier Giraud, Denis Ullmo, Juliette Billy, David Guéry-Odelin, Bertrand Georgeot, and Gabriel Lemarié, Chaos-Assisted Long-Range Tunneling for Quantum Simulation, Phys. Rev. Lett. 126, 174102 (2021).

[49] Jens Koch, V. Manucharyan, M. H. Devoret, and L. I. Glazman, Charging Effects in the Inductively Shunted Josephson Junction, Phys. Rev. Lett. 103, 217004 (2009).

[50] A. Mouchet, C. Miniatura, R. Kaiser, B. Grémaud, and D. Delande, Chaos-assisted tunneling with cold atoms, Phys. Rev. E 64, 016221 (2001).

[51] M. Arnal, G. Chatelain, M. Martinez, N. Dupont, O. Giraud, D. Ullmo, B. Georgeot, G. Lemarié, J. Billy, and D. Guéry-Odelin, Chaos-assisted tunneling resonances in a synthetic Floquet superlattice, Sci. Adv. 6, eabc4886 (2020).

[52] Ziwen Huang, Pranav S. Mundada, András Gyenis, David I. Schuster, Andrew A. Houck, and Jens Koch, Engineering

Dynamical Sweet Spots to Protect Qubits from 1/f Noise, Phys. Rev. Appl. 15, 034065 (2021).

[53] G. Ithier, E. Collin, P. Joyez, P. J. Meeson, D. Vion, D. Esteve, F. Chiarello, A. Shnirman, Y. Makhlin, J. Schriefl, and G. Schön, Decoherence in a superconducting quantum bit circuit, Phys. Rev. B 72, 134519 (2005).

[54] B. G. Christensen, C. D. Wilen, A. Opremcak, J. Nelson, F. Schlenker, C. H. Zimonick, L. Faoro, L. B. Ioffe, Y. J. Rosen, J. L. DuBois, B. L. T. Plourde, and R. McDermott, Anomalous charge noise in superconducting qubits, Phys. Rev. B 100, 140503 (2019).

[55] C. D. Wilen, S. Abdullah, N. A. Kurinsky, C. Stanford, L. Cardani, G. D'Imperio, C. Tomei, L. Faoro, L. B. Ioffe, C. H. Liu, A. Opremcak, B. G. Christensen, J. L. DuBois, and R. McDermott, Correlated charge noise and relaxation errors in superconducting qubits, Nature 594, 369 (2021).

[56] Maxime Boissonneault, J. M. Gambetta, and Alexandre Blais, Improved Superconducting Qubit Readout by Qubit-Induced Nonlinearities, Phys. Rev. Lett. 105, 100504 (2010).

[57] Vera Gramich, Simone Gasparinetti, Paolo Solinas, and Joachim Ankerhold, Lamb-Shift Enhancement and Detection in Strongly Driven Superconducting Circuits, Phys. Rev. Lett. 113, 027001 (2014).

[58] Moein Malekakhlagh, Easwar Magesan, and David C. McKay, First-principles analysis of cross-resonance gate operation, Phys. Rev. A 102, 042605 (2020).

[59] Yaxing Zhang, Brian J. Lester, Yvonne Y. Gao, Liang Jiang, R. J. Schoelkopf, and S. M. Girvin, Engineering bilinear mode coupling in circuit QED: Theory and experiment, Phys. Rev. A 99, 012314 (2019).

[60] Ryo Hanai, Alexander McDonald, and Aashish Clerk, Intrinsic mechanisms for drive-dependent Purcell decay in superconducting quantum circuits, Phys. Rev. Res. 3, 043228 (2021).

[61] M. Houzet and L. I. Glazman, Critical Fluorescence of a Transmon at the Schmid Transition, Phys. Rev. Lett. 125, 267701 (2020).

[62] Michael Sheinman, Shmuel Fishman, Italo Guarneri, and Laura Rebuzzini, Decay of quantum accelerator modes, Phys. Rev. A 73, 052110 (2006).

[63] Long B. Nguyen, Yen-Hsiang Lin, Aaron Somoroff, Raymond Mencia, Nicholas Grabon, and Vladimir E. Manucharyan, High-Coherence Fluxonium Qubit, Phys. Rev. X 9, 041041 (2019).

[64] Michiel Burgelman, Pierre Rouchon, Alain Sarlette, and Mazyar Mirrahimi, Structurally Stable Subharmonic Regime of a Driven Quantum Josephson Circuit, Phys. Rev. Appl. 18, 064044 (2022).

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

[66] H. P. Breuer and M. Holthaus, Adiabatic processes in the ionization of highly excited hydrogen atoms, Z. Phys. D: At. Mol. Clusters 11, 1 (1989).