Rui Lin

Paolo Molignini

Axel U J Lode

R Chitra

Pathway to chaos through hierarchical superfluidity in a cavity-BEC system

We explore the role of atomic correlations in a harmonically trapped Bose-Einstein condensate coupled to a dissipative cavity, where both the atoms and the cavity are blue detuned. Using a genuine many-body approach, we extract density distributions and many-body correlations to unveil a pathway to chaos at large pump power through a hierarchical self-organization of the atoms. The atoms transition from a single-well optical lattice to a double-well optical lattice. Correlated states of the atoms emerge and are characterized by local superfluid correlations in phases which are globally superfluid or Mott insulating. Local superfluid-Mott transitions are precluded by a dynamical instability to chaos which occurs via strange attractors. Our findings thus explain the mechanism behind the dynamical instabilities observed in experiments.

Introduction -Experimental advances in the past decade have heralded a new era in light-matter hybrid systems, where quantum light is used to engineer correlated phases of matter. In the solid state realm, coupling to light has been used to activate phases of matter like ferroelectricity [1] and superconductivity [2]. In the quantum engineering domain, cavity-QED systems with their tunable light matter couplings provide a versatile platform to realize hybrid correlated quantum fluids like polaritons [3,4], permit the encoding of qubits through photons [5], and generate entangled quantum many-body states for quantum computation [6].

A landmark example of a light induced phase is superradiance [7][8][9][10] in strongly coupled cavity-Bose-Einstein condensates (BEC) systems where the atoms in the BEC self-organize onto a lattice dynamically generated by the cavity field [11][12][13]. Cavity-BEC systems also host complex phases like Mott insulators [14][15][16][17][18], supersolids [19][20][21], and spin textures [22][23][24][25][26]. They additionally permit the simulation of many-body Hamiltonians having no solid-state counterparts, like spin models with both short and long range interactions [2,27], and the realization of exotic collective magnetic phenomena [28].

Though the red-detuned regime of the cavity is well understood, the blue-detuned regime remains germinal. Theoretical studies highlight the potential of blue detuning to stabilize dynamical states like phase slippage [30], limit cycles and chaos [30,31], and topology in a gas of fermions [32]. A preliminary observation of dynamical instabilities was reported recently in a first study of a blue-detuned cavity with blue-detuned atoms [33]. Motivated by this study, we explore the interplay between correlations and dynamical instabilities in a realistic harmonically trapped and weakly interacting BEC coupled to a blue-detuned cavity [see Fig. 1(a)]. Beyond the Gross-Pitaevskii mean-field limit, this interplay results in a spontaneous deformation of the emergent superradiant optical lattice into a double-well optical lattice, generating new correlated phases of the atoms. Our results are summarized in the rich phase diagram as shown in Fig. 1(b).

Model and method -We consider a cavity-BEC sys- Phase diagram. For cavity detuning 0 < δ < 1/2, the system transitions from the normal phase to the dynamically unstable region via superradiance with increasing pump rate A ∝ η.

The strongly correlated phases are self-organized superradiant (SSF), a self-organized Mott insulator (SMI), a self-organized dimerized superfluid (SDSF), and a self-organized secondorder superfluid (2-SSF) phase. The orange line delineates superfluid phases and globally Mott insulating phases, while the green line marks the deformation to a double-well potential. Pronounced sensitivity to the ramping protocol is seen in the hatched dark green region. At higher η, the system is dynamically unstable to the formation of strange attractors followed by chaos. In presence of a trap and atomic fluctuations, the strange attractors only exist in the region represented by the thick gray line, while the thin dashed section represents a direct transition to the chaos. The dimensionless detuning δ and the potential strength A are normalized with respect to N U0 and √ ωR, respectively. (c) Sketch of the SSF, SMI, SDSF, and 2-SSF phases.

tem with N bosonic atoms of mass m. From a computational perspective, since the physics of the system (g,j), the color code follows the functionln(1g (1) ). The double-well splitting is seen in the central lattice site but not in the other two lattice sites, because only M = 4 orbitals are used in the numerical simulations [29].

does not qualitatively depend on the dimensionality, we study a one-dimensional version of the model and later discuss the validity of the results obtained for the twodimensional system. In the rotating frame of the pump laser, the system is described by the following coupled equations of motion for the cavity expectation value α and the atomic field operators Ψ( †) (x) [11,34]:

i∂ t Ψ(x) = - ∂ 2 x 2m + g | Ψ † (x) Ψ(x)| 2 + 1 V trap (x) (1a) + U 0 cos 2 (k c x)|α| 2 + η cos(k c x)(α + α * ) Ψ(x). ∂ t α = [i∆ c -iN U 0 B -κ] α -iηN θ, (1b)
Here, g is the weak interatomic interaction, U 0 > 0 is the single photon light shift for blue-detuned atoms, η is the effective pump rate, ∆ c is the cavity detuning, k c is the wave vector of the cavity field, and κ is the cavity dissipation rate. The atoms are confined by a harmonic trap

V trap (x) = 1 2 mω 2 x x 2 .
Since we are interested in regimes far from the normal phase-superradiant phase boundary where the cavity is in a coherent state and the fluctuations are negligible [18,[35][36][37], the cavity field can be replaced by a complex number α. The variables θ and B in Eq.(1b) are the order parameters associated with superradiance:

θ = 1 N dx Ψ † (x) Ψ(x) cos(k c x), (2a)
B = 1 N dx Ψ † (x) Ψ(x) cos 2 (k c x). (2b)
In the limit κ ω R ≡ k 2 c /2m where ω R is the recoil energy, the cavity field can be adiabatically eliminated via ∂ t α = 0. A further mean-field treatment of the atoms where both contact interactions between the atoms and the harmonic trap are neglected reduces Eq.( 1) to an effective single-particle Hamiltonian for the atoms [11] Ĥ

(1) = - ∂ 2 x 2m + A 2 ω R [2(δ -B)θ cos(k c x) + θ 2 cos 2 (k c x)], (3)
with

A = ηN √ U 0 / [(∆ c -N U 0 B) 2 + κ 2
]ω R the dimensionless overall effective potential strength and δ = ∆ c /N U 0 the dimensionless cavity detuning. The cavity dynamically creates an optical lattice potential comprising two terms: a cos 2 (k c x) term describing the standing wave of the cavity field and a cos(k c x) term stemming from the interference between the cavity and pump fields. The amplitudes of both sinusoidal terms are determined by the instantaneous state of the atoms via θ and B.

Mean-field analysis predicts that the self-organization of atoms takes place in the regime ∆ c < N U 0 /2 [11,12,38]. For blue-detuned atoms, U 0 > 0, self-organization occurs for both red-and blue-detuned cavities. When the cavity is red-detuned, ∆ c < 0, the atoms are localized at the lattice sites x n = nπ/k c with all n ∈ N either even or odd [11,12,38]. However, this lattice is not always stable in a blue-detuned cavity even when 0 < ∆ c < N U 0 /2. To see this, consider the large pump rate limit A 1, where |θ| = B = 1 and the effective potential is positively proportional to 2sgn(θ)(δ -1) cos(k c x) + cos 2 (k c x). For 0 < δ < 1/2, the standing wave term [cos 2 (k c x)] becomes relevant and spontaneously forms a double-well optical lattice, where each lattice minimum x n is replaced by a pair of degenerate adjacent minima. It has been pro-posed to use this double-well lattice to realize the Su-Schrieffer-Heeger model in a cavity-fermion system [32]. This analysis provides the first glimpse of intrinsically different physics in the blue-detuned region.

Beyond the mean-field limit, the combination of the double-well optical lattice and weak atomic interactions leads to a rich phase diagram comprising interesting correlated phases of matter. We investigate the full phase diagram of the blue-detuned system described by Eq. ( 1) using the Multi-Configurational Time-Dependent Hartree method for indistinguishable particles (MCTDH-X) [39][40][41][START_REF] Lode | MCTDH-X: the time-dependent multiconfigurational hartree for indistinguishable particles software[END_REF]. The simulated N = 50 atoms are initialized in a state which is almost a Thomas-Fermi cloud, and the many-body state of the atoms coupled to the cavity is propagated in real time. The pump rate η is linearly ramped up at fixed detunings to reach its desired value. The simulation parameters correspond to those realized experimentally in Ref. [12] and are given in detail in Ref. [29]. The phase diagram is extracted from the observables: θ, B, the position space density distribution ρ(x) = Ψ † (x) Ψ(x) /N , the momentum space density distribution ρ(k) = Ψ † (k) Ψ(k) /N , and the Glauber one-body correlation function g (1) [START_REF] Glauber | [END_REF]44].

(x, x ) = Ψ † (x) Ψ(x ) / Ψ † (x) Ψ(x) Ψ † (x ) Ψ(x )
Results -Our results for the blue-detuned cavity-BEC are summarized in the phase diagram [Fig. 1(b)] along with a schematic representation of the different phases [Fig. 1(c)]. For cavity detuning 0 < δ < 1/2, the system self-organizes and enters the superradiant phase above a critical pump rate η c , which is roughly consistent with the mean-field prediction A(η c ) = 1/ √ 1 -2δ [11]. We plot ρ(x) and ρ(k) at two representative detunings in Fig. 2(a)-(d). The self-organized superfluid (SSF) phase is characterized by a continuous density distribution ρ(x) with pronounced peaks at the sites of the emergent lattice with spacing λ c = 2π/k c . The corresponding ρ(k), which can be measured by time-of-flight experiments, is characterized by a principal peak at the center k = 0 straddled by two satellite peaks at k = ±k c stemming from the superfluid correlations between the atoms at different lattice sites [14,17,18]. At lower detunings δ < 0.2 and larger pump rate, the system transitions from the superfluid into a self-organized Mott insulator (SMI) phase. This phase is characterized by the disappearance of the peaks at k = ±k c and the broadening of the central peak at k = 0 in ρ(k) [Fig. 1(c)] [14,18,45,46]. The superfluid and Mott insulating phases are analogue to the ones in a standard Bose-Hubbard model [18].

As opposed to the red-detuned case, new correlated states of the atoms emerge as the detunings and pump rates are varied. As the system transitions from the single-well to the double-well optical lattice, the atoms collectively enter a new stage of self-organization. Depending on the degree of correlations between the atoms at different sites, we obtain either a self-organized dimerized superfluid (SDSF) where global superfluid correlations persist across the double-well dimers, or a self-organized second-order superfluid (2-SSF) phase where superfluid correlations exist only within each double-well dimer. The signatures of these two states lie in the distributions ρ(x) and ρ(k) and the correlation function g (1) (x, x ) as shown in Figs. 2(e)-(j). The double-well optical lattice is confirmed by the two-humped density distribution in ρ(x) at each lattice site [Fig. 2(e,h)], and the concomitant reduction of the one-body correlation from unity within one lattice site [Fig. 2(g,j)]. Within each double-well dimer, local superfluidity exists and manifests itself as two peaks in ρ(k) appearing at k = ±k * [Fig. 2(f,i)], where

k * = πk c arccos[(B -δ)/|θ|] (4)
corresponds to the distance between the minima of the double-well potential. As the pump rate increases, k * approaches 2k c from above and the peak height increases as the double well becomes deeper.

In the SDSF (2-SSF) phase, global superfluid correlations between different pairs of double wells is present (absent). This corresponds to the presence (absence) of the peaks at k c in ρ(k) [Fig. 2(f,i)], and a finite (vanishing) correlation in g (1) between different lattice sites [Fig. 2(g,j)]. In a 2-SSF state, the correlation length of superfluidity has a completely different length scale from the SSF and SDSF states, since the correlations exist only locally within each double-well dimer. These two new phases realize a variant of the Bose-Hubbard model with degenerate double-well lattices with Hamiltonian

ĤBH = - i (t 1 ĉ † i,L ĉi,R + t 2 ĉ † i,R ĉi+1,L + H.c.) + i,σ=L,R U 2 (ĉ † i,σ ĉi,σ ) 2 + µ i ĉ † i,σ ĉi,σ , (5)
where L, R denote the subsites [47]. We find that the single-well optical lattice smoothly deforms to the double-well lattice and we have not been able to numerically establish if the system transitions or crossovers to the SDSF and 2-SSF phases [29]. Nonetheless, as shown in detail in Ref. [29], clear hysteretic behavior is seen across the boundary between SDSF and 2-SSF phases [see the hatched dark green region in Fig. 1(b)]. This implies that the transition between the globally superfluid phase and the globally Mott insulating phase is of first order. At higher pump rates, dynamical instabilities appear which preclude the superfluid-Mott transition within a double-well dimer. These were first discussed in Ref. [30], where in the absence of interactions and a trap, the limit cycle phase was predicted at the mean-field level. The instabilities result from a Hopf bifurcation [START_REF] Strogatz | Nonlinear Dynamics and Chaos[END_REF] (see Ref. [29] for detail) and manifest themselves as strange attractors in the region denoted by the thick solid gray line and as chaos in the pink region in Fig. 1(b). In Fig. 3(a), we compare the sensitivity of the strange attractor to three different ramping protocols: the trajectory always converges to the same strange attractor in the (B(t), θ(t)) phase space with roughly the same amplitude and frequency profile (not strictly periodic in time) in all three cases [Fig. 3(b)]. The magnitude of the oscillations increases with the pump rate, eventually leading to temporary vanishing of both θ and the optical lattice. Consequently, higher momentum modes are excited and the system becomes thermalized. We emphasize that in the presence of atomic interactions and the harmonic trap, strange attractors are the norm rather than limit cycles predicted in Refs. [30,31], and tightly trapped systems are more prone to thermalization (see Ref.

[29] for detail).

Extension to 2D systems -Our phase diagram for the 1D system can be straightforwardly extended to the experimentally relevant 2D systems. In this case, the atoms are subject to three light-induced effective potentials which we require to be:

V light (x, y) = η 2 /U 0 cos 2 (k c y) + U 0 |α| 2 cos 2 (k c x) + η(α + α * ) cos(k c x) sin(k c y). (6)
Compared to the 1D version [Eq. (1a)], the first term is an additional term stemming from the standing wave induced by the transverse pumping laser. The cavity field α follows the equation of motion Eq. (1b), where the lattice order parameter is generalized to θ = dxdy Ψ † (x, y) Ψ(x, y) cos(k c x) sin(k c y)/N . Compared to the effective potential realized very recently in experiment [33], the essential difference in the above system [Eq. ( 6)] is an extra phase shift of π/2 in the interference (last) term along the pump (y) axis. This phase shift is necessary to observe our phase diagram and can straightforwardly be implemented in the experimental setup.

The 2D system is expected to qualitatively possess the same phase diagram shown in Fig. 1(b). All the phases should be accessible for reasonable values of the detuning and pump rate. In Figs. 3(c,d), we plot two representative examples of spatial density distributions based on mean-field simulations, i) the standard self-organization on a checkerboard lattice corresponding to an SSF or SMI, and ii) self-organization on the double-well optical lattice structures. The latter would lead to extra peaks in momentum space at k = (±k * , 0) superposed over the underlying Mott insulating or superfluid momentum distribution, which serves as the smoking-gun evidence of the SDSF and 2-SSF phases in experiment. Note that oscillatory phases, reminiscent of our strange attractors have already been observed experimentally [33]. Despite the phase shift between the experimental setup and the one discussed in this work, we expect the mechanism of the dynamical instabilities to be the same. To clearly observe the strange attractor in an experiment, a loose trap along the cavity axis with ω x ∼ 10 -3 ω R is suggested.

Our work illustrates the potential of blue detuning to realize exotic phases of matter and should stimulate future studies of complex cavity-cold atom platforms, both bosonic and fermionic. harmonic trap confining the atoms has a trapping frequency of ω x = 0.136ω R = 2π ×504Hz in most of the simulations, but looser traps with ω x = 0.068ω R = 2π × 252Hz and ω x = 0.0068ω R = 2π × 25.2Hz are also used to analyze the strange attractor phase.

III. CROSSOVER OR TRANSITION TO SDSF AND 2-SSF PHASES

To investigate if the spontaneous deformation to the double-well lattice signifies a transition or a crossover, we choose the height of the peak at k = k * in momentum space as the relevant order parameter,

ξ = 0, B -δ ≥ |θ| ρ(k * ), B -δ < |θ|, (S5)
where

k * = πk c arccos[(B -δ)/|θ|] . (S6)
This ad hoc order parameter ξ as a function of A at two different detunings δ = 0.14 and δ = 0.28 is shown in Fig. S1. As discussed in the main text, for δ = 0.14 the system enters the normal phase (NP), the self-organized superfluid (SSF) phase, the self-organized Mott insulator (SMI) phase and the self-organized second-order superfluid phase (2-SSF) phase sequentially. In the latter case, the system enters the NP, SSF phase and self-organized dimerized superfluid (SDSF) phase sequentially.

In both cases, ξ gains a tiny finite value as soon as θ > Bδ, as the effective potential [Eq. (3) in the main text] at each lattice site evolves from a single well into a double well. However, at δ = 0.14, ξ increases slowly at first in the SSF and SMI phases and then it increases rapidly after entering the 2-SSF phase. The onset of the rapid increase (A ≈ 9) roughly traces the boundary between the SMI and 2-SSF phases. On the contrary, at δ = 0.28, the increase in ξ is almost negligible. It is more likely that the system undergoes a crossover between the SSF and SDSF phases and there is no clear boundary between these two phases.

Tracing through the green line indicating the dimerization transition in the phase diagram [Fig. 1(b) in the main text], we have seen two different transition behaviors depending on the existence of the global superfluidity. When the detuning is large and the system is globally superfluid, a crossover is likely to take place; while when the detuning is small and the system is globally Mott insulating, a second order transition is likely to take place. A firm conclusion can be drawn only after more detailed examination, including the dependence of the transition behaviors on atom number and orbital number. FIG. S2: The hysteretic behaviors across the SDSF -2SSF boundary. The cavity detuning is ramped up (red lines) and down (blue lines) between δ = 0.211 (∆c = 45) and δ = 0.234 (∆c = 50) at different pump rates (a) A = 2.9, ramping within SSF, (b) A = 4.4, and (c) A = 5.9. The ramping is performed at different rates, where the fastest ramp is performed within a time interval of t = 13.6/ωR (solid lines), the medium ramp within t = 40.7/ωR (dashed lines), and the slowest ramp within t = 67.8/ωR (dotted lines). The dimensionless cavity detuning δ is given with respect to N U0.

IV. HYSTERESIS IN ATOMIC CORRELATIONS BETWEEN THE TWO DIMERIZED PHASES

We now analyze the hatched dark green region separating the self-organized dimerized superfluid (SDSF) phase and the self-organized second-order superfluid (2-SSF) phase in Fig. 1(b) in the main text. A systematic analysis shows that the correlations of the system hardly converge close to the boundary between the SDSF and 2-SSF phases. The height of the central peak in the momentum space density distribution ρ(k = 0) has been shown in the literature to serve as a simple and useful indicator of the correlations between the atoms [9][10][11]. In Fig. S2, we fix the pump rate at three different values and ramp the cavity detuning back and forth across the SDSF -2SSF boundary at different ramping rates. At small pump rate A = 2.9 [Fig. S2(a)], the system is still in the SSF phase and no hysteresis exists even for the fastest ramp. At higher pump rate A = 4.4 [Fig. S2(b)], the system enters the SDSF phase and the 2SSF phase at high and low detunings, respectively. A hysteresis in ρ(k = 0) can be clearly seen in a fast ramping. However, the hysteresis area becomes much smaller as the ramping slows down. At an even higher pump rate A = 5.9 [Fig. S2(c)] where the system is approaching the thermalized phase, the hysteresis is prominent, and the hysteresis area does not vanish even for the slowest ramping. Our results indicate that the SDSF-2-SSF transition is a first-order transition.

V. SENSITIVITY OF THE STRANGE ATTRACTOR TO THE HARMONIC TRAP AND THE ATOMIC FLUCTUATIONS

The presence of trapping and atomic interactions have large impact on the transition to chaos. The strange attractor is highly susceptible to the harmonic trap and atomic fluctuations. To quantify the dynamical instability in the system, we propagate the system to the desired pump rate using different ramp-up times. We then compute the magnitude of the system's oscillation through the variance of θ, i.e.,

∆θ = θ 2 -θ 2 , ( S7
)
where the bar means average in time. For comparison, these propagations are performed with different values of the harmonic trapping frequency, both in the mean-field limit [Fig. S3(a)] and beyond mean-field [Fig. S3(b)]. We find that a looser trap shifts the boundary between the stable and unstable region towards higher pump rates and also widens the SA region. In the mean-field limit, as the trap loosens, the width of the SA region (∆A) becomes wider and approaches the limit ∆A ≈ 0.8 predicted by Ref. [12] in the absence of a trap. The SA region shrinks even further as we turn on the atomic fluctuations. Particularly, with the tight trapping, ω x /ω R = 0.136, the SA region between the SDSF phase and the thermalized region vanishes completely. In conclusion, the strange attractors are extremely susceptible to the atomic fluctuations and the harmonic trap frequency. A tight trap can drive the strange attractors into chaos. In experiments, a loose trap with ω x ∼ 10 -3 ω R is recommended for the observation of the strange attractors.

VI. EVOLUTION OF THE PARAMETERS θ AND B AND THE FIXED POINTS OF THE SYSTEM

The mean-field dynamical behaviors of the blue-detuned system have already been investigated in Ref. [12], and Ref. [13] has used free energy arguments to show that steady state no longer exists at large pump rate in the system. However, the mechanism of the dynamical instabilities was unclear. In this section, we will map the system into discretized time and show that the dynamical instabilities originate from Neimark-Sacker bifurcation [14,15], which is the equivalence of the Hopf bifurcation [16,17] in a continuous system. The dynamical behavior of the system can already be captured by the evolution of the parameters θ and B. These two parameters are chosen because they enter the effective Hamiltonian [cf. Eq.(3) in the main text] and they quantify the self-organization. We derive the discrete dynamical map governing the evolution of these two parameters and show how the system evolves from a stable steady state to a limit cycle and finally to chaos with increasing pump rate.

In the bad-cavity limit, the cavity is always in its steady state and the atomic field evolves according to the following mean-field steady-state Hamiltonian Eq.(3) in the main text:

Ĥ(1) = - ∂ 2 x 2m + A 2 ω R [2(δ -B)θ cos(k c x) + θ 2 cos 2 (k c x)]. (S8)
Here, the dimensionless overall effective potential strength A and the dimensionless cavity detuning δ are related to the original system parameters through

A = ηN √ U 0 / [(∆ c -N U 0 B) 2 + κ 2
]ω R > 0 and δ = ∆ c /N U 0 , while θ and B are given by the instantaneous atomic density θ = dxρ(x) cos(k c x) and B = dxρ(x) cos 2 (k c x), with ρ(x) = Ψ † (x) Ψ(x) /N the density distribution. We note that the self-consistency of the sign of θ requires that the coefficient δ -B should be negative. For simplicity, we have ignored the atom-atom interaction and the harmonic trap, and treat the system in the mean-field limit in this section.

Due to the driven-dissipative nature of the system, the atoms are subject to an effective potential which is determined by the atomic distribution itself through θ and B. On the other hand, once the instantaneous effective potential is known, we are able to solve the ground state of the Hamiltonian, which will give a new pair of parameters θ and B. In this way, we are able to obtain discretized evolution equations for these two parameters. For this purpose, we investigate the effective potential, which is given by the second and third terms of Eq. (S8).

When Bδ ≥ |θ| > 0, the minima of the effective potential lie at x n = 2nπ/k c [x n = (2n + 1)π/k c ], n ∈ N, for positive (negative) θ. In the vicinity of the minima, the effective potential can be expanded into a quadratic form In all these panels, different values of A/ωR are chosen, and δ is always chosen as δ = 0.188. (a,f) For A = 0.5ωR, there is a stable fixed point at (θ, B) = (0, 1/2), and the cavity-BEC system is in the normal phase; (b,g) for A = 2ωR, there is a stable fixed point with 0 < θ < Bδ ≈ 0.5δ, and the cavity-BEC system self-organizes into a lattice without lattice site dimerization; (c,h) for A = 10ωR, there is a stable fixed point at 0 < Bδ < θ, and the lattice sites of the cavity-BEC system dimerize; (d,i) for A = 50ωR, the fixed point becomes unstable, and the trajectory converges to a limit cycle surrounding the fixed point; (e,j) for A = 90ωR, the fixed point is unstable, and θ switches sign after a few steps. In the cavity-BEC system, this indicates that the effective potential vanishes and the system will become thermalized easily. The dimensionless detuning δ and the potential strength A are normalized with respect to N U0 and √ ωR, respectively.

V (x = x n + δx) = ω R A 2 |θ|(B -δ -|θ|)k 2 c δx 2 ≡ m 2 Ω 2 1 δx 2 . ( S9
In the absence of atom-atom interaction, the ground state of the effective potential can be approximated by a sum of Gaussians,

ρ(x) ∝ n exp - mΩ 1 (x -x n ) 2 . ( S10
)
Suppose that at time t the density distribution of the atoms is given by ρ(x; t), which is further determined by the instantaneous θ t and B t .. At the next moment, the parameters θ and B are evolving to θ t+1 = ∞ -∞ cos(k c x)ρ(x; t) and B t+1 = ∞ -∞ cos 2 (k c x)ρ(x; t), respectively. According to the aforementioned ansatz of the density distribution, these are given explicitly by

θ t+1 = sgn(θ t ) exp - k 2 c 4mΩ 1 = sgn(θ t ) exp - 1 4 A 2 |θ t |(B t -δ -|θ t |)/2 (S11a)
B t+1 = 1 2 + 1 2 exp - k 2 c mΩ 1 = 1 2 + 1 2 exp - 1
A 2 |θ t |(B t -δ -|θ t |)/2 (S11b)
The fixed point of the system is given by (θ t+1 , B t+1 ) = (θ t , B t ), and (θ, B) = (0, 1/2) is always a fixed point for any value of A. However, a second fixed point might exist depending on the value of A. When A is large enough, the following equation, For small values of A, it is smaller than one |λ1| < 1 and the fixed points are stable (solid dots). It becomes greater than one |λ1| > 1 for A 2 > 45 and the fixed points become unstable (empty dots). The system undergoes a Neimark-Sacker bifurcation at roughly A 2 = 45. (b) The trajectory of the limit cycles on the (Bδ, θ) phase space with three different values of A 2 = 50, 60, 90. For the smallest value of A 2 = 50 (red dots), the system oscillates among many points on the phase space. For A 2 = 60 (blue dashed line), the system oscillates among only 5 points. For A 2 = 90 (orange dotted line), θ switches sign in every few steps, which eventually causes thermalization in the cavity-BEC system. (c) The fixed points on the (Bδ, θ) phase space as a function of A. The solid black dots are stable fixed points, while the empty black dots are unstable fixed points. A blue line is superposed to show the trajectory in MCTDH-X simulation, where the same detuning δ = 0.188 is used and the pump rate A ∝ η is ramped up linearly in time. The results for 2.5 < A 2 < 5 are not shown, because in these cases the fixed point is close to the line θ = Bδ, and the results become diverging and unreliable. In panels (a-c), the detuning takes the value δ = 0.188. (d) The mean-field phase diagram obtained by (solid lines) the discretized time analysis Eq. (S11) and (S14), and (dashed lines) MCTDH-X simulations. The simulation results are reproduced from Fig. 1(b) of the main text. The blue lines show the boundary between the normal phase (NP) and the superradiant phase (SP), while the red lines separate the dynamically stable and unstable regions. For small detunings δ < 0.3, the system behaviors are qualitatively the same between the discretized time analysis and the simulations. For large detunings δ > 0.3, the discretized time analysis predicts that the system will become dynamically unstable as soon as it enters the SP. The dimensionless detuning δ and the potential strength A are normalized with respect to N U0 and √ ωR, respectively.

8θ ln 2 θ 1 2 -δ + 1 2 θ 4 -θ = 1 A 2 , ( S12
When |θ| > Bδ, the minima x n of the effective potential satisfy cos(k c x n ) = (Bδ)/θ. In the vicinity of the minima, the effective potential can be expanded as

V (x = x n + δx) = ω R A 2 k 2 [θ 2 -(δ -B) 2 ]k 2 c δx 2 ≡ m 2 Ω 2 2 δx 2 . (S13)
With the Gaussian ansatz of the ground state, the evolution of the parameters θ and B are given by

θ t+1 = exp - k 2 c 4mΩ 2 cos(k c x n ) = exp - 1 4 A 2 [θ 2 t -(δ -B t ) 2 ]/2 B t -δ θ t (S14a) B t+1 = 1 2 + 1 2 exp - k 2 c mΩ 2 (2 cos 2 k c x n -1) = 1 2 + 1 2 exp - 1 A 2 [θ 2 t -(δ -B t ) 2 ]/2 2 B t -δ θ t 2 -1 (S14b)
Fig. S4(a-e) shows the flow of the parameter pair (θ, Bδ) according to Eqs. (S11) and (S14), with δ = 0.188 and different values of A. For a small potential strength A, the fixed point lies at (θ, Bδ) = (0, 1/2δ) corresponding to the normal phase (NP). As the potential strength increases, the fixed point moves to a finite θ with B remaining roughly 1/2. The system has entered the superradiant phase (SP) and is in either the self-organized superfluid (SSF) phase or the self-organized Mott insulator (SMI) phase, since Bδ > |θ| and the minima of the potential lie at x n = 2nπ/k c . As A increases even further, the fixed point moves to the lower right part of the phase space with B -δ < |θ|. As discussed in Sec. III, in this region the effective potential becomes a double well at each lattice site, but the system has not necessarily already entered the self-organized dimerized superfluid (SDSF) or the self-organized second-order superfluid (2-SSF) phases. The SDSF and 2-SSF phases come from correlation effects which cannot captured by the mean-field model discussed in this section. Exemplary trajectories of the system with different values of A and a randomly chosen initial point are shown in Fig. S4(f-h).

The flows for A 2 = 10, 50 and 90 look similar except that the fixed point gradually moves towards smaller values of θ and B. However, the trajectories of the system with these three different values of A show completely different behaviors [Fig. S4(h-j)]. The trajectory converges to a stable fixed point when A 2 = 10, but to a limit cycle when A 2 = 50. At A 2 = 90, the order parameter θ flips sign after several steps of evolution. In the cavity-BEC system, this means that the effective potential vanishes at some point. As a result, higher momentum modes will be excited and the system will become thermalized and chaotic.

To investigate the change in the dynamical behavior of the system, we linearize the system in the vicinity of the fixed point and calculate the Jacobian

J =    ∂θ t+1 ∂θ t ∂θ t+1 ∂B t ∂B t+1 ∂θ t ∂B t+1 ∂ B t    . (S15)
We denote the eigenvalues of the Jacobian as λ 1 and λ 2 . For small A such that the fixed point lies in the upper left part of the phase space (Bδ > θ), the Jacobian has a non-zero 0 < |λ 1 | < 1 and a zero λ 2 = 0 eigenvalues. For large A such that the fixed point lies in the lower right part of the phase space (Bδ < θ), we have two complex eigenvalues which are conjugates of each other λ 1 = λ * 2 . The modulus of the eigenvalues |λ 1 | = |λ 2 | as a function of A is shown in Fig. S5(a). As A increases, the modulus also increases gradually while the real parts of λ 1 and λ 2 remain positive (not shown in the figure). After reaching the threshold A 2 ≈ 37, the modulus becomes greater than one. The system undergoes a Neimark-Sacker bifurcation [14,15]. A Neimark-Sacker bifurcation is the counterpart of the Hopf bifurcation in discretized time [16,17]. Its prerequisite is that the Jacobian for the fixed point (S15) has a pair of conjugate complex eigenvalues; it occurs when the moduli of these eigenvalues become larger than one; and it turns the stable fixed point into an unstable one with a limit cycle bifurcating from it. As A increases further, the area covered by the limit cycle continues to increase. The limit cycles for different values of A are shown in Fig. S5(b). Unlike in a period-doubling bifurcation, at a smaller A 2 = 40, this system oscillates among extremely many parameter points, but then at larger A 2 = 60 and A 2 = 90, it oscillates among much fewer parameters points. After exceeding A 2 = 82, the limit cycle expands to the phase space with negative θ. As discussed previously, this implies that the cavity-BEC system becomes thermalized and is driven into chaos.

Finally, we compare the discretized system [Eq. (S11), (S14)] and the simulation results of the cavity-BEC system by MCTDH-X. In Fig. S5(c), we compare the fixed point of the system obtained from the discretized system and the trajectory of the system in MCTDH-X simulations, which agree with each other qualitatively. In Fig. S5(d), we compare the normal phase -superradiant phase boundary and the stable -unstable boundary obtained from both methods, and they also agree qualitatively.

In conclusion, we have mapped the system into the (θ, B) phase space, and investigated its dynamics in discretized time. The mapping reveals the dynamical behaviors of the system qualitatively. It can not only capture the mean-field transition from the normal phase via the undimerized superradiant phase to the dimerized superradiant phase, but also show that the dynamical instabilities, including limit cycles (strange attractors in the trapped system) and chaos, arise from Hopf bifurcation in the original cavity-BEC system.

FIG. 1 :

 1 FIG.1:(a) Schematic setup of a trapped cavity-BEC. (b) Phase diagram. For cavity detuning 0 < δ < 1/2, the system transitions from the normal phase to the dynamically unstable region via superradiance with increasing pump rate A ∝ η. The strongly correlated phases are self-organized superradiant (SSF), a self-organized Mott insulator (SMI), a self-organized dimerized superfluid (SDSF), and a self-organized secondorder superfluid (2-SSF) phase. The orange line delineates superfluid phases and globally Mott insulating phases, while the green line marks the deformation to a double-well potential. Pronounced sensitivity to the ramping protocol is seen in the hatched dark green region. At higher η, the system is dynamically unstable to the formation of strange attractors followed by chaos. In presence of a trap and atomic fluctuations, the strange attractors only exist in the region represented by the thick gray line, while the thin dashed section represents a direct transition to the chaos. The dimensionless detuning δ and the potential strength A are normalized with respect to N U0 and √ ωR, respectively. (c) Sketch of the SSF,

1 FIG. 2 :

 12 FIG.2: (a-d) The momentum and position space density distributions ρ(k) and ρ(x) as a function of pump rate A ∝ η at two detunings δ = 0.14 and δ = 0.28. At lower detuning, δ = 0.14, the system starts from the normal phase and then enters the SSF phase at A = 1.4, the SMI phase at A = 2.7, and the 2-SSF phase sequentially. At higher detuning δ = 0.28, the system starts from the normal phase and then enters the SSF phase at A = 1.9, and the SDSF phase sequentially. The dotted lines are guides to the eye. The dimensionless potential strength A is normalized with respect to √ ωR. (e-j) The position and momentum space density distributions and the Glauber one-body correlation function of an SDSF state (first row) and a 2-SSF state (second row). In Figs.(g,j), the color code follows the functionln(1g(1) ). The double-well splitting is seen in the central lattice site but not in the other two lattice sites, because only M = 4 orbitals are used in the numerical simulations[29].

3 FIG. 3 :

 33 FIG. 3: (a)The order parameter θ as a function of time for three different ramping protocols from the stable region into the strange attractor phase (δ = 0.14, A = 11.8). In the blue and green trajectories, the detuning is fixed at δ = 0.14 and the pump rate A ∝ η is ramped up linearly at rates of (blue) dA/dt = 4 × 10 -4 ωR and (green) dA/dt = 8 × 10 -5 ωR.In the orange trajectory, the pump rate is fixed at A = 11.8 and the detuning is ramped up linearly at a rate of dδ/dt = 1.5 × 10 -4 ωR. (b) For all three cases, the system converges to the same strange attractor in the (Bθ) phase space. (c,d) The position space density distribution ρ(x, y) of (c) a single-well state and (d) a double-well lattice state in a twodimensional cavity-BEC system as described by Eq. (6).

28 FIG

 28 FIG. S1:The dimerization order parameter ξ [Eq. (S5)] and θ -B + δ as functions of pump rate at different detunings (a) δ = 0.14 and (b) δ = 0.28. The dimerization order parameter ξ is multiplied by a factor of N = 50 such that it is in the same order of magnitude as θ -B + δ. The dimensionless potential strength A is given with respect to √ ωR.

3 A

 3 FIG.S3: Measurement of the variance of the order parameter ∆θ [Eq. (S7)] over time as the system is propagated at different pump rate at δ = 0.4. The propagation is done with different trapping frequencies (ωx/ωR = 0.136, 0.068, and 0.0068) both (a) in the mean-field limit and (b) beyond. The crosses mark the largest value of pump rate where the system is not thermalized. The dimensionless potential strength A is given with respect to √ ωR.

)

  FIG.S4: (a-e) The flow (θt+1θt, Bt+1 -Bt) on the (Bδ, θ) phase space according to Eq. (S11) for Bδ > |θ| and Eq. (S14) for Bδ < |θ|. The arrows have unit length and only show the directions. (f-i) The trajectory of the parameter pair (θ, B) on the phase space with randomly chosen initial values θt=0 = 0.3, Bt=0δ = 0.1. The initial points are marked by crosses in the figures. In all these panels, different values of A/ωR are chosen, and δ is always chosen as δ = 0.188. (a,f) For A = 0.5ωR, there is a stable fixed point at (θ, B) = (0, 1/2), and the cavity-BEC system is in the normal phase; (b,g) for A = 2ωR, there is a stable fixed point with 0 < θ < Bδ ≈ 0.5δ, and the cavity-BEC system self-organizes into a lattice without lattice site dimerization; (c,h) for A = 10ωR, there is a stable fixed point at 0 < Bδ < θ, and the lattice sites of the cavity-BEC system dimerize; (d,i) for A = 50ωR, the fixed point becomes unstable, and the trajectory converges to a limit cycle surrounding the fixed point; (e,j) for A = 90ωR, the fixed point is unstable, and θ switches sign after a few steps. In the cavity-BEC system, this indicates that the effective potential vanishes and the system will become thermalized easily. The dimensionless detuning δ and the potential strength A are normalized with respect to N U0 and √ ωR, respectively.

  ) has a root in the interval θ ∈ [0, 1]. The normal-superradiant phase boundary predicted in this manner is compared to the numerical results from MCTDH-X in Fig.S5(d).

  FIG. S5: (a)The modulus of the eigenvalue of the Jacobian |λ1| as a function of A, when the fixed point lie on the lower right part of the phase space (A 2 > 5, θ > Bδ). For small values of A, it is smaller than one |λ1| < 1 and the fixed points are stable (solid dots). It becomes greater than one |λ1| > 1 for A 2 > 45 and the fixed points become unstable (empty dots). The system undergoes a Neimark-Sacker bifurcation at roughly A 2 = 45. (b) The trajectory of the limit cycles on the (Bδ, θ) phase space with three different values of A 2 = 50, 60, 90. For the smallest value of A 2 = 50 (red dots), the system oscillates among many points on the phase space. For A 2 = 60 (blue dashed line), the system oscillates among only 5 points. For A 2 = 90 (orange dotted line), θ switches sign in every few steps, which eventually causes thermalization in the cavity-BEC system. (c) The fixed points on the (Bδ, θ) phase space as a function of A. The solid black dots are stable fixed points, while the empty black dots are unstable fixed points. A blue line is superposed to show the trajectory in MCTDH-X simulation, where the same detuning δ = 0.188 is used and the pump rate A ∝ η is ramped up linearly in time. The results for 2.5 < A 2 < 5 are not shown, because in these cases the fixed point is close to the line θ = Bδ, and the results become diverging and unreliable. In panels (a-c), the detuning takes the value δ = 0.188. (d) The mean-field phase diagram obtained by (solid lines) the discretized time analysis Eq. (S11) and (S14), and (dashed lines) MCTDH-X simulations. The simulation results are reproduced from Fig.1(b) of the main text. The blue lines show the boundary between the normal phase (NP) and the superradiant phase (SP), while the red lines separate the dynamically stable and unstable regions. For small detunings δ < 0.3, the system behaviors are qualitatively the same between the discretized time analysis and the simulations. For large detunings δ > 0.3, the discretized time analysis predicts that the system will become dynamically unstable as soon as it enters the SP. The dimensionless detuning δ and the potential strength A are normalized with respect to N U0 and √ ωR, respectively.

We acknowledge the financial support from the Swiss National Science Foundation (SNSF), the ETH Grants, Mr. Giulio Anderheggen, the Austrian Science Foundation (FWF) under grant P32033. We also acknowledge the computation time on the ETH Euler cluster.

Supplementary Material for Pathway to chaos through hierarchical superfluidity in a cavity-BEC system Rui Lin, 1 Paolo Molignini, 1 Axel U. J. Lode, 2 and R. Chitra In this work, the Multiconfigurational Time-Dependent Hartree method for indistinguishable particles [1][2][3] (encoded in the software MCTDH-X [4]) is used to evolve the system's state. This algorithm is able to solve a Hamiltonian with a one-body potential V (x) and a two-body interaction W (x, x ),

MCTDH-X is based on the following ansatz for the many-body wave function

where N is the number of atoms, M is the number of single-particle wave functions (orbitals) and n = (n 1 , n 2 , ..., n k ) gives the number of atoms in each orbital with the constraint

where ψ i (x; t) are the wave functions of the individual working orbitals. The state is evolved by propagating the coefficients C n (t) and the working orbitals ψ i (x; t) using the time-dependent variational principle [5] In MCTDH-X, the number of orbitals M is important to obtain converged results, especially of correlation functions [2,3,6]. With a single orbital M = 1, the method reproduces the Gross-Pitaevskii mean-field limit; while with infinitely many orbitals M → ∞, the method is numerically exact. To simulate a Mott insulating state, the number of required orbitals is usually the same as the number of lattice sites [7]. In our case, considering that there are three lattice sites and each site splits into two sub-sites, M = 6 orbitals are needed in total to fully capture the atomic correlations. However, due to the amount of computational time needed, the simulations are performed with M = 3 orbitals or M = 4 orbitals when specified. In the latter case, we are able to observe the effect of the correlation induced dimerization on the one-body correlation function inside the central lattice site, as shown in Fig. 2(b) in the main text. We expect that with M = 6 orbitals, reduction of correlation from unity can also be seen in the other two lattice sites.

II. SYSTEM PARAMETERS

The parameters used in this work are inspired by the experimental setup in the ETH Zürich experiments [8]. We simulate a total number of N = 50 87 Rb atoms with mass m = 1.44 × 10 -25 kg and contact interaction N g = 1.0 × 10 -17 eV • m. The wavelength of the laser pump and the cavity field is chosen as λ c = 784nm, and therefore the recoil frequency is ω R = 2π × 3734Hz. The collective photon light shift is given by N U 0 = 28.9ω R = 2π × 108kHz, and the cavity detuning takes the values between ∆ c = 0 and ∆ c = N U 0 /2 = 2π × 54kHz. The maximum pump rate is η = 2π × 825kHz, this corresponds to an overall effective potential strength of A = 12 [cf. Eq.(3) in the main text]. The cavity dissipation rate is chosen to be large enough κ = 700ω R = 2π × 2.60MHz such that the system operates in the bad-cavity limit. We note that the numerical simulations can also be done for less dissipative cavities. The arXiv:1910.01143v1 [cond-mat.quant-gas] 2 Oct 2019