Signatures of quantum chaos of Rydberg dressed bosons in a triple-well potential

We study signatures of quantum chaos in dynamics of Rydberg dressed bosonic atoms held in a one dimensional triple-well potential. Long-range nearest-neighbor and next-nearest-neighbor interactions, induced by laser dressing atoms to strongly interacting Rydberg states, affect drastically mean field and quantum many-body dynamics. By analyzing the mean field dynamics, classical chaos regions with positive and large Lyapunov exponents are identified as a function of the potential well tilting and dressed interactions. In the quantum regime, it is found that level statistics of the eigen-energies gains a Wigner-Dyson distribution when the Lyapunov exponents are large, giving rise to signatures of strong quantum chaos. We find that both the time averaged entanglement entropy and survival probability of the initial state have distinctively large values in the quantum chaos regime. We further show that population variances could be used as an indicator of the emergence of quantum chaos. This might provide a way to directly probe quantum chaotic dynamics through analyzing population dynamics in individual potential wells.


Introduction
Understanding dynamics of quantum many-body systems has been a lucrative field of study that allows for the exploration of new physics and finding quantum technological applications. Among many experimental platforms, ultracold atoms trapped in optical lattices, due to their high controllability, provide a versatile toolbox for studying quantum many-body phases [1]. One is typically keen to achieve long-range interactions, which allow us to create and probe exotic many-body dynamics. A paradigmatic example is the extended Bose-Hubbard model (EBHM), where rich phases (e.g., superfluid, supersolid and checkboard phases) are obtained [2][3][4][5][6]. Dipolar atoms provide relatively weak longrange interactions [7]. Recent experiments have shown that much stronger long-range interactions can be achieved in electronically high-lying Rydberg states [8]. Long coherence times are realized through Rydberg dressing, where ground-state atoms are weakly coupled to Rydberg states using far off-resonant lasers [9]. Such Rydberg dressing leads to a softcore-shaped long-range interaction potential between ground-state atoms, where the radius of the soft-core potential is in the order of a few micrometers [9][10][11][12][13][14][15][16].
Existing studies have focused on static and dynamical properties of Rydberg-dressed atoms confined in traps [17][18][19][20][21][22][23] and optical lattices [24][25][26][27][28][29][30][31]. Interaction effects due to the Rydberg dressing have been experimentally demonstrated in optical tweezers [32], optical lattices [33][34][35] and harmonic traps [36]. In optical lattice potentials, such a large soft-core radius is typically much larger than the lattice spacing, as depicted in Figure 1a. The dynamics of Rydberg-dressed bosons are described by an EBHM [14]. In Refs. [31,37], we studied nonlinear and chaotic dynamics of Rydberg-dressed bosons in finite potential wells in the semiclassical regime. In the quantum regime, it is known that the BHM is non-integrable in general, where quantum chaos can be triggered by the competition between coherent hopping and on-site interactions [38][39][40][41][42]. In the presence of long-range interactions, it has been shown that quantum chaos emerges in EBHMs realized with dipolar bosons, which feature nearest-neighbor interactions [43,44]. The respective chaotic properties are characterized by quantities such as level statistics [44], Shannon entropy, and survival probabilities of the initial state [43]. (a) Three-well trapping potential and soft-core-shaped interaction between Rydberg-dressed atoms. The soft-core radius R c is larger than the distance between neighboring sites. By tuning R C , both the nearest-neighbor and next-nearest-neighbor interaction can be made to be strong. In this work, we investigated signatures of the quantum chaos of Rydberg-dressed bosons in finite one-dimensional traps. Because of the large soft-core radius, our setup allowed us to focus on the dynamics in a region where both nearest-neighbor and nextnearest-neighbor interactions are equally strong. In the semiclassical regime, Lyapunov exponents of the mean-field equations were evaluated numerically. Positive and large-valued Lyapunov exponents were found, signifying the emergence of chaos. The dependence of the Lyapunov exponents on the dressed interaction and tilting of the potential was studied. By diagonalizing quantum many-body Hamiltonians with finite atom numbers, we show that nearest-neighbor level statistics gain a Wigner-Dyson distribution with parameters when the Lyapunov exponents are large. We further characterized signatures of quantum chaos through analyzing quantum many-body dynamics. It was found that the time-averaged entanglement entropy shows a similar dependence on the parameters to that of the level statistics and Lyapunov exponents. Features of quantum chaos can be captured by variances of the population in different potential wells, which might provide a way to directly probe the quantum chaos.
The structure of the paper is as follows. In Section 2, we introduce the system and model Hamiltonian. Section 3 is divided into three sub-sections. In Section 3.1, the level statistics of the system are studied. Parameters corresponding to Wigner-Dyson distribution are found. In Section 3.2, dynamics of the entanglement entropy are investigated. The time-averaged entanglement entropy gains larger values in the quantum chaos regime. In Section 3.3, we evaluate the survival probability of the initial state and particle numbers in different potential wells. Population variance in different regimes is examined. Conclusions are given in Section 4.

Model
We considered N bosonic atoms confined in a one-dimensional chain with L sites. As depicted in Figure 1a, the Rydberg-dressed soft-core potential induces interactions between atoms in different sites. The dynamics of this model are described by an extended Bose-Hubbard model [14] (h = 1): whereâ j andâ † j are bosonic annihilation and creation operators acting on site j, with the corresponding number operatorn j =â † jâ j . Here, ⟨i, j⟩ denotes the pair of nearest-neighbour indices. Parameter Γ j = −[j − 1 − ⌊L/2⌋]γ describes a local tilt of the potential well, in which γ represents the level bias between neighboring sites. J denotes the hopping rate of atoms between neighboring sites. The s-wave interaction is characterized by g = 4πa s /m, with a s and m being the s-wave scattering length and mass of the atom [45]. The long-range soft-core interaction induced by the Rydberg dressing is given by Λ i,j = C 6 /[|i − j| 6 d 6 + R 6 c ], where C 6 , d, and R c are an effective dispersion coefficient, lattice constant and soft-core radius, respectively.
Hilbert space of the EBHM is given by L = (N + L − 1)!/N!(L − 1)!. When N ∼ L, it increases rapidly when N is large. For example, L = 92378 when N = L = 10. Resulting numerical calculations become difficult in unit-filling situations. In this work, we focused on large filling regime (N ≫ L) in a minimal trapping setting, i.e., Rydberg-dressed bosons trapped in triple-well potential (L = 3). Here, the dimension of the Hilbert space was L = (N + 1)(N + 2)/2, which is quadratic with N. This allowed us to deal with large particle numbers N ∼ 100. Due to the high filling N/L ≫ 1, we could also directly compare the quantum many-body calculations with mean-field results.
With the above consideration, we obtained the Hamiltonian of the three-well system, where we defined W = Ng + 2NΛ j,j , U = NΛ j,j±1 and V = NΛ j,j±2 . We neglected constant terms that will not affect many-body dynamics. In following numerical analysis, we scaled the Hamiltonian with respect to J. For concreteness, we chose parameters such that the effective on-site interaction vanishes but nearest-neighbor (NN) interaction is strong, i.e., W = 0 and U = 2V [31]. The latter can be achieved by tuning the s-wave scattering length through Feshbach resonance [45]. This allows us to concentrate on effects due to the long-range interactions.

Results
Before discussing many-body calculations of the EBHM, we show the signature of chaos of this model in the semiclassical limit [31,37]. In the semiclassical (mean-field) calculation, one replaces the bosonic operators a i (a † i ) with complex classical fields ψ i (ψ * i ). The mean-field calculation is relatively simple, as the number of complex fields is equal to the number of site L. More details of the semiclassical calculations can be found in Refs. [31,37]. A key figure of merit in the mean-field calculation is Lyapunov exponents. Positive and large Lyapunov exponents indicate the emergence of classical chaos. This means that perturbations to initial states will grow exponentially with time. Examples of Lyapunov exponents are shown in Figure 1b,c. In Figure 1b, Lyapunov exponent λ is close to zero when U/J < 3.5 (fixing γ/J = 2.5), where dynamics are linear. Positive λ is found at around U/J ≥ 3.5, and it reaches its maximum at around U/J = 7.0. Focusing on this strong interaction region (U/J = 7.0), we calculated Lyapunov exponents by varying the tilting γ. We found that λ increases with an increasing γ, and achieves its maximum at around γ/J = 2.5.
We now turn to analyzing static and dynamics properties of the three-well system with finite N. By numerically diagonalizing Hamiltonian (2) eigenstates, |j⟩ ≡ ∑ n C j n |n⟩ and eigenvalue E j were obtained. Here, C j n is the probability amplitude of Fock state |n⟩ ≡ |n 1 , n 2 , n 3 ⟩ in which n i denotes the number of particles at the i-th site. Statistical distributions of E j provide valuable information on whether the system is integrable or chaotic [46]. If the system is integrable, corresponding energy levels are not correlated, and are not prohibited from direct crossing when varying parameters [47][48][49], whereas, in the quantum chaotic region, energy levels are correlated and crossings are avoided [46,50]. In the weak interaction limit U → 0.0 (i.e., Figure 2a), atom hopping dominates dynamics such that the system is integrable [51]. As we increase U/J to 7.0, we can see that anticrossings start to appear, which indicate the emergence of quantum chaos (Figure 2b). When fixing U/J = 7.0, we show eigen-energies in the interval 0.0 < γ < 1.0 in Figure 3a and 2.0 < γ < 3.0 in Figure 3b. The difference is that both direct and avoided crossings are found in Figure 3a, whereas only avoided crossings are encountered when γ is large, as shown in Figure 3b.  γ/J is fixed to be 0, 2.5 and 7, respectively. The dashed blue and solid red lines in (d-f) represent the Poissonian and WD distribution, respectively. When γ/J = 2.5, the levels approach the WD distribution. Other parameters are U/J = 7 and N = 120 in (d-f).

Level Statistics
To understand the different behaviors of the eigen-energies, we analyzed statistics of the nearest-neighbor spacing s between energy levels. The probability distribution P(s) of the spacing s can be used to distinguish integrable and chaotic dynamics [46]. For integrable systems, P(s) follows the Poissonian distribution, For chaotic systems, their level spacing corresponds to the Gaussian orthogonal ensemble (GOE) [52]. Distribution P(s) changes to the Wigner-Dyson (WD) distribution [53], In realistic situations, the level spacing will typically not follow the Poissonian or WD distribution exactly. P(s) can be fitted by the Brody distribution [54], where β is a chaos indicator and b depends on β, where Γ(x) is the gamma function. Here, 0 ≤ β ≤ 1 measures the repulsion of the levels, also known as the level repulsion exponent. For chaotic systems, β ∼ 1, where the Brody distribution recovers the WD profile, whereas β ∼ 0 in integrable systems, where the Brody distribution becomes the Poissonian distribution. After carrying out the standard unfolding procedure [46], the level statistics were calculated and fitted with Equation (5). Values of β were extracted and are shown in Figure 2c and Figure 3c, respectively. Fixing γ/J = 2.5, it was found that β is small when U/J → 0.0. In this region, the level statistics follow the Poissonian distribution, as illustrated in Figure 2d. Then, β keeps increasing with am increasing U/J until U/J ∼ 7.0, where the level statistics display the WD distribution i.e., Figure 2e. Hence, the system enters the quantum chaos regime. Further increasing U, β decreases gradually, where the statistics follow the Poissonian distribution again, i.e., Figure 2f. The results show that our model is largely integrable in the weakly and strongly interacting regimes, which is similar to the Bose-Hubbard model [40]. The trend of β shown in Figure 2c is consistent with the Lyapunov exponent shown in Figure 1b.
In Figure 3c, we plot β as a function of γ/J. Similar to the Lyapunov exponents shown in Figure 1c

Entanglement Entropy
In this section, we analyze the entanglement entropy of our model. The three wells were divided into subsystems of the left well (part A) and the other two wells (part B). Entanglement entropy with respect to subsystem A is given by [55] S EE = −Tr(ρ A ln ρ A ), where ρ A = Tr B (ρ) is the reduced density matrix by tracing out the subsystem B from the density matrix ρ = |ψ⟩⟨ψ| of the whole system. Entanglement entropy S EE measures how much subsystem A and B are intertwined with each other. For the three-well bosonic system, it has a maximal value S max EE = ln(N + 1), which only depends on the total number of atoms of the whole system (see Appendix A for calculations).
To evaluate the entanglement entropy numerically, we started with an initial state |ψ 0 ⟩ = ∑ n a n |n⟩. The time-dependent many-body state |ψ(t)⟩ = e −iĤt |ψ 0 ⟩ was expanded in the Fock basis, in which we defined α j ≡ ∑ n (C j n ) * a n . Given the time-dependent many-body state in the Fock basis, the entanglement entropy S EE (t) can be evaluated.
In Figure 5a, the dynamical evolution of S EE (t) is shown for U/J = 3.0, 7.0 and 13.0, with initial state |N/3, N/3, N/3⟩ and N = 120. The entropy initially increases monotonically for a short period of time and then fluctuates around a constant value with small amplitudes. Here, the saturation values of S EE (t) depend on U/J. If one increases U, amplitudes of the fluctuation will be further reduced. To explore the dependence of the entropy on parameters, we numerically calculated time-averaged and normalized entanglement entropy,S EE = 1/(S max EE ∆t) The normalization could partially eliminate the dependence ofS EE on N. In the numerical calculation, we chose t i = 4/J to avoid the influence due to the initial stage of the entropy. We set t f = 20/J in calculations, as S EE saturates when t f ≤ 20/J. It is important to note that long-time behaviors of S EE are independent of the initial configurations, though details of the dynamics depend on the configuration. By focusing on |N/3, N/3, N/3⟩ initially, we can compare dynamics of various parameters. This is particularly useful when γ = 0. The results ofS EE as a function of U are shown in Figure 5b, where γ/J was fixed to be 2.5. It is small when U/J ≪ 7.0.S EE then increases with U and gains its peak value at U/J ≈ 7.0. The peak value is close to the upper bound S max EE when N is large (e.g., N > 120). Further increasing U,S EE decreases gradually. This leads to a similar pattern to β, i.e., Figure 2c. Importantly, we note that the maximalS EE appears at the same U/J as that of β. Hence, both predict the position of quantum chaos as indicated by level statistics, i.e., Figure 2. We note thatS EE is shifted upwards (increased) globally when N is increased. Such a dependence is different from β, where increasing N will suppress the fluctuation and lead to a better fitting.
In Figure 5c,d, we show the dependence of the entanglement entropy on γ while fixing U/J = 7.0. The entanglement entropy first increases rapidly with time and then saturates when Jt ≳ 4.0, as depicted in Figure 5c. The saturation varies with γ. To explore this dependence, we evaluated the normalized average entropyS EE , and show the results in Figure 5d.S EE first increases with γ, reaches peak values around γ/J = 2.5 and then decreases with an increasing γ. Such a pattern is similar to the dependence of β on γ, as seen in Figure 3c. Increasing atom number N, the profile ofS EE is largely unchanged, but is shifted globally upwards.

Survival Probability and Variance of Populations
The survival probability P s (t) = |⟨ψ 0 |ψ(t)⟩| 2 gives the probability of finding initial state |ψ 0 ⟩ at time t. In the long time limit, P s (t) will stay in close proximity to that of GOE matrices if the dynamics are chaotic. The survival probability of GOEs saturates at P GOE s ≈ 3/D [56] at the long time limit, which depends only on the dimension D of the Hilbert space. This provides a way to characterize chaotic quantum dynamics. Considering initial state |n 0 ⟩ = |N/3, N/3, N/3⟩ (N = 90), the respective survival probability is where C j n 0 denotes the probability amplitude of the initial state |n 0 ⟩. In Figure 4, we plot the moving average of the survival probability within temporal windows of constant size together with the saturation value of the survival probability of GOE matrices with D = 4186. From Figure 4a, it can be seen that the long time limit of the survival probability at γ/J = 2.5 is in close proximity to the GOE. We can also see that the system significantly deviates from GOE when γ/J is away from 2.5. The evolution of the survival probability for different U is shown in Figure 4b. In the long time limit, the survival probability approaches that of the GOE when U ≈ 7.0. This seemingly contradicts the prediction based on β and entropy. However, we note that the system is already in the chaotic regime with the two parameters. It is not so surprising that the long-time survival probability is so close. On the other hand, the survival probabilities are generally close to each other in all the cases (note the logarithmic scale in Figure 4), where values of the survival probability are very small. These observation means that they are difficult to observe in cold atom experiments. In the following, we will show that variances of populations in different potential wells can be used to identify chaotic dynamics. We will show that the population variance exhibits different behaviors in the chaotic region by varying parameter U and γ. Changes in the variance are sizable, which might provide a plausible way to probe signatures of quantum chaos in the triple-well system. We numerically evaluated population ⟨n i ⟩ in the i-th site with initial state |N/3, N/3, N/3⟩. The time-averaged variance σ 2 of the populations in the three wells is defined as where ∆t = t f − t 0 . Figure 6 shows dependences of the variance on γ and with U/J = 7.0.
Here, σ 2 first increases, gains a peak value at around γ ≈ 2.5 and then decreases with an increasing γ. This dependence is similar to that of the Lyapunov exponent, i.e., Figure 1c, chaos indicator, i.e., Figure 3c, and entanglement entropy, i.e., Figure 4d. Hence, the similar profiles indicate that one could measure the variance to identify chaotic dynamics. Furthermore, we show the variance as a function of U/J in Figure 7. We find that the variance first increases with an increasing U/J and then saturates at around 90. Hence, the variance becomes insensitive if we further increase U. This could be attributed to the fact that strong interactions suppress hopping and hence reduce the population variance. Though the profile is different from other chaos indicators in this case, we notice that the saturation starts from U/J = 7. This corresponds to the parameter where other chaos indicators arrive at their maximal values.

Discussion and Conclusions
In this work, we studied quantum dynamics of Rydberg-dressed bosons within a triple-well setup. In the semiclassical regime, we identified parameter regions where Lyapunov exponents are large. Through diagonalizing the Hamiltonian, we showed that the level statistics fulfil the Wigner-Dyson distribution in the parameter region where the Lyapunov exponents are maximal, which indicates the emergence of quantum chaos in this system. By analyzing the quantum many-body dynamics, we showed that time-averaged entanglement entropy has a similar dependence on U and γ. We showed that the survival probability exhibits a similar dependence on U and γ to that of the level statistics and Lyapunov exponents. Our numerical calculations show that the time-averaged variance of the population could be used to provide chaotic dynamics. Acknowledgments: T.Y. would like to thank T. Hamlyn for useful discussions and critical reading of the draft. We thank Lea Santos for a critical comment on the saturation value in Figure 4.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Saturation Values of S EE
In this appendix, we derive the maximal entanglement entropy following Ref. [57]. Given a bosonic system bipartited into two subsystems A and B, the full Hilbert space H can be written in the form in which n A and n B are numbers of particles in subsystem A and B.
A quantum state |ψ⟩ can be expanded using the Fock basis, in which d , respectively, and form the computational basis of H For the sake of simplicity, we can project |ψ⟩ into a certain sector of n A forming |ψ n A ⟩ so that we do not need to consider the first summation in Equation (A2) and we have |ψ n A ⟩ in the form By applying the singular value decomposition (SVD) onto the d in which U is a d diagonal matrix with non-negative real numbers and V is a d semi-unitary matrix. After applying Schmidt decomposition, Equation (A4) can be rewritten as in which d (n A ) = min{d are the first d (n A ) columns of matrices U and V, respectively. The entanglement entropy S (n A ) EE with respect to n A -sector can therefore be easily calculated: The total entanglement entropy S EE can be computed by summing over all possible sectors: By applying Jensen's theorem onto Equation (A7), the upper bound of S EE can be found [57]: (A8) For a bosonic system with L sites in total and L A sites in subsystem A, Equation (A8) can be rewritten: In our triple-well setup, subsystem A only contains the first site. Therefore the upper bound of entanglement entropy only depends on the total number of particles: