Landau–Zener without a qubit: multiphoton sidebands interaction and signatures of dissipative quantum chaos

Experimental system and model

We aim to investigate all regimes of nonlinearity and dissipation, namely, qubit, Kerr, Duffing, and linear, as shown in Fig. 1a–d, respectively. To this extent, we design and fabricate two frequency-tunable nonlinear resonators that can operate in these different regimes according to the driving amplitude. These are superconducting SQUID arrays^65,66,67, galvanically connected to the ground on one side, and capacitively shunted to the ground on the other side, as shown in Fig. 2a, b, and f. A detailed summary of their parameters is reported in Table 1.

a Equivalent lumped electrical circuit of the device composed of N SQUIDs in series. The static flux is controlled via an external coil and a microwave signal can be applied to the flux line (purple) for fast frequency modulation. The resonator is coupled to a feedline in a notch configuration. b, f Optical micrograph of the two SQUID array resonators studied in this work with N = 10 (red) and N = 32 (blue) SQUIDs, with details on the Josephson junction shown in the inset where the white scale bar represents 5 μm. c, d For the Kerr case (N = 10), the magnitude and phase of the transmission coefficient through the feedline as a function of drive power at the sample. The black dashed curve indicates where the phase of S₂₁ is zero according to numerical simulations, highlighting the position of the multiphoton resonances. e Selected traces of the data reported in (d). Experimental data are shown with circle markers whose colors correspond to the ticks in (d). Numerical fits to a full quantum model are shown in black solid lines. The gray dotted vertical lines indicate the position of the multiphoton resonances obtained from numerical simulations. g–i Same measurement as in (c–e), but for the Duffing case. The dashed curves in (g, h) indicate the minima of ∣S₂₁∣ obtained from a full quantum simulation.

The frequency of the resonators can be tuned by a dedicated flux line that uniformly threads the magnetic fields in each SQUID loop (in purple) and an external superconducting coil. The two resonators differ in the number of SQUIDs in each array, as highlighted by red and blue false colors, determining the two orders of magnitude difference in their Kerr nonlinearity χ. Increasing the number of SQUIDs in an array leads to a decrease in nonlinearity due to the reduced Josephson inductance required to maintain a fixed frequency and the diminished phase fluctuations across each junction⁶⁸. Hereafter, the red and blue color schemes will always indicate measurements of the devices in the Kerr/qubit and Duffing/linear regimes, respectively. Each resonator is also capacitively coupled to a feedline (in green) in a notch configuration, resulting in an external coupling κ_ext close to the internal dissipation rate κ_int (critically coupled regime). We define the total dissipation rate as κ = κ_ext + κ_int.

Each device is thermally and mechanically anchored at the mixing chamber plate of a dilution refrigerator, reaching an average base temperature of 15 mK. The devices are probed by a coherent drive with amplitude F at the sample, and injected in the feedline through highly attenuated coaxial lines. The drive amplitude is related to the input power P_in by \(F=\sqrt{{P}_{{\rm{in}}}{\kappa }_{{\rm{ext}}}/\hslash {\omega }_{{\rm{d}}}}\), where ω_d is the drive frequency. Although the frequency of the untuned cavity is ω_c, through the paper, the frequency working point of the resonators, ω_wp, is set by a static flux generated by direct current through an external superconducting coil. The flux operating point is Φ_wp = 0.45Φ₀ for the N = 10 device, while Φ_wp = 0.32Φ₀ for the N = 32 one, where Φ₀ = h/2e is the magnetic flux quantum. Finally, driving the fluxline at a frequency Ω periodically modulates the frequency of the resonator, approximately between ω_wp ± ζ, with ζ representing the strength of the modulation. The single-tone spectroscopy of the resonators at low-driving power as a function of the external magnetic flux is reported in the Supplementary Information.

As explained in the Supplementary Information, both devices can be modeled as a bosonic mode with the following time-dependent Hamiltonian:

where Δ = ω_d − ω_wp is the detuning between the working point of the devices (ω_wp) and the drive frequency (ω_d). Beyond the total dissipation rate κ, the system is subject to dephasing with rate κ_ϕ. We include them in the time evolution of the density matrix \(\hat{\rho }\) using the Lindblad master equation

Here, \({\mathcal{D}}[\hat{L}]\hat{\rho }\equiv \hat{L}\hat{\rho }{\hat{L}}^{\dagger }-\{{\hat{L}}^{\dagger }\hat{L},\hat{\rho }\}/2\) is the Lindblad dissipator⁶⁹.

In Fig. 2 we characterize the coherent response of the resonators in the absence of modulation (i.e., ζ = 0) and use it to determine the parameters of the two devices. In Fig. 2c, d we report the magnitude and phase of the transmission coefficient S₂₁ (see “Methods”) as a function of the driving power P_in ∝ F² in the Kerr regime. At low power, only a single dip around Δ = 0 is visible, representing the transition to the first excited state (noted as \(\left\vert 0\right\rangle \to \left\vert 1\right\rangle\), where \(\left\vert n\right\rangle\) is the photon number state of the resonator). At larger values of the drive, several dips appear, representing the so-called Kerr multiphoton transitions between the ground state and the higher-excited levels (\(\left\vert 0\right\rangle \to \left\vert n\right\rangle\)), highlighted in the single traces shown in Fig. 2e. According to Eq. (1), all the dips should appear at Δ ≃ (n − 1)χ. Small deviations from this prediction are due to higher-order nonlinearities (see “Methods”). For even larger powers, the nonlinearity suppresses the intracavity photon number with respect to the input power, resulting in an almost unitary transmission S₂₁.

We report the same measurements for the resonator in the Duffing regime in Fig. 2g–i. In this case, dissipation smears the multiphoton resonances, resulting in an indistinguishable level structure. Increasing the drive, the single dip of S₂₁ originally at Δ = 0 moves to negative detunings, indicating that the drive is exciting higher levels. Scanning the detuning from negative to positive values, as done in Fig. 2i, reveals the presence of a sharp jump, where the resonator passes from a highly- to a lowly-populated phase. This behavior is associated with optical bistability, i.e., the presence of two metastable states that require a long time to decay to the steady state^42,43. This phenomenon gives rise to hysteresis^40,41 and makes it difficult to properly resolve the exact detuning where the transition occurs.

Linear and qubit LZSM interference

We can investigate the linear and qubit regimes using the N = 10 and N = 32 resonators described above. Indeed, for the N = 10 resonator, the second-excited level is not significantly populated if F² ≪ ∣χ∣κ (see Eq. (21) “Methods”). For the N = 10 device parameters and the drive F/2π ≃ 1.6 MHz considered here, the third level is predicted to be populated less than 0.03%. For the values used in this first part of the experiment, the system effectively behaves like an ideal qubit subject to dissipation and dephasing. We report the experimental data in Fig. 3b, c. In Fig. 3b we show the norm of the scattering coefficient S₂₁ sweeping the detuning Δ, for a fixed modulation frequency Ω, and varying the modulation strength ζ. One observes the LZSM pattern emerging, with populated regions at Δ = mΩ, for integer m. Fixing ζ and scanning Ω, in Fig. 3c we observe again the interference pattern at Δ = mΩ. We thus confirm the presence of LZSM interference and the control over the modulation of the resonator in the Kerr regime.

a–c (red) Analysis of the N = 10 device, with input power P_in = −138.8 dBm, ensuring that we are in the qubit regime. b The transmission coefficient ∣S₂₁∣ as a function of the detuning Δ and the modulation strength ζ, for fixed modulation frequency Ω/2π = 30 MHz [see Eq. (1)]. a Comparison of the experimental and theoretical data for ∣S₂₁∣. Solid lines represent the results of the numerical simulations of the full quantum model obtained at Δ = 0, Δ = −Ω, and Δ = −2Ω (see Supplementary Information). The circles are the experimental data obtained from panel (b), in which Δ is slightly re-scaled to account for the nonlinear flux-dependency of the resonator frequency (see Supplementary Information). c ∣S₂₁∣ as a function of Δ and Ω. d–f (blue) As in (a–c), but for the N = 32 device, with P_in = −133.3 dBm to ensure that the system is in the linear regime. From these plots, the two regimes appear almost indistinguishable. g The photon number vs Δ and ζ is obtained from a simulation using the effective model of Eq. (3) that reproduces the interference pattern in (b, e). h–j Depiction of the time evolution of the energy level \(\left\vert 1\right\rangle\), in the frame rotating at the drive frequency ω_d, if F = 0. A finite drive F opens gaps at each crossing between \(\left\vert 0\right\rangle\) and \(\left\vert 1\right\rangle\), allowing a non-adiabatic passage between the two. The parameters Δ and ζ are indicated by green markers in (g). h At Δ = 0, the level \(\left\vert 1\right\rangle\) becomes resonant with \(\left\vert 0\right\rangle\) (they form a level crossing, see the inset). The values of ζ, F, and κ then determine the probability of transitioning out of the vacuum. i For non-zero detuning (e.g., ∣Δ∣ = Ω) and small modulation (ζ ≪ Ω), the level \(\left\vert 1\right\rangle\) is never resonant with \(\left\vert 0\right\rangle\) and it cannot be populated. j For strong enough modulation ζ > ∣Δ∣, the level \(\left\vert 1\right\rangle\) can form again an avoided level crossing, and constructive interference is possible again.

We now consider the N = 32 device. In this case, the oscillator approximately behaves as a purely linear resonator if \(F < \sqrt{{\kappa }^{3}/| \chi | }\) (see Eq. (24) in Methods). For the N = 32 device parameters and the drive F/2π ≃ 3 MHz considered here, we estimate a relative photon-number deviation from a completely linear resonator of less than 3%. Within this regime, we repeat the previous measurements and report them in Fig. 3e, f. Surprisingly, we observe the same interference pattern emerging, with no distinguishable differences between the qubit and the completely linear case. This feature indicates that only the energy difference between \(\left\vert 0\right\rangle\) and \(\left\vert 1\right\rangle\) determines the interference pattern in both the qubit and the linear regimes (c.f. Fig. 3g–j). This similarity may be expected from linear response theory; indeed, weakly driven nonlinear oscillators should behave similarly even if they have widely different anharmonicities. The presence of LZSM interference seems even more general, as it should be observable even in a purely linear cavity and for an arbitrarily large number of photons. We remark that the frequency modulation of almost linear oscillators has been used to perform mode-conversion⁷⁰, parametric amplification⁷¹, and squeezing⁷². However, to the best of our knowledge, LZSM interferences were not studied in linear oscillators, and this discussion is missing from recent reviews on the topic^7,73.

To provide more quantitative reasoning, we choose \(\bar{m}\) minimizing \(\Delta -\bar{m}\Omega\) and, following the procedure derived in the Supplementary Information, and passing in the frame rotating at the frequency \(\bar{m}\Omega\), we have

where the renormalized detuning \({\Delta }_{\bar{m}}\) and renormalized drive \({F}_{\bar{m}}\) are

with \({J}_{\bar{m}}\left(\zeta /\Omega \right)\) indicating the Bessel function of the first kind. All dissipative terms maintain their form as in (2). In other words, when we can single out a single relevant frequency \({\Delta }_{\bar{m}}\) for each of the LZSM interference dips, the devices behave as a collection of independent nonlinear resonators, whose driving amplitudes \({F}_{\bar{m}}\) are modulated via Bessel functions. For the parameters we consider here, and if we also assume a weak enough drive to be in the linear and qubit regime (see Eq. (21), Eq. (24) “Methods”), we obtain

with β = 1 for a linear resonator regime and β = 4 in the weakly driven qubit limit. Namely, the two regimes have identical interference patterns, only slightly modulated by the dephasing rate κ_ϕ. To further demonstrate the validity of these results, additional LZSM interference patterns are reported in the Supplementary Information, highlighting the precise control of the number and frequency spacing of modes over a broad range of modulation strengths and frequencies.

The approximation of the effective model correctly captures the value of the photon number, but not that of the field \(\hat{a}\) (and thus cannot be used to quantitatively study S₂₁). As is discussed in the Supplementary Information, to correctly capture this feature, one has to resort to a full quantum simulation of the Floquet model. This is shown in Fig. 3a, d, where we plot ∣S₂₁∣ of the first three LZSM lobes, comparing the experimental data with the theoretical predictions both for the qubit and linear regimes. In both cases, we find an excellent agreement between theory and experiments. We note that the maxima and the minima of ∣S₂₁∣ of the \(\bar{m}{\rm{th}}\) LZSM mode coincide with the extremes of the associated Bessel function \({J}_{\bar{m}}\), showing the qualitative validity of Eq. (3) in describing also S₂₁.

We remark here that the Hamiltonian in Eq. (3) could be obtained by approximating the response of an array of nonlinear resonators, each at a frequency \({\Delta }_{\bar{m}}\). Therefore, we can interpret each of the LZSM dips as the response of a different Floquet synthetic mode. As we show below, by increasing the drive, these initially non-interacting modes will begin to interact.

LZSM beyond the qubit approximation: Kerr regime

We now focus on those phenomena emerging due to the simultaneous presence of the multilevel structure of nonlinear resonators and the modulation of their eigenenergies, studying the devices beyond their qubit and linear regimes.

In the Kerr regime and for strong enough drives to probe the multiphoton transitions (see Eq. (21)), we investigate how the frequency and amplitude of the modulation modifies the multiphoton resonances.

The system’s behavior around the multiphoton resonance \(\left\vert 0\right\rangle \to \left\vert n\right\rangle\) occurring for Δ ≃ χ(n − 1) can be described by a 2 × 2 matrix. For instance, the \(\left\vert 0\right\rangle \to \left\vert 2\right\rangle\) multiphoton transition can be described as

where G⁽²⁾ represents the effective drive between the vacuum and the state \(\left\vert 2\right\rangle\). For ζ = 0, one has G⁽²⁾ = F²/χ for Δ = χ. This formula can be generalized to obtain G⁽ⁿ⁾ for arbitrary \(\left\vert 0\right\rangle \to \left\vert n\right\rangle\) transitions⁷⁴. The dissipation maintains its form, instead. As we show below, the fundamental parameter to describe these phenomena is the ratio between the modulation frequency Ω, determining the position of the LZSM sidebands, and the nonlinearity χ, determining the position of the multiphoton resonance.

Strong modulation case

We first choose Ω ≫ ∣χ∣ (strongly modulated case). In Fig. 4a–c we report the scattering coefficient ∣S₂₁∣ as a function of the detuning Δ and the strength ζ of the modulation. As the drive amplitude F is increased, several additional dips appear, signaling the transitions between the photon number states \(\left\vert 0\right\rangle\) and \(\left\vert n\right\rangle\) of the resonator. These dips occur at a frequency lower than each main LZSM dip associated with the transition \(\left\vert 0\right\rangle \to \left\vert 1\right\rangle\). Each new additional dip is detuned by the same frequency as the unmodulated multiphoton resonances shown in Fig. 2c–e. Within a first approximation, this effect is due to the interplay between the modulation in Eq. (3) and the nonlinearity of the system, as shown in Fig. 4d reporting the result of a numerical simulation.

a–c The magnitude of S₂₁ is measured vs Δ and ζ for fixed modulation frequency Ω/2π = 150 MHz. As the drive power P_in is increased, Kerr multiphoton resonances from \(\left\vert 0\right\rangle\) to \(\left\vert n\right\rangle\) appear detuned by (n − 1)χ on the left of bare LZSM resonances. For large ζ, notice the shift of the pattern to negative detuning, due to the nonlinear dependence of the SQUID array frequency on the flux, as explained in the Supplementary Information. d Photon-number simulation using the effective model of Eq. (3) for the same parameters as in (b), recovering the same interference pattern. e–g In the drive frame, energy vs time for different values of ζ and Δ including the first three levels of an undriven Kerr resonator (F = 0). Green markers indicate the corresponding value of Δ and ζ in (d). e For Δ = 0, although multiple levels cross with \(\left\vert 0\right\rangle\), only the level \(\left\vert 1\right\rangle\) forms a constructive interference. f For Δ = χ, the second level \(\left\vert 2\right\rangle\) crosses \(\left\vert 0\right\rangle\), and an appropriate choice of parameters leads to constructive interference. g For Δ = χ + Ω, similar LZSM interference can be constructive again and the level \(\left\vert 2\right\rangle\) can be populated. We verified that both the data and full numerical simulations recover that the interference patterns are fully constructive at Δ = Ω and ζ ≈ 1.84Ω, where the Bessel function J₁(ζ/Ω) is at a maximum, confirming the prediction of the effective model in Eq. (3).

To explain this behavior, we can assume that, around each of the LZSM dips, we again have a drive of the same form as Eq. (3). When we then match the condition for a multiphoton resonance, it is this effective drive that leads to the excitation of the state \(\left\vert 2\right\rangle\). One then obtains

where⁷⁴

This formula can be generalized to arbitrary n-photon resonances with \({\Delta }_{\bar{m}}^{(n)}=\Delta -\bar{m}\Omega\) and \({G}_{\bar{m}}^{(n)}\propto {({J}_{\bar{m}}\left(\zeta /\Omega \right))}^{n}\). We conclude that when the rescaled detuning matches the condition for the nth multiphoton resonance, and if the rescaled drive \({F}_{\bar{m}}\) is strong enough, an additional dip appears. Therefore, we can treat each of the multiphoton resonances for each LZSM dip as a yet separate phenomenon.

As sketched in Fig. 4e–g, at the multiphoton resonance, i.e., at Δ = mΩ + (n − 1)χ, the states \(\left\vert 0\right\rangle\) and \(\left\vert 2\right\rangle\) can satisfy the conditions for the development of constructive interference. In other words, around each of the main LZSM dips, and for large enough drive amplitude, several multiphoton resonances emerge with the same characteristics as those shown in Fig. 2c–e.

Weak modulation case

When Ω ≪ ∣χ∣ (weakly modulated case), instead, one can capture the system’s behavior around the second multiphoton resonance via the Hamiltonian in Eq. (6), with G⁽²⁾ = F²/χ representing the effective drive between the vacuum and the state \(\left\vert 2\right\rangle\) if ζ = 0. Removing the modulation using the same approximation as in Eq. (3) leads to an equation identical to Eq. (7), where now

This formula can be generalized to arbitrary n-photon resonances, with \({\Delta }_{\bar{m}}^{(n)}=\Delta -\bar{m}\Omega /n\) and \({G}_{\bar{m}}^{(n)}\propto {J}_{\bar{m}}\left(n\zeta /\Omega \right)\). Thus, for detunings close to the nth multiphoton transition, a new LZSM interference pattern should emerge, characterized by an effective modulation frequency Ω/n. It is this scaling that differentiates the weakly and strongly modulated cases, c.f. Figs. 5a, e and 4d. While previously, for Ω ≫ ∣χ∣, the standard LZSM sidebands were dressed by Kerr multiphoton resonances, we now find that, for Ω ≪ ∣χ∣, each Kerr n-th multiphoton resonance is dressed by LZSM sidebands with effective modulation frequency Ω/n.

Through the figure, we set the drive input power to P_in = −128 dBm. a–d Both modulation strength ζ and frequency Ω are swept together to maintain a constant ratio ζ/Ω. This choice ensures that the effective drives in Eqs. (3) and (7) are kept constant. This allows enhancing the visibility of the transition between the strongly- and weakly-modulated cases. a Sketch of the results of Eq. (7) for the transition \(\left\vert 0\right\rangle \to \left\vert 1\right\rangle\) (purple, labeled \(\left\vert 1\right\rangle\)) and \(\left\vert 0\right\rangle \to \left\vert 2\right\rangle\) (green, labeled \(\left\vert 2\right\rangle\)). For \(\left\vert 1\right\rangle\), the pattern radiates from Δ = 0 with frequency modulation Ω. For the multiphoton transition to \(\left\vert 2\right\rangle\), the LZSM interference pattern is centered at Δ = χ and scales with Ω/2. b Simulation of ∣S₂₁∣ as a function of Δ and Ω with the full quantum model in Eq. (2) (see Supplementary Information for details on the simulation method), having fixed ζ/Ω = 0.86. The corresponding measurement is shown in (c) and perfectly overlaps with the results of the numerical simulation. The black arrows in (a–d) mark the position of two avoided crossings, where the “bare levels” in (a) interact and hybridize in (b–d). The crossings are further highlighted by the solid (associated with \(\left\vert 1\right\rangle\)) and dashed (\(\left\vert 2\right\rangle\)) lines in (b). The amplitude of the different avoided crossings can be controlled by modulating the Bessel functions \({J}_{\bar{m}}(n\zeta /\Omega )\) as shown in (d), where a larger ratio ζ/Ω is chosen. e As in panel (a), the sketch of the results of Eq. (7) for Ω = ∣χ∣ and as a function of Δ and ζ. In this “bare picture”, the two independent LZSM interference patterns scale with Ω and Ω/2 for \(\left\vert 1\right\rangle\) and \(\left\vert 2\right\rangle\), respectively. f Full quantum simulation and g corresponding measurement of ∣S₂₁∣ for Ω = ∣χ∣. The position of some LZSM resonances in the bare picture is superimposed in (f) as a guideline for the eye. h Repeating the measurement for Ω = 1.5∣χ∣, we observe line splittings, indicating a modulation of the coupling between different Floquet states.

For the device under consideration, accessing the weakly modulated case would require κ ≪ Ω/n to distinguish between the different LZSM dips. To better resolve this feature, we propose the following driving scheme. We fix the ratio ζ/Ω to have a constant effective drive, according to both effective theories in Eqs. (3) and (7). We then increase Ω and ζ, to cross from the weakly modulated ∣χ∣ > Ω to the strongly modulated case ∣χ∣ < Ω. This is shown in Fig. 5b–d where, for small Ω, we distinctly see the expected LZSM dips associated with the second multiphoton resonance \(\left\vert 0\right\rangle \to \left\vert 2\right\rangle\) and with a slope Ω/2. We note that such two-photon LZSM transitions were recently reported in a linearly-modulated three-level system⁷⁵, with a similar factor two in the LZSM velocity compared to regular single-photon LZSM transitions.

Non-perturbative regime

The weak- and strong-modulation regimes have very different scales from each other [c.f. the effective models in Eqs. (8) and (9)]. We thus expect that there is a non-perturbative passage from weak- to strong modulation through some effective interaction, and the transition between these two regimes cannot be explained using any of the two effective theories alone.

Particularly interesting are the values of Ω ≃ n∣χ∣, where the system passes from the weak- to the strong-modulated case for a specific state \(\left\vert n\right\rangle\). At these values, it is possible for a n-photon resonance to exactly match the LZSM dips of a different m-photon resonance. We observe the signatures of avoided level crossings between resonances in Fig. 5b, c, indicating that the \(\left\vert 0\right\rangle \to \left\vert 1\right\rangle\) and \(\left\vert 0\right\rangle \to \left\vert 2\right\rangle\) resonances interact through the action of an effective emergent coupling. In this sense, these different resonances constitute a controllable synthetic Floquet space, where changing Ω and ζ allows selecting an effective interaction between these multiphoton resonances. This is also evident in Fig. 5d, where the ratio Ω/ζ is changed, leading both to different interference patterns and different splittings between the Floquet states.

To further highlight an example of these non-perturbative effects, in Fig. 5f, g we fix Ω = ∣χ∣. First, we numerically simulate the interplay of these effects in Fig. 5f. We predict a partial overlap between the second multiphoton transition \(\left\vert 0\right\rangle \to \left\vert 2\right\rangle\) with the first LZSM dip associated with the \(\left\vert 0\right\rangle \to \left\vert 1\right\rangle\) transition at Δ = −Ω. For increasing modulation strength ζ, the LZSM structure predicted by Eq. (7) is observed, although strongly deformed compared to the prediction of the effective model due to the presence of the LZSM lobe associated with the \(\left\vert 0\right\rangle \to \left\vert 1\right\rangle\) resonance. These theoretical predictions are completely recovered in the data in Fig. 5g. Finally, in Fig. 5h we fix Ω = 1.5∣χ∣, and we observe a line splitting of several resonances, indicating again the merging and interaction between \(\left\vert 0\right\rangle \to \left\vert 2\right\rangle\) and \(\left\vert 0\right\rangle \to \left\vert 1\right\rangle\) transitions. For larger drive amplitudes (not shown), the system shows an extremely rich structure that cannot be simply assigned to any of these original phenomena. Note also the asymmetric nature of the interference pattern, determined by the negative sign of the Kerr nonlinearity.

LZSM beyond the qubit approximation: Duffing regime

Finally, we investigate the Duffing regime κ > ∣χ∣ for a drive amplitude sufficiently large to deviate from the linear regime (see Eq. (24)). For the intermediate drive amplitudes shown in Fig. 6a, b, the various dips are well separated despite showing an asymmetric bending of ∣S₂₁∣. When compared with Fig. 2g, h, we observe a similar deformation of the transmission dips. Therefore, we assign this feature to the emergence of bistability triggered by the competition between detuning and Kerr nonlinearity. For these parameters, we find that the formula in Eq. (3) captures the deformation of the dips, as discussed more in detail in the Supplementary Information. Thus, the system behaves as a collection of independent Duffing oscillators and the overall effect of the modulation is to rescale the drive amplitude F of each sideband.

a–c Measured magnitude of S₂₁ vs Δ and ζ for increasing drive power P_in. The dashed color lines refer to the values of ζ chosen for panels (d–f). The star in (c) indicates the value where the analysis of chaos is performed in Fig. 7d. d–f Measured S₂₁ as a function of detuning and for three specific values of ζ/2π: 41.3 MHz (green curve), 101.1 MHz (orange curve), and 161.2 MHz (purple curve). The black superimposed curves are the results of the numerical simulation of the full quantum model for the parameters in Table 1 (detailed in “Methods”). The modulation frequency is set to Ω/2π = 30 MHz. The systematic discrepancy in the position of the dips between theory and experiments is due to the nonlinear dependence of the modulation of the flux amplitude discussed in the Supplementary Information.

When the driving power is further increased in Fig. 6c, several of the neighboring LZSM dips eventually overlap. This case cannot be simply captured as separated LZSM interferences, and it is qualitatively different from all the previously studied cases. The simplified picture of Eq. (3) thus breaks down, and the system becomes multimodal and behaves as a set of interacting nonlinear cavities. Nonetheless, the full simulation of the quantum Floquet model matches the data in all regimes, as shown in Fig. 6d–f.

As detailed in the methods section, the merging of several modes and the qualitative change in the system’s behavior can be associated with the emergence of dissipative quantum chaos. At weak pump power, the LZSM dips correspond to distinct Fourier modes, each characterized by its own frequency. As the pump power increases, these modes begin to interact and merge, analogous to phenomena observed in strongly driven resonators^45,46, thereby suggesting the onset of chaos in the Floquet system. Classical chaos is characterized by a system’s sensitive dependence on initial conditions, often quantified by a positive Lyapunov exponent⁷⁶. On the other hand, the characterization of quantum chaos often relies on the spectral properties predicted by random matrix theory^77,78,79,80. In open quantum systems, quantum chaos can be extended through the analysis of the Liouvillian superoperator, which governs the dynamics of the density matrix and provides insights into integrability and chaos. The complex spacing ratio is an efficient criterion for distinguishing between integrable and chaotic regimes, as it assesses the distribution of spacings between eigenvalues⁸¹. However, as shown in the Methods section, applying the usual complex spacing ratio criterion to Floquet systems fails to capture the nuances of dissipative quantum chaos. Instead, we generalize the spectral statistics of quantum trajectories (SSQT) criteria introduced in ref. ⁴⁵ to Floquet systems, demonstrating its relevance and correctness in identifying chaotic behavior (see Fig. 9). This refined approach allows for a precise analysis of the system’s dynamics, accurately pinpointing the transition to chaotic phases.

Merging of sidebands: signature of dissipative quantum chaos

Utilizing the novel criterion developed in the methods for Floquet dissipative quantum systems, we demonstrate that the parameters predicting chaos in the model correspond precisely with those where the sidebands merge.

We investigate LZSM interference for three values of modulation strength ζ as a function of the driving power P_in, as shown in Fig. 7a–c. As in all the other experiments presented in this work, measurements are done in a steady state and are thus independent of the initial conditions of the system. At low driving power, we find a linear regime where m dips appear separated by the frequency Ω/2π = 30 MHz and with visibility given by Bessel functions \({J}_{\bar{m}}(\zeta /\Omega )\). This regime is remarkably similar to that of several nonlinear modes separated by the same frequency Ω. As the driving power increases, each of these dips initially follows the typical Duffing behavior of a single resonator, as already mentioned. For high enough input power, however, these individual dips disappear and merge, leading to a very broad response of the system. At this point, one completely loses the notion of individual synthetic modes and their bistability.

a–c Measurement of ∣S₂₁∣ for increasing drive power, fixed frequency modulation Ω/2π = 30 MHz and increasing ratios ζ/Ω ≈ 1.4 (a), ζ/Ω ≈ 2.6 (b), and ζ/Ω ≈ 3.8 (c). For low input powers the system behaves as a collection of noninteracting nonlinear modes, each one well separated from the others. For larger values of P_in, the system enters a phase characterized by a single broad response where the notion of isolated mode is lost. Such a response can be observed in multimode nonlinear systems and has been associated with a transition from integrability to dissipative quantum chaos⁴⁵. To show that this is indeed a dissipative quantum chaotic phase we plot in (d) the histogram of the probability density p(s) of the level spacings s obtained by diagonalizing the Floquet Liouvillian in the broad-response region indicated by the star in (a). Parameters are set to ζ/2π = 41.3 MHz, Ω/2π = 30 MHz, Δ = −1.1Ω and F/2π = 49.5 MHz (P_in ≈ −105 dBm). The cutoff in the Hilbert space is set to 90. The solid black (orange) curve represents the ideal Poisson (Ginibre) distribution given by Eqs.(15) and (16) associated with integrability (chaos).

The merging and broadening of the dips of the scattering coefficient ∣S₂₁∣ in Fig. 7a occur for an input power P_in ≈ −108 dBm, which coincides with the point where dissipative quantum chaos emerges according to our theory, as reported in Fig. 9 of the Methods section. In Fig. 7d, we plot the probability density of the level spacings obtained by diagonalizing the Floquet Liouvillian for the parameters indicated by a star in Fig. 7a where LZSM dips have merged. It has been shown that integrable systems exhibit Poisson-distributed level spacings, indicating no level repulsion, while chaotic systems follow Ginibre statistics, characterized by level repulsion and non-Hermitian random matrix behavior^45,80. We find that, upon the merging of the Floquet modes, the Floquet Liouvillian level statistics conform to the Ginibre distribution [see Eq. (16)], indicating a clear transition to the dissipative quantum chaos regime.

The onset of the chaotic phase can also be controlled by tuning the spacing between LZSM resonances through the modulation frequency Ω, as shown in the Supplementary Information. For instance, the separated bistable regions of Fig. 6b would start overlapping by decreasing Ω, potentially resulting in a chaotic state.