Classical and quantum chaos in chirally-driven, dissipative Bose-Hubbard systems

We study the dissipative Bose-Hubbard model on a small ring of sites in the presence of a chiral drive and explore its long-time dynamical structure using the mean-field equations and by simulating the quantum master equation. Remarkably, for large enough drivings, we find that the system admits, in a wide range of parameters, a chaotic attractor at the mean-field level, which manifests as a complex Wigner function on the quantum level. The latter is shown to have the largest weight around the approximate region of phase space occupied by the chaotic attractor. We demonstrate that this behavior could be revealed via measurement of various bosonic correlation functions. In particular, we employ open system methods to calculate the out-of-time-ordered correlator, whose exponential growth signifies a positive quantum Lyapunov exponent in our system. This can open a pathway to the study of chaotic dynamics in interacting systems of photons.


INTRODUCTION
Methods to confine and control photons using solid-state materials have opened a pathway to the generation of novel liquids of light 1 , exhibiting a wide range of phenomena, from superfluidity 2 to topological states 3 . Coupled cavity arrays evolved into a valuable experimental framework for the study of lattices of photons in the presence of interactions, drive, and dissipation [4][5][6] . Such systems benefit from the exquisite fabrication techniques of solid-state materials as well as from the high coherence properties of the photons, making them an attractive venue for future applications 7 .
Even a small number of driven coupled cavities can exhibit rich behavior: a single cavity gives rise to regions of bistability and to critical slowing down 8 ; cavity dimers demonstrate Josephson oscillations 9 and parametric instabilities 10 and could be driven into a more exotic quantum limit-cycle if the dissipation is further engineered 11 ; the levels of a six cavity ring were demonstrated to be individually excitable, generating chirality, with their population exhibiting hysteretic behavior, indicative of a bistability 12 ; and many more. These phenomena explore at most a quasi-two-dimensional manifold of phase space. Yet, on the classical, mean-field level such coupled cavities are described as driven-dissipative coupled nonlinear oscillators, that are in principle capable of exhibiting chaotic behavior [13][14][15] . However, it remains unclear what are the key ingredients which are required to generate chaotic behavior in such systems. Further, once the interactions increase and the system crosses into a more quantum regime, what implications does the emergence of classical chaos carry for the quantum behavior of the system?
In this paper, we consider a ring of three coupled cavities, driven by finite orbital momentum light [16][17][18] (illustrated in Fig. 1). We demonstrate that as a direct result of the chirality of the applied drive, the system exhibits robust chaotic dynamics over a large part of its mean-field phase diagram. On this level, the chaos is manifested as a dense, compact attractor in phase space. For the quantum problem, a full quantum trajectory approach 19 reveals a complex Wigner function, indicative of the underlying chaos. The g (1) correlator places the system on the chaotic attractor, while the g (2) correlator reveals signatures of the chaotic behavior. We find that, most prominently, the out-of-time-ordered correlator (OTOC) reveals the onset of chaos via an exponential sensitivity to the application of a weak perturbation. These predictions could pave the way to the study of the interplay of classical and quantum chaotic dynamics and their properties in coupled cavity arrays of photons.

Proposed setup and model
We consider the Bose-Hubbard Hamiltonian for M sites (cavities) on a ring in the presence of drive and dissipation, with a driving term that is coherent with respect to both frequency and momentum, j¼0 e ipxjb j þ e Àipxjb y j : (1) It prescribes the non-Hermitian effective HamiltonianĤ eff ¼Ĥ À j¼0b y jb j which is the deterministic part of the Lindblad 1 evolution; as well as the so-called "recycling" term, Γ P jb jρb y j . In the following we use the rescaled parameters u = U/J, d = D/J, ω = Ω/J, and γ = Γ/J.

Mean-field phase diagram
We start by analysing the mean-field equations of motion, derived from the hamiltonianĤ eff . We analyze the fixed points using the ansatz that only a single-mode participates in the dynamics, j e ipxjb j , with the operators replaced by c-numbers. The regimes of existence and stability of the resulting fixed points are delineated by black lines in Fig. 2a, b (note that the shading has different meaning as we explain in the next paragraphs). These include 'dim' (denoted by I and II, depending on the existence and stability of the other fixed points) and 'bright' (denoted by III) regimes which are each dominated by the effect of a single fixed point with low and high bosonic occupation, respectively, as well as a bistable regime (denoted 'Bi') where both fixed points are stable, and a region where no fixed point is stable (denoted 'NFP').
The difference between a uniform and a chiral drive is already exposed in this mean-field level. These two types of drive lead to the two dramatically different phase diagrams, presented in Fig. 2, with panel (a) describing the uniform case and (b) the chiral case. In between the dim and bright regimes, for the uniform drive and for the chiral drive at low amplitudes, there exist bistability regimes. This picture changes considerably for the chiral drive at higher amplitudes, where a large regime where no fixed point remains stable appears in the phase diagram.
To gain further intuition into the nature of this regime, we note that the actual steady state the system attains is decided by the  initial condition. When we start from the vacuum state 0 j i with zero bosonic occupation, the resulting steady state is depicted by the shading in Fig. 2a, b. In particular, in most of the bistable regime, the system reaches the dim fixed point. In the region NFP it turns out that the system reaches no particular fixed point, but the long-time dynamics remains confined in a certain dense region of phase space. In order to analyze it, we have to go beyond the single-mode ansatz (see the "Methods" section).
We next relax the requirement that we start from the 0 j i state and consider the two-dimensional phase plane of initial conditions associated with b p . Each initial condition is color-coded according to the attractor the system eventually reaches. The result is described in Fig. 2c-h, and defines the basins of attraction of the relevant attractors. In addition, the long-time trajectories of b p are superimposed in red.
The spiral shape of the basins was highlighted in ref. 20 for the case of a single cavity. Our case is richer since the higher dimensional phase space allows for the possibility of classical chaos, which in our system is accessible due to the operation of a chiral drive. The onset of chaos is marked by an exponential sensitivity to initial conditions and is characterized by a positive (classical) Lyapunov exponent. For the NFP regime, we find that the system ends up in either a limit cycle or in a strange attractor (SA), the latter depicted by the red shape in Fig. 2e-h. In fact, the SA also exists beyond the region marked by NFP, and when the vacuum is contained in its basin of attraction the system will flow into this steady state, see e.g., Fig. 2d. Technically, we identify the limit cycles by the occurrence of finite frequencies in the Fourier transform of the time-dependence of the occupation combined with a vanishing maximum Lyapunov exponent.
Once the interaction in the system increases the system will characteristically hold lower bosonic occupations and can cross into the quantum regime. In the next section, we will investigate the quantum behavior of the system within and in the vicinity of the NFP region, and explore its properties.

Quantum approach to chaos
We now turn to study the full quantum evolution of Eq. (2) using the method of quantum trajectories 19 , and compare with the classical equation of motion analysis. The quantum trajectories are visualized in Fig. 3 using a histogram, which depicts the formation and interplay of the various fixed points and their exploration by the quantum dynamics. In Fig. 3a-c, the trajectories explore the vicinity of the bright fixed point but as the frequency increases a larger portion of the chaotic basin is traversed quantum mechanically. Figure 3c, d depicts a regime where the full classical chaotic basin is explored. Curiously, in Fig. 3d the trajectories initially tend towards the newly nucleated dim fixed point but eventually flow past it and end up in the chaotic attractor. Figure 3e, f shows the coexistence of two attractors, a chaotic one and a dim fixed point, as the weight of the former decreases with increasing frequency. The state the system attains at long times can be characterized by the Wigner function Wðβ; tÞ in the momentum-space coherent state representation β j i, of which we display only the twodimensional cut with finite β p (and β k = 0 for k ≠ p) in Fig. 4a-e (see the "Methods" section for details of the calculation). When the system flows into a fixed point it generates a Wigner function that is localized around this fixed point, forming an approximate coherent state, see e.g. Fig. 4e. In contrast, in the NFP regime, the Wigner function exhibits a pattern of islands of positive and negative values which aggregate around the SA, see Fig. 4b-d. When projected onto the resonant plane, we introduce the Wigner function which characterizes the associated reduced density matrix, which becomes strictly positive, see Fig. 4f-j. Remarkably, it admits the greatest support around the region of the classical SA in phase space, Fig. 4h, i 21 .
Next, we consider photonic correlation measurements and weigh their relevance in unveiling the presence of quantum chaos in our system, starting with the standard photonic correlators g (1) and g (2) , obtained by the quantum trajectories procedure (see "Methods" for details).
First, the correlator g ð1Þ p ðtÞ ¼ hb y p ðtÞb p ðtÞi is associated with the photonic occupation and is presented in the inset of Fig. 5a. As the frequency increases, it is shown to reach a steady-state value that progressively deviates from the value of the bright fixed point towards a mean value associated with the chaotic attractor (see the main panel of Fig. 5a). When the frequency further increases, the value of the g (1) correlator rapidly drops into the dim fixed point.  correlator becomes low, i.e., at low and high frequencies.
As it is proportional to the uncertainty of the occupation ΔN 2 p ðtÞ ¼ hb y p ðtÞb p ðtÞb y p ðtÞb p ðtÞi À hb y p ðtÞb p ðtÞi 2 , we plot the latter instead at long times for the resonating momentum p in the main panel of Fig. 5b. In addition, the fluctuations of the quantum trajectories around the mean value could indicate the presence of the SA if the quantum trajectories are individually measured 22 .
The most prominent indicator of quantum chaos in our system turns out to be the OTOC 23,24 , defined by, Oðt; τÞ ¼ Àh½Qðt þ τÞ;PðtÞ 2 i; are the dimensionless position and momentum operators associated withb p . The OTOC generalizes the concept of a Lyapunov exponent to the quantum regime, i.e., lim t!1 Oðt; τÞ $ e λqτ where λ q is the quantum Lyapunov exponent and τ is in some intermediate regime: for small enough τ the commutator is set by its equal-time value while for large enough τ it saturates 25 . The saturation value of lim t;τ!1 Oðt; τÞ is related to the size of the available phase space in the chaotic basin 26 and its relation to chaos was demonstrated for the periodically kicked Dicke model 27 . A version of the OTOC was measured in a trapped ion quantum magnet 28 . In the "Methods" we detail a procedure to calculate it for a dissipative system.
The result, presented in the inset of Fig. 5c, exhibits the anticipated exponential increase within the chaotic regime.  Remarkably, we find a non-monotonic behavior of the saturation value of the OTOC as function of the drive frequency ω (see the main panel of Fig. 5c). It is interesting to note that a positive quantum Lyapunov exponent is detected for lower values of ω as compared to the classical phase diagram. At higher frequencies, once the frequency crosses into the classical dim fixed point regime, a dramatic reduction in the quantum Lyapunov exponent occurs. It therefore turns out that the quantum chaotic regime extends beyond the classical one into the lower frequency regime; see main panel of Fig. 5c where the classical chaotic regime is depicted by blue shading. The higher frequency regime is marked by a coexistence of the dim and chaotic fixed points. The main panel of Fig. 5d displays, using colored points, the entanglement entropy associated with the resonant mode calculated for individual trajectories, S p ðtÞ ¼ ÀTrρ p ðtÞlogρ p ðtÞ Â Ã , withρ p ðtÞ ¼ Tr k≠p ψðtÞ j i ψðtÞ h j the reduced density matrix. For each trajectory S p (t) is well defined, and the average over trajectories S p ðtÞ is displayed as a black line. In addition, the variance of the entanglement entropy over trajectories, quantified by δS p ¼ ðS

DISCUSSION
We have demonstrated that the driven-dissipative Bose-Hubbard model with a small number of sites on a ring, once chirally-driven, could lead to the generation of chaotic regimes in both the meanfield analysis and at the full quantum level. To characterize the quantum chaos, we have calculated the Wigner function on a cut through the 6-dimensional phase space and projected it onto the 2-dimensional resonant plane, where it is shown to spread over the approximate region of the classical chaotic phase space. Next, we extended the Monte-Carlo method to calculate the OTOC in driven-dissipative systems (see "Methods" section), and demonstrated its sensitivity to the onset of quantum chaos in our model. Finally, we calculated the measurable correlators g (1) and g (2) and extracted the relevant signatures of the chaotic attractor. The results could prove relevant for photons in coupled cavity arrays and may lead to a study of chaos in such systems.

Mean-field analysis
To study the mean-field fixed points we employ above the ansatz that only the single mode that shares its momentum with the drive is macroscopically occupied, i.e. b k = 0 for k ≠ p. The equation of motion for b p is given in this case by Denoting the steady-state solution as b p ðt ! 1Þ ¼ ffiffiffiffiffi ffi N 0 p e Àiθ we arrive at the following equation for the occupation of the condensate, n 0 = N 0 /M, and its phase θ ¼ sin À1 ðγ ffiffiffiffiffi n 0 p =2dÞ. Equation Here L q ¼ Àϵ 0 ðqÞ þ iσ y ϵ y þ σ z ϵ z ðqÞ, where ϵ 0 ðqÞ ¼ sinðqÞ sinðpÞ, ϵ y = n 0 u, ϵ z ðqÞ ¼ 2n 0 u À cosðqÞ cosðpÞ À ω, σ i are the Pauli matrices and μ q is the Bogoliubov dispersion, given by μ q ¼ Àϵ 0 ðqÞ ± ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ϵ 2 z ðqÞ À ϵ 2 y q . For a solution to be stable we require δb j to decay with time, thus stability is contingent upon |Imμ q | ≤ γ/2 for all values of q ¼ 2π M m, m ∈ {0, 1, . . , M − 1}. This is an extension 1 of the more familiar case without loss (γ = 0).
In the main text, we describe the resulting phase diagram, which is presented in Fig. 2a, b. Here we describe the phase diagram in terms of existence and stability of the fixed points. In region I (III) there exists only a single fixed point that is also stable: the dim (bright), ν = 1 (ν = 3), is the low (high) occupation fixed point (the intermediate fixed point, ν = 2, is always unstable). Whereas, in region II there exist three fixed points of which only the dim (ν = 1) one is stable. In addition, there are the bistable regime (denoted by "Bi") in which both ν = 1 and ν = 3 are stable, and the regime where no fixed point is stable (denoted by "NFP") in which the system flows into the chaotic attractor, which is not revealed by the linear analysis above.

Quantum analysis
We calculate the Wigner function displayed in Fig. 4 using the formula where ρ β;η ðtÞ ¼ β À η 2 jρðtÞjβ þ η 2 . Here β j i represents a product of momentum coherent states, i.e., β j i¼ N k β k j i; satisfyingb k β j i ¼ β k β j i. Expressing the density matrix using momentum occupation states n j i¼ N k n k j i, we get the following transform where F n;n 0 ðξÞ ¼ ffiffiffiffiffiffiffiffiffi n!n 0 ! p ξ n 0 Àn U Àn; n 0 À n þ 1; jξj 2 ; with U(a, b, c) denoting the confluent hypergeometric function. The density matrix ρ n;n 0 ðtÞ is then calculated employing 2000 quantum trajectories. In the main text we considered two types of Wigner functions. The first was the full 6-dimensional Wigner function of Eq. (8) of which we displayed only the 2-dimensional cut representing the resonant plane (β k = 0 for k ≠ p) at time Jt m (see Fig. 4a-e). Second, we consider the reduced density matrix relevant to the resonant momentum, i.e, the density matrix in Eq. (8) is replaced withρ p ðtÞ ¼ Tr k≠pρ ðtÞ (see Fig. 4f-j). All quantum trajectory simulations are performed by truncating each of the three bosonic modes to 30 occupations states. Using the same quantum trajectory Monte-Carlo formalism, we calculate the different photonic correlation functions as we now describe.
To calculate the g ð2Þ p ðt; τÞ ¼ hb y p ðtÞb y p ðt þ τÞb p ðt þ τÞb p ðtÞi=hb y p ðtÞb p ðtÞi 2 correlator, we employ the usual method of first enforcing a quantum a jump at t and then finding the evolution along τ following the same Monte-Carlo procedure.
Calculating the OTOC in a dissipative-driven system is more challenging. The OTOC can be expanded into a sum of four two-time correlators with both forward and backward time evolutions. For each of these correlators, let t 1 be the time associated with the rightmost operator, and let t 2 be the other time. We proceed according to the following steps: (i), we evolve the wavefunction to t 1 (i.e. either t or t + τ); (ii), we apply the operators associated with t 1 on both the evolved ket and bra or, where this is not possible, we define four helper states following Mølmer et al. 19 ; next, (iii), we apply either forward or backward evolution in time from t 1 to t 2 , as necessary (i.e., forward evolution from t 1 = t to t 2 = t + τ or backward evolution from t 1 = t + τ to t 2 = t). The latter is done according toĤ ! ÀĤ but with the dissipative part of the Lindblad evolution remaining the same 30 ; finally, (iv), we find the expectation value of the operator or product of operators associated with time t 2 .

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.