Next Article in Journal
A PSO-SVM for Burst Header Packet Flooding Attacks Detection in Optical Burst Switching Networks
Next Article in Special Issue
Dynamical Collective Excitations and Entanglement of Two Strongly Correlated Rydberg Superatoms
Previous Article in Journal
On the Outage Capacity of Transdermal Optical Wireless Links with Stochastic Spatial Jitter and Skin-Induced Attenuation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Hyperchaos in a Bose-Hubbard Chain with Rydberg-Dressed Interactions

1
School of Physics and Astronomy, Centre for the Mathematics and Theoretical Physics of Quantum Non-Equilibrium Systems, University of Nottingham, Nottingham NG7 2RD, UK
2
Department of Physics, Indian Institute of Science Education and Research, Pune 411 008, India
*
Authors to whom correspondence should be addressed.
Photonics 2021, 8(12), 554; https://doi.org/10.3390/photonics8120554
Submission received: 9 November 2021 / Revised: 30 November 2021 / Accepted: 2 December 2021 / Published: 5 December 2021
(This article belongs to the Special Issue Optical Quantum Manipulation of Rydberg Atoms)

Abstract

:
We study the chaos and hyperchaos of Rydberg-dressed Bose–Einstein condensates (BECs) in a one-dimensional optical lattice. Due to the long-range, soft-core interaction between the dressed atoms, the dynamics of the BECs are described by the extended Bose-Hubbard model. In the mean-field regime, we analyze the dynamical stability of the BEC by focusing on the ground state and localized state configurations. Lyapunov exponents of the two configurations are calculated by varying the soft-core interaction strength, potential bias, and length of the lattice. Both configurations can have multiple positive Lyapunov exponents, exhibiting hyperchaotic dynamics. We show the dependence of the number of the positive Lyapunov exponents and the largest Lyapunov exponent on the length of the optical lattice. The largest Lyapunov exponent is directly proportional to areas of phase space encompassed by the associated Poincaré sections. We demonstrate that linear and hysteresis quenches of the lattice potential and the dressed interaction lead to distinct dynamics due to the chaos and hyperchaos. Our work is relevant to current research on chaos as well as collective and emergent nonlinear dynamics of BECs with long-range interactions.

1. Introduction

Over the past two decades, Bose–Einstein condensates (BECs) of ultracold atomic gases have become an ideal system for studying both quantum and nonlinear dynamics due to the high controllability of the two-body interactions [1], trapping potentials [2], and spatial dimensions [3,4], along with long coherence times. The emerging nonlinear phenomena depend strongly on the two-body interactions between atoms. In the presence of s-wave interactions, BECs can form dark and bright soliton [5,6,7,8,9,10,11,12] and exhibit Newton’s cradle behavior [13], which are paradigmatic examples in nonlinear physics. In trap array and optical lattice settings, the self-trapping of the BEC emerges due to strong repulsive interactions [14,15,16,17,18,19,20,21,22,23], where the BEC is localized within a single site. This is in contrast to the homogeneous, superfluid state, which form the ground state of an infinite lattice when the interaction is weak [24,25,26]. Both the homogeneous and self-trapped states correspond to solutions, i.e., fixed points, of the discrete Gross–Pitaevskii (GP) equation [27], which is a nonlinear Schrödinger equation that governs the mean-field dynamics. The stability of these fixed points depend on various parameters such as the s-wave interaction. It has been shown that the self-trapped state in a double-well potential can only be stable when the onsite interaction strength is much stronger than the tunneling strength [20]. Nonetheless, the homogeneous state can be disturbed by the s-wave interaction and external potentials, giving rise to chaotic dynamics [28,29]. Under strong periodic modulation of the hopping, extended chaotic regions are found in phase space [30].
On the other hand, long-range interactions play important roles in determining the dynamical stability of BECs. Solitons may occur in BECs in the presence of dipolar interactions [31,32,33,34,35]. The competition between s-wave and dipolar interactions [36,37,38] leads to bifurcations of the eigenspectra and chaotic dynamics when confined in a harmonic trap [39,40]. The self-trapping of dipolar BECs in double-well [41,42,43,44] and triple-well potentials [45,46] have been examined theoretically. Besides the dipolar interaction, one can laser couple ground state atoms to high-lying Rydberg states [47,48,49,50,51,52,53,54], which induces a long-range, soft-core interaction between two dressed atoms (with a distance r). The soft-core interaction is constant when r is within the soft-core radius R, typically of the order of several micrometers [48]. For r > R , the interaction decreases rapidly as r 6 , shown in Figure 1a. Various theoretical studies on the static and dynamical properties of Rydberg-dressed atoms confined in harmonic traps [55,56,57,58] and optical lattices [59,60,61,62,63,64,65] have been conducted in the past decade. Rydberg-dressed interactions have been experimentally demonstrated in optical tweezers [66], optical lattices [67,68,69], and harmonic traps [70], using various atomic species such as Rb [67,68], Cs [66,70], and Li [69] atoms. In Ref. [71], we showed that the self-trapping dynamics of Rydberg-dressed BECs can be controlled in a triple-well potential through mean-field and quantum mechanical analyses.
In this work, we investigate chaotic properties of Rydberg-dressed BECs in a one-dimensional (1D) optical lattice in which the dressed interaction leads to a multi-site density–density interaction. In the semiclassical regime, the nonlinear dynamics of the Bose-Hubbard model is captured by a discrete, coupled GP equation. Nonlinear eigenenergies, Bogoliubov spectra as well as Lyapunov exponents of the dressed BEC in the lattice are investigated. We then explore the dynamical stability of the ground state and localized state, where the dependence of the largest and the total number of positive Lyapunov exponents [72,73] on the dressed interaction and system size are explored. We show that at least two positive Lyapunov exponents can be found when the interaction is strong and when the number of sites is equal to or larger than four. This leads to the emergence of hyperchaos in the dynamics. We probe the chaotic dynamics by employing both a linear and a hysteresis quench of the potential bias and dressed interaction [74,75,76].
The paper is organized as follows. In Section 2, the Hamiltonian of the Bose-Hubbard Chain is introduced. The corresponding mean-field approximation and GP equations are given. Methods on calculating the eigenenergy, Bogoliubov spectra, and Lyapunov exponents are briefly introduced. Quench schemes of the potential bias and nonlinear interaction are explained. We explore the static (eigenenergies and Bogoliubov spectra) and dynamical properties (Lyapunov exponents) of the ground state and localized state configurations in Section 3 and Section 4, respectively. Dynamics driven by both the linear and hysteresis quenching parameters are explored with different initial states. In Section 5, we examine the scaling of the Lyapunov exponents with the system size for the two different configurations. We demonstrate, through numerical calculations, that areas of the Poincaré sections are almost linearly proportional to the largest Lyapunov exponent. We conclude our work in Section 6.

2. Model and Method

2.1. Extended Bose-Hubbard Model in the Semiclassical Limit

Our setting consists of a gas of N bosonic atoms at zero temperature confined in a one-dimensional lattice, with the lattice constant d, as depicted in Figure 1a. The Rydberg dressing induces long-range interactions between atoms at different sites. Taking into account the hopping between nearest-neighbor sites, we obtain an extended Bose-Hubbard Hamiltonian of L sites [52] ( = 1 ):
H ^ = J i , j L a ^ i a ^ j + j L Γ j n ^ j + 1 2 H ^ int ,
where a ^ j ( a ^ j ) is the bosonic annihilation (creation) operator at site j. The tunneling strength J acts only on nearest-neighbor sites, denoted by · in the summation. Here, n ^ j = a ^ j a ^ j is the number operator, while Γ j is the local tilting potential. The titling is given by Γ j = γ j 1 L / 2 , where · and γ are the floor function and level bias between neighboring sites, respectively. The onsite and long-range interactions are described by H ^ int = g j L n ^ j n ^ j 1 + i , j L Λ i , j n ^ i n ^ j . The onsite interaction g = 4 π a s / m [2] depends on the s-wave scattering length a s and mass m, where the former can be adjusted by Feshbach resonances [2]. The soft-core shaped long-range interaction is given by Λ i , j = C 6 / [ | i j | 6 d 6 + R 6 ] , with C 6 being the dispersion coefficient (Figure 1a). Both the soft-core radius R and C 6 can be tuned by laser parameters [48]. In this work, we will restrict our investigation to the onsite, nearest-neighbor ( Λ j , j ± 1 ) and next-nearest-neighbor ( Λ j , j ± 2 ) interactions only, where R d . This approximation is valid as the soft-core interaction decays rapidly when the separation between sites is larger than the soft-core radius.
In the semiclassical limit N 1 , we employ the mean-field approximation where the bosonic operator is described by a classical field ψ j , i.e., a ^ j ψ j N and a ^ j ψ j * N , with the normalization condition j | ψ j | 2 = 1 . This yields the semiclassical Hamiltonian H H ^ / N :
H = j L Γ j | ψ j | 2 J j L ψ j + 1 * ψ j + ψ j 1 * ψ j + N 2 i , j L | ψ j | 2 g | ψ j | 2 1 + Λ i , j | ψ i | 2 .
The dynamics of the classical field ψ j is obtained via the canonical equation i d ψ j / d t = H / ψ j * , yielding the coupled GP equations:
i ψ ˙ j = J ψ j + 1 + ψ j 1 + [ Γ j + W | ψ j | 2 + U | ψ j + 1 | 2 + | ψ j 1 | 2 ) + V ( | ψ j + 2 | 2 + | ψ j 2 | 2 ] ψ j ,
where we have defined W = N ( Λ j , j + g ) , U = N Λ j , j ± 1 , and V = N Λ j , j ± 2 to be the onsite, nearest-neighbor, and next-nearest-neighbor interaction strengths. The onsite interaction W takes into account contributions from both the s-wave and soft-core interactions. We will assume a vanishing onsite interaction, i.e., W = 0 , which allows us to focus on effects induced by the long-range interaction part. To be concrete, we will fix the nearest-neighbor and next-nearest-neighbor interactions at U = 2 V in the following discussion. Time and energy will be scaled with respect to 1 / J and J in what follows.
It is convenient to examine the real ( R j = Re [ ψ j ] ) and imaginary components ( I j = Im [ ψ j ] ) of ψ j :
R ˙ j = J I j + 1 + I j 1 + [ Γ j + W | ψ j | 2 + U | ψ j + 1 | 2 + | ψ j 1 | 2 ) + V ( | ψ j + 2 | 2 + | ψ j 2 | 2 ] I j ,
I ˙ j = J R j + 1 + R j 1 [ Γ j + W | ψ j | 2 + U | ψ j + 1 | 2 + | ψ j 1 | 2 ) + V ( | ψ j + 2 | 2 + | ψ j 2 | 2 ] R j ,
with | ψ j | 2 = R j 2 + I j 2 . Both R j and I j are real-valued functions of time. We will calculate Lyapunov exponents and the Poincaré sections based on these real functions. Note that R j and I j represent mean values of the quadrature of the operator a ^ j . The quadrature fulfills the commutation relation similar to the position and momentum operator [77]. Hence, the mean values of the quadrature allow us to obtain useful information on the dynamics of the system in phase space. For small systems, L = 2 or 3, one can also describe the classical field with the canonical phase and particle number decomposition [15,78].

2.2. Nonlinear Eigenenergies and Bogoliubov Spectra

Though the Hamiltonian (2) is Hermitian, the density-dependent nonlinearity prevents us from calculating the eigenenergy through conventional diagonalization. To overcome this, a shooting method will be employed to numerically evaluate the eigenstate Ψ ¯ j = [ ψ ¯ 1 , ψ ¯ 2 , , ψ ¯ L ] and the corresponding eigenenergy ε j self-consistently [71]. A trial solution is seeded into the semiclassical Hamiltonian. It is then diagonalized, leading to a new eigenstate and eigenenergy. This process is iterated until the resulting eigenstate and eigenenergy are obtained self-consistently.
For interacting systems, one can analyze the Bogoliubov spectra ϵ B to understand the stability of the eigenstate. This is achieved by linearizing around a given state Ψ ¯ (e.g., a fixed point of the semiclassical system), where each component is given by ψ j = ψ ¯ j + u j e i ϵ B t v j * e i ϵ B t , with u j and v j being the probability amplitudes of the Bogoliubov quasiparticles [2]. The dynamics of u j and v j are described by the Bogoliubov equations [79,80]:
L N N L u v = ϵ B u v
where L = H ˜ 0 + 2 U P μ and N = U P . H ˜ 0 and P are L × L block matrices. From Equation (3), we obtain the matrix elements ψ j | H ˜ 0 | ψ j = Γ j , ψ j | H ˜ 0 | ψ j ± 1 = J , ψ j | P | ψ j = | ψ ¯ j + 1 | 2 + | ψ ¯ j 1 | 2 + ( | ψ ¯ j + 2 | 2 + | ψ ¯ j 2 | 2 ) / 2 , ψ j | P | ψ j ± 1 = 2 ψ ¯ j ± 1 ψ ¯ j and ψ j | P | ψ j ± 2 = ψ ¯ j ± 2 ψ ¯ j , while other matrix elements are zero. If the Bogoliubov spectra are complex numbers, the state is then dynamically unstable as Bogoliubov quasiparticles grow (decay) exponentially with time, the rate of which is determined by the imaginary part of the spectra.

2.3. Poincaré Sections and Lyapunov Exponents

The emergence of chaos in the dynamics can be characterized by the Poincaré sections and Lyapunov exponents. For L sites, the possible trajectories are the complete set of { R 1 , , R L , I 1 , , I L } . Due to the normalization condition, we need to solve a 2 L 1 dimensional system to obtain the dynamics. It is difficult to comprehend the stability of the trajectories in such a high-dimensional phase space. Instead, we project the dynamics to a two-dimensional (2D) Poincaré section to identify the dynamical properties. To calculate the 2D Poincaré section, we record trajectories of selected variables ( R j , I j ) as they cut through the U k plane ( j k ), provided that U ˙ k > 0 . These intersecting points form the 2D Poincaré section. To be specific, we will evaluate the Poincaré section of the variable ( R 2 , I 2 ) on the U 1 plane.
The strength of the chaos can be measured by the Lyapunov exponents associated with the equations of motion [72,73]. The Lyapunov exponents give the rate of separation between trajectories for a given initial state. As the Lyapunov exponents depend on the initial state, we initially consider both the ground state and a localized state. In a localized state, nearly all the condensate sits in a single site, which can be stable (i.e., the self-trapping state) when the nonlinear interaction is strong. In this work, the Lyapunov exponents λ j ( j = 1 , , 2 L ) are calculated via DynamicalSystems.jl, a fast and reliable Julia library, to determine the dynamics of nonlinear systems [81]. We confirmed with the method in Ref. [72] that it gives consistent data. When there exists at least one positive Lyapunov exponent, the trajectories will separate exponentially, leading to chaotic dynamics. The dynamics is hyperchaotic when there are more than two positive Lyapunov exponents [73].

2.4. Quenching Schemes

In Section 3 and Section 4 we will explore the dynamics of the system with time-dependent parameters via the following quenching schemes:
Scheme I : First, we consider a linear quench of the potential bias [71]. The bias between two neighboring sites is given by the following function:
γ L = γ i + α t ,
where γ i and α are the initial value and quench rate, respectively. With γ i < 0 , the quench takes place from t = 0 to t = 2 γ f / α with γ f = γ i , depicted by the solid curve in Figure 1b.
Scheme II : Alternatively, we consider a hysteresis quench [74,75,76] where the system begins at γ i and then evolves to γ f . At time τ = γ f / α , the potential bias is quenched back towards γ i . The function describing this scheme is as follows:
γ H = γ f + ( γ i γ f ) | τ t | τ .
The corresponding scheme is shown by the solid and dashed curves in Figure 1b.
Scheme III : In addition to quenching the level bias, we also change the two-body interaction strength through a linear ramp:
U L = U i + β t
where U i is the initial interaction strength and β is the quench rate. This is shown by the solid curve in Figure 1c. Note that the next-nearest-neighbor interaction V also depends on time due to the relation U = 2 V .
Scheme IV : The hysteresis counterpart of the interaction quench is given by the following equation:
U H = U f + ( U i U f ) | τ t | τ ,
where U f is the final interaction strength, with τ = U f / β .

3. Stability of the Ground State

3.1. Eigenenergies, Bogoliubov Spectra, and Lyapunov Exponents

Without the nonlinearity, the number of eigenenergies N ϵ is identical to L, the dimension of the semiclassical system. The number of eigenenergies can be larger than L when the interaction is strong. As an example, eigenenergies for L = 3 as a function of the bias γ are shown in Figure 2a. We find that N ϵ > L when | γ | U , where the nonlinearity dominates. Loops and crossings appear in the eigenenergies, except the highest energy level.
For a given state of the nonlinear system, one obtains 2 L Bogoliubov spectra, whose values depend on the specific eigenstate and nonlinear interaction strength. The Bogoliubov modes are stable for all γ when the system is in the ground state, i.e., the Bogoliubov spectra ϵ B are real, as shown in Figure 2b. This is in contrast to excited eigenstates, whose Bogoliubov spectra have imaginary components. As an example, the Bogoliubov spectra of the first excited state are shown in Figure 2c. The corresponding Bogoliubov mode will decay (grow) exponentially, when the imaginary part is negative (positive).
The chaotic dynamics of the system is characterized by positive Lyapunov exponents. In Figure 2d, Lyapunov exponents are shown for the ground state of the system. When increasing γ , negative and positive Lyapunov exponents are found in regions where the eigenenergies show loops. The negative and positive Lyapunov exponents appear in pairs with the same absolute values as our system is conservative. In this example, one positive Lyapunov exponent can be found when | γ | < 1 , indicating the presence of chaos. This means that small fluctuations on the ground state could gain exponential growth, and hence drives the system away from the ground state.
To further understand the roles played by the nonlinearity, we calculate eigenenergies as a function of the interaction strength U shown in Figure 3a,b for L = 3 and L = 5 , respectively. It can be seen that new branches are generated when the nonlinear interaction U is large enough. Lyapunov exponents of the ground state of the nonlinear system are shown in Figure 3c,d. Positive Lyapunov exponents are found in the strongly interacting region, whose values increase with increasing U. Larger Lyapunov exponents mean that the exponential growth of the instability can be even faster. More importantly, the number of Lyapunov exponents now depends on L. For L = 3 , one obtains a single positive Lyapunov exponent when U 4 . When L = 5 , there are three positive Lyapunov exponents. This indicates that the system enters the so-called hyperchaos regime [82,83,84], where more than one positive Lyapunov exponents can be found in the dynamics. In the two examples, we obtain maximally L 2 positive Lyapunov exponents as the energy and particle number are conserved in the Bose-Hubbard Chain.

3.2. Quench Dynamics

In the linear regime, dynamics of the system will follow the eigenstate adiabatically when slowly quenching the tilt potential. However, the dynamics may deviate from the adiabatic eigenstate in the nonlinear regime, especially when positive Lyapunov exponents are found. This will be illustrated through the quenching of the tilt potential and interaction strength given by Equations (7)–(10). To trigger the instability in the dynamics, we consider a thermal mixed state Ψ ¯ j = [ ψ ¯ 1 e i θ 1 , ψ ¯ 2 e i θ 2 , , ψ ¯ L e i θ L ] around a given state Ψ ¯ j (the ground state), where θ j is a random phase distributed uniformly between 0 and 2 π [76]. In numerical simulations, we typically consider an ensemble of M = 100 realizations with a given set of parameters.
We first examine a linear quench of the bias γ when the system is prepared in the ground state at γ i = 10 and L = 3 . The majority of the condensate is located on the first (leftmost) site [ n 1 ( 0 ) = | ψ 1 | 2 1 ] initially (Figure 4a). In the adiabatic limit and without the nonlinear interaction, the condensate will move to the third well, n 3 ( τ ) 1 , after the quench [71]. The population at this adiabatic limit is shown with a black dot in each panel. The quench dynamics, however, depend on the finite quench rate and the interaction strength. When the interaction is weak, the condensate can be in any of the three sites since the tunneling strength between neighboring sites plays the dominant role. The distribution of the final population is affected by the noise on the initial state and also depends on the final time in the simulations. By increasing U, the population is distributed into a larger region of the phase space, i.e., it occupies larger areas in the n 3 n 1 plane. By fixing the interaction U, our numerical simulation shows that the smaller α is, the closer the population distribution is to the adiabatic limit.
For the hysteresis quench given by Equation (8), we see that even for U = 5 (meaning the eigenstate exhibits complicated level crossings) the density mostly returns to their initial state, at least when α 1 (Figure 4b). Here, the hysteresis quench has allowed for a large level of reversibility in the dynamics [76] as the chaotic regions have not been triggered. By increasing the quench rate α , the population distributions cluster around much smaller regions in the phase space than the one shown in panel (a).
In Figure 4c we quench the interaction according to Equation (9). The initial states depend on the value of γ . For example, the ground state is Ψ ¯ = [ 0.5 , 1 / 2 , 0.5 ] for γ = 0 . The final states are highly dependent on the initial conditions due to the chaos in the dynamics (see the crossing energy levels in Figure 3a and Lyapunov exponent in Figure 3c). We have verified that by increasing γ , the associated randomness with the final states decreases since the number of crossings in the eigenenergy also decreases.
In case of the hysteresis quench of U, we find that the results (Figure 4d) are similar to the linear quench. When looking at γ = 0 , the final states do not return to the initial value. As shown in Figure 3c, the Lyapunov exponent of the ground state becomes positive when U 4 , which causes the final state to be more random, i.e., to have a broader distribution of the densities. As the tilt γ increases, we have verified that chaos is gradually suppressed as the population that localizes in the trap corresponds to the lowest energy state. In order to trigger chaotic dynamics in the tilted case, stronger interactions are generally needed.

4. Stability of the Localized State

4.1. Bogoliubov Spectra and Lyapunov Exponents

In this section, we will explore stability of a situation where the condensate is trapped in a single site. When localized at one end of the lattice, it corresponds to the ground state if the lattice potential is strongly tilted | γ | 1 . We will examine dynamics of localized states even in the balanced case ( γ = 0 ), partially motivated by the fact that the self-trapped state can be stabilized by strong nonlinear interactions. We will show that dynamical instabilities of localized states will depend strongly on the long-range interaction. To be concrete, we will consider a scenario where the condensate is confined in the second trap from the left of the lattice, i.e., Ψ ¯ = [ 0 , 1 , , 0 ] . In the numerical simulations of the dynamics, uniform density fluctuations are applied to the lattice to trigger the hopping dynamics. This modifies the initial state into Ψ ¯ = [ ε / L e ϕ 1 , 1 ε e ϕ 2 , , ε / L e ϕ L ] with ϵ 1 and ϕ j being a random phase. Furthermore, this choice insures that the energy of different initial states are almost identical.
In Figure 5a,b, dynamical unstable regions in the Bogoliubov spectra for L = 3 and L = 5 are shown (highlighted with a dark red color). In the unstable region, ϵ B develops imaginary components, which depend on U, γ and L. In the case of L = 3 , the condensate is localized in the middle site initially, meaning the Bogoliubov spectra are symmetric with respect to γ . Figure 5a shows that the system is dynamically unstable when U is small, in particular when the lattice is balanced ( | γ | is small). This is not surprising since the localized state is not the ground state, nor does the system support the self-trapped state. By increasing the interaction strength, we note that the localized state returns to a stable configuration when | γ | is small. This means that the localized state becomes a stable, self-trapped state [71]. When L = 5 , the dynamical stability now depends heavily on tilt γ . When γ > 0 , there is a much broader range of unstable regions. This feature is largely due to the fact that the nonsymmetric initial state has higher energies. Therefore, we expect to see qualitatively different dynamics from the various quenching schemes.
The Lyapunov exponents exhibit sensitive dependence on the system size. As shown in Figure 5c, the Lyapunov exponents for L = 3 show an unusually symmetric shape when U 4 . The exponents are a smooth function of U and reach maximal value around U = 5 . By further increasing U, the positive Lyapunov exponents decrease. This indicates that the localized configuration could exhibit chaotic dynamics for a large U. For L = 5 we notice that positive Lyapunov exponents can be found when U is relatively small. A key difference is that there are multiple positive Lyapunov exponents (Figure 5d) where the nonlinear dynamics enters the hyperchaotic regime.
To understand the maximal Lyapunov exponents, we slightly alter the initial state so that we have Ψ ¯ = ε , 1 2 ε 2 , ε , , 0 , where ε is a small perturbation to the wave function of the traps on either side of the localized site, with 0 < ε < 0.01 . In Figure 5e,f (corresponding to L = 3 and 5), the largest Lyapunov exponent λ m (red) and Lyapunov exponents obtained with modified initial states (black) are shown (only the positive branch). It shows that a minor change to the initial state will change Lyapunov exponents significantly. However, λ m gives an approximate upper bound for all the Lyapunov exponents.

4.2. Quench Dynamics

For U = 5 and L = 3 , a linear quench (Figure 6a) from γ i = 10 to γ f = 10 shows strong, self-trapping behavior in the rightmost potential. Ideally, we would expect that by performing a hysteresis quench back towards γ i , the population would localize in the leftmost site again. However, from Figure 6b, we see that the final state is rather chaotic. Due to the dynamical instability and chaos near | γ | < 1 , the final state deviates from the initial state. In panels (c) and (d), we quench according to Equations (9) and (10), respectively. The dynamics shows that in both cases, the localized initial state loses population to the outer potential wells in an approximately equal manner for both the linear and hysteresis quenches. The strong, nonlinear interactions in the initial localized trap repel the condensate symmetrically between the two neighboring traps. Additionally, we notice that in panel (d), the population could be n 1 = n 3 0 , meaning that the final state is exactly equal to the initial state. We have achieved full reversibility with the hysteresis dynamics in these simulations. As shown in panel (c), this is not the case where the populations are always n 1 n 2 > 0 , implying that the strong, two-body interactions prevent the complete localization of the condensate in a single site.
We now move on to examine the dynamics for the five-site system. Without two-body interactions, linearly quenching from γ i to γ f will force the atoms towards the rightmost trap. However, from Figure 6e, we see that the occupation is never very much greater than n 5 0.5 , even for the slowest quenching rates considered in the simulation. When the quench rate is fast ( α 1 ), we find less occupation in both the first and last sites, implying that the occupation has been spread amongst the remaining sites. In panel (f), the hysteresis counterpart is shown. Now the population should tend towards n 1 1 . However, this is not what is found in the numerical simulations. The populations distribute randomly in all sites. In panels (g) and (h), the dynamics is qualitatively different from the L = 3 scenario. The symmetry between the densities of the two outermost sites is completely lost and is replaced with a chaotic distribution, largely due to the presence of hyperchaos (see Figure 5d).

5. Discussion

In the following, we will discuss how the maximal and total number of Lyapunov exponents depend on the system size and initial state, focusing on parameter regimes where the nonlinear interaction cannot be neglected, i.e., chaos and hyperchaos are expected in the dynamics. In general, Lyapunov exponents depend on the input state of the calculation [72]. Two different initial states, i.e., the ground state and the localized state, will be examined in detail.
In Figure 7a the largest Lyapunov exponent λ m for the ground state configuration is shown. When 2 L 4 , the values of λ m are generally small. This is due to the fact that chaos has not been triggered (see Figure 5b,c). When L > 4 , the situation changes as chaos is already found with the given U. We find that λ m decreases gradually when U = 3 and U = 5 for larger L. On the other hand, the total number of positive Lyapunov exponents η is seen to increase almost linearly with L when U = 3 and U = 5 , depicted in Figure 7b. More importantly, η > 2 when L > 4 for both U = 3 and U = 5 , i.e., the dynamics is hyperchaotic. On the other hand, η decreases and deviates from the linear dependence on L when L is large, e.g., at L = 10 when U = 3 and L = 14 when U = 5 . In general, the linear relation holds up to a larger L for larger U. It has recently been shown that the largest Lyapunov exponents in the BH model can be obtained from the echo dynamics of the condensate [85]. A similar technique could be applied to extract the largest Lyapunov exponents studied here.
Figure 7c,d show both λ m and η for the the localized state. In this case, λ m is largest when L = 5 and decreases with increasing L for U = 3 and U = 5 . Compared to the ground state, a visible difference is that λ m 0 when U = 1 for the localized state. Their values, however, are smaller than the one for U = 3 and U = 5 . This implies that it will be difficult to observe chaotic dynamics with this level of nonlinear interactions. On the other hand, η increases with increasing U. When L > 10 , η still increases with L, slightly deviating from the linear scaling with L. A similar dependence is also found for stronger nonlinear interactions, as demonstrated with U = 5 in panel (d). For such states, η > 1 can be seen even with a relatively weak interaction (e.g., U = 1 ), leading to more pronounced hyperchaotic dynamics.
The total number of nonlinear differential equations is 2 L (the real and imaginary parts of ψ j ). For conservative systems, the number of positive and negative Lyapunov exponents are the same, and the sum of the Lyapunov exponents is zero. These features can be seen, e.g., in Figure 2d. Our numerical simulations show that the maximal number of positive Lyapunov exponents is L 2 (see Figure 7b when U = 5 and L 11 and Figure 7d when U = 5 and L 14 ). As the extended Bose-Hubbard model is a Hamiltonian system, not only does the sum of the Lyapunov exponents vanish, but conservative quantities such as the energy and particle number are also found in the dynamics. This indicates that the maximal number of the Lyapunov exponents is L 2 and not L. For a sufficiently large L, the total number of positive Lyapunov exponents is smaller than L 2 as the nonlinear interaction becomes smaller. For the ground state, one can estimate the interaction energy for a given site to be approximately 2 ( U + V ) / L 2 , i.e., the mean local interaction energy decreases with increasing L.
The chaotic dynamics depends strongly on the largest Lyapunov exponents λ m , which is considered to be an indication of chaos in the dynamics. To illustrate this, the Poincaré section on the U 1 plane for different system sizes is shown in Figure 8, showing that profiles of the Poincaré section depend on the system size and the initial state. The area is largest when L = 5 (Figure 8a) and decreases with increasing L in the case of the ground state (Figure 8b,c). For different L, the profile of the Poincaré section is largely symmetric with respect to R 2 = 0 and I 2 = 0 . In the case of the localized state, similar dependence on L is found, as depicted in Figure 8d–f. We note two differences when compared to the ground state. First, the profile of the Poincaré section displays symmetry with respect to I 2 = 0 but not R 2 = 0 . Second, the areas of the Poincaré section in the localized state are slightly larger as the corresponding λ m is larger (see Figure 7a,c).
The area is largely determined by the largest Lyapunov exponent. To verify this, we find the area of the Poincaré section approximately by numerically fitting the Poincaré section, as shown in Figure 9. For the ground state, the dependence of the fitted area and λ m on L agrees well when U = 5 . For U = 3 , a good agreement is also found when L 8 . When L = 5 and L = 7 , the fitted areas differ largely from the corresponding λ m . This discrepancy might be caused by the fact that the relatively weak nonlinear interaction leads to uncertainties in calculating the Lyapunov exponent. For the localized state, the agreement is improved in general for both U = 3 and U = 5 . This suggests that the discrepancy in the ground state could be a boundary effect when L is small, as the localized state suffers less from the boundary effect.

6. Conclusions and Outlook

We investigated the chaotic and hyperchaotic dynamics of a one-dimensional Bose-Hubbard Chain of Rydberg-dressed BECs in the semiclassical regime. We showed that both the ground state and localized state could have positive Lyapunov exponents even though the corresponding Bogoliubov spectra were real-valued. As a result, small perturbations to these states led to large fluctuations, which were captured by the quench dynamics. The vastly different population distribution after the quench provided evidence of chaos where the former could be measured experimentally. We found that hyperchaos emerges in both the ground and localized states when the nonlinear interaction is strong and L is large. The total number of positive Lyapunov exponents, η , is bound by L 2 ( L 3 ). We showed that η grows with the system size L when U is large. So far, our investigations have been focused on the semiclassical regime. There has been exploration into the relationships between chaos and quantum entanglement [86]. Moreover, quantum chaos can be seen by analyzing the statistics of the eigenspectra on the Bose-Hubbard model with onsite [87] and long-range interactions [88,89]. It is therefore worthwhile to explore chaos and hyperchaos as well as the thermalization of the Rydberg atoms in the quantum regime.

Author Contributions

Conceptualization, G.M. and W.L.; methodology, G.M.; formal analysis, G.M. and W.L.; writing—original draft preparation, G.M.; writing—review and editing, W.L. and R.N. All authors have read and agreed to the published version of the manuscript.

Funding

The research leading to these results received funding from EPSRC Grant No. EP/R04340X/1 via the QuantERA project “ERyQSenS”, UKIERI-UGC Thematic Partnership No. IND/CONT/G/16-17/73, and the Royal Society through International Exchanges Cost Share Award No. IEC∖NSFC∖181078. R.N. acknowledges DST-SERB for a Swarnajayanti fellowship (File No. SB/SJF/2020-21/19).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Robert, P.S.; Zoran, H. Effects of Interactions on Bose-Einstein Condensation of an Atomic Gas. In Physics of Quantum Fluids: New Trends and Hot Topics in Atomic and Polariton Condensates; Alberto, B., Michele, M., Eds.; Springer: Berlin/Heidelberg, Germany, 2013; pp. 341–359. [Google Scholar]
  2. Pethick, C.J.; Smith, H. Microscopic theory of the Bose gas. In Bose–Einstein Condensation in Dilute Gases; Cambridge University Press: Cambridge, UK, 2008. [Google Scholar]
  3. Leonardo, F.; Chiara, F.; Jessica, E.L.; Massimo, I. Bose-Einstein condensate in an optical lattice with tunable spacing: Transport and static properties. Opt. Express 2005, 13, 4303–4313. [Google Scholar]
  4. Smerzi, A.; Fantoni, S.; Giovanazzi, S.; Shenoy, S.R. Quantum coherent atomic tunneling between two trapped bose-einstein condensates. Phys. Rev. Lett. 1997, 79, 4950–4953. [Google Scholar] [CrossRef] [Green Version]
  5. Manjun, M.R.; Navarro, R.; Carretero-González, R. Solitons riding on solitons and the quantum newton’s cradle. Phys. Rev. E 2016, 93, 022202. [Google Scholar]
  6. Anderson, B.P.; Haljan, P.C.; Regal, C.A.; Feder, D.L.; Collins, L.A.; Clark, C.W.; Cornell, E.A. Watching dark solitons decay into vortex rings in a bose-einstein condensate. Phys. Rev. Lett. 2001, 86, 2926–2929. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Zachary, D.; Michael, B.; Christopher, S.; Lene Vestergaard, H. Observation of quantum shock waves created with ultra- compressed slow light pulses in a bose-einstein condensate. Science 2001, 293, 663–668. [Google Scholar]
  8. Denschlag, J.; Simsarian, J.E.; Feder, D.L.; Clark, C.W.; Collins, L.A.; Cubizolles, J.; Deng, L.; Hagley, E.W.; Helmerson, K.; Phillips, W.D.; et al. Generating solitons by phase engineering of a bose-einstein condensate. Science 2000, 287, 97–101. [Google Scholar] [CrossRef] [PubMed]
  9. Burger, S.; Bongs, K.; Dettmer, S.; Ertmer, W.; Sengstock, K.; Sanpera, A.; Shlyapnikov, G.V.; Lewenstein, M. Dark solitons in bose-einstein condensates. Phys. Rev. Lett. 1999, 83, 5198–5201. [Google Scholar] [CrossRef] [Green Version]
  10. Cornish, S.L.; Thompson, S.T.; Wieman, C.E. Formation of bright matter-wave solitons during the collapse of attractive bose-einstein condensates. Phys. Rev. Lett. 2006, 96, 170401. [Google Scholar] [CrossRef] [Green Version]
  11. Khaykovich, L.; Schreck, F.; Ferrari, G.; Bourdel, T.; Cubizolles, J.; Carr, L.D.; Castin, Y.; Salomon, C. Formation of a matter-wave bright soliton. Science 2002, 296, 1290–1293. [Google Scholar] [CrossRef] [Green Version]
  12. Strecker, K.E.; Partridge, G.B.; Truscott, A.G.; Hulet, R.G. Formation and propagation of matter-wave soliton trains. Nature 2002, 417, 150–153. [Google Scholar] [CrossRef]
  13. Kinoshita, T.; Wenger, T.; Weiss, D.S. A quantum newton’s cradle. Nature 2006, 440, 900–903. [Google Scholar] [CrossRef]
  14. Xia, B.; Hai, W.; Chong, G. Stability and chaotic behavior of a two-component Bose-Einstein condensate. Phys. Lett. Sect. A Gen. At. Solid State Phys. 2006, 351, 136–142. [Google Scholar] [CrossRef]
  15. Liu, B.; Fu, L.B.; Yang, S.P.; Liu, J. Josephson oscillation and transition to self-trapping for Bose-Einstein condensates in a triple-well trap. Phys. Rev. A 2007, 75, 033601. [Google Scholar] [CrossRef] [Green Version]
  16. Graefe, E.M.; Korsch, H.J.; Witthaut, D. Mean-field dynamics of a Bose-Einstein condensate in a time-dependent triple-well trap: Nonlinear eigenstates, Landau-Zener models, and stimulated Raman adiabatic passage. Phys. Rev. A 2006, 73, 013617. [Google Scholar] [CrossRef] [Green Version]
  17. Viscondi, T.F.; Furuya, K. Dynamics of a Bose-Einstein condensate in a symmetric triple-well trap. J. Phys. A Math. Theor. 2011, 44, 175301. [Google Scholar] [CrossRef] [Green Version]
  18. Li, L.; Wang, B.; Lü, X.Y.; Wu, Y. Chaos-related Localization in Modulated Lattice Array. Annalen der Physik 2018, 530, 1700218. [Google Scholar] [CrossRef]
  19. Chong, G.; Hai, W.; Xie, Q. Controlling chaos in a weakly coupled array of Bose-Einstein condensates. Phys. Rev. E 2005, 71, 016202. [Google Scholar] [CrossRef]
  20. Liu, J.; Fu, L.; Ou, B.Y.; Chen, S.G.; Choi, D.I.; Wu, B.; Niu, Q. Theory of nonlinear Landau-Zener tunneling. Phys. Rev. A 2002, 66, 023404. [Google Scholar] [CrossRef] [Green Version]
  21. Liu, J.; Wu, B.; Niu, Q. Nonlinear Evolution of Quantum States in the Adiabatic Regime. Phys. Rev. Lett. 2003, 90, 170404. [Google Scholar] [CrossRef] [Green Version]
  22. Albiez, M.; Gati, R.; Fölling, J.; Hunsmann, S.; Cristiani, M.; Oberthaler, M.K. Direct observation of tunneling and nonlinear self-trapping in a single bosonic josephson junction. Phys. Rev. Lett. 2005, 95, 010402. [Google Scholar] [CrossRef] [Green Version]
  23. Zibold, T.; Nicklas, E.; Gross, C.; Oberthaler, M.K. Classical bifurcation at the transition from rabi to Josephson dynamics. Phys. Rev. Lett. 2010, 105, 204101. [Google Scholar] [CrossRef]
  24. Gotlibovych, I.; Schmidutz, T.F.; Gaunt, A.L.; Navon, N.; Smith, R.P.; Hadzibabic, Z. Observing properties of an interacting homogeneous bose-einstein condensate: Heisenberg-limited momentum spread, interaction energy, and free-expansion dynamics. Phys. Rev. A 2014, 89, 061604. [Google Scholar] [CrossRef] [Green Version]
  25. Gaunt, A.L.; Schmidutz, T.F.; Gotlibovych, I.; Smith, R.P.; Hadzibabic, Z. Bose-einstein condensation of atoms in a uniform potential. Phys. Rev. Lett. 2013, 110, 200406. [Google Scholar] [CrossRef] [PubMed]
  26. Schmidutz, T.F.; Gotlibovych, I.; Gaunt, A.L.; Smith, R.P.; Navon, N.; Hadzibabic, Z. Quantum joule-thomson effect in a saturated homogeneous bose gas. Phys. Rev. Lett. 2014, 112, 040403. [Google Scholar] [CrossRef] [Green Version]
  27. Buonsante, P.; Penna, V. Some remarks on the coherent-state variational approach to nonlinear boson models. J. Phys. A Math. Theor. 2008, 41, 175301. [Google Scholar] [CrossRef]
  28. Hai, W.; Rong, S.; Zhu, Q. Discrete chaotic states of a Bose-Einstein condensate. Phys. Rev. E 2008, 78, 066214. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Sinha, S.; Sinha, S. Chaos and Quantum Scars in Bose-Josephson Junction Coupled to a Bosonic Mode. Phys. Rev. Lett. 2020, 125, 134101. [Google Scholar] [CrossRef] [PubMed]
  30. Boukobza, E.; Moore, M.G.; Cohen, D.; Vardi, A. Nonlinear phase dynamics in a driven bosonic josephson junction. Phys. Rev. Lett. 2010, 104, 240402. [Google Scholar] [CrossRef] [Green Version]
  31. Pedri, P.; Santos, L. Two-Dimensional Bright Solitons in Dipolar Bose-Einstein Condensates. Phys. Rev. Lett. 2005, 95, 200404. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Tikhonenkov, I.; Malomed, B.A.; Vardi, A. Anisotropic Solitons in Dipolar Bose-Einstein Condensates. Phys. Rev. Lett. 2008, 100, 090406. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Nath, R.; Pedri, P.; Santos, L. Phonon Instability with Respect to Soliton Formation in Two-Dimensional Dipolar Bose-Einstein Condensates. Phys. Rev. Lett. 2009, 102, 050401. [Google Scholar] [CrossRef] [Green Version]
  34. Cuevas, J.; Malomed, B.A.; Kevrekidis, P.G.; Frantzeskakis, D.J. Solitons in quasi-one-dimensional Bose-Einstein condensates with competing dipolar and local interactions. Phys. Rev. A 2009, 79, 053608. [Google Scholar] [CrossRef] [Green Version]
  35. Young-S, L.E.; Muruganandam, P.; Adhikari, S.K. Dynamics of quasi-one-dimensional bright and vortex solitons of a dipolar Bose–Einstein condensate with repulsive atomic interaction. J. Phys. B At. Mol. Opt. Phys. 2011, 44, 101001. [Google Scholar] [CrossRef] [Green Version]
  36. Lahaye, T.; Pfau, T.; Santos, L. Mesoscopic ensembles of polar bosons in triple-well potentials. Phys. Rev. Lett. 2010, 104, 170404. [Google Scholar] [CrossRef] [Green Version]
  37. Xiong, B.; Fischer, U.R. Interaction-induced coherence among polar bosons stored in triple-well potentials. Phys. Rev. A 2013, 88, 063608. [Google Scholar] [CrossRef] [Green Version]
  38. Gallemí, A.; Guilleumas, M.; Mayol, R.; Sanpera, A. Role of anisotropy in dipolar bosons in triple-well potentials. Phys. Rev. A 2013, 88, 063645. [Google Scholar] [CrossRef] [Green Version]
  39. Köberle, P.; Cartarius, H.; Fabčič, T.; Main, J.; Wunner, G. Bifurcations, order and chaos in the bose–einstein condensation of dipolar gases. New J. Phys. 2009, 11, 023017. [Google Scholar] [CrossRef]
  40. Andreev, P.A. Quantum hydrodynamic theory of quantum fluctuations in dipolar bose–einstein condensate. Chaos An Interdiscip. J. Nonlinear Sci. 2021, 31, 023120. [Google Scholar] [CrossRef]
  41. Xiong, B.; Gong, J.; Pu, H.; Bao, W.; Li, B. Symmetry breaking and self-trapping of a dipolar Bose-Einstein condensate in a double-well potential. Phys. Rev. A 2009, 79, 013626. [Google Scholar] [CrossRef] [Green Version]
  42. Abad, M.; Guilleumas, M.; Mayol, R.; Pi, M.; Jezek, D.M. A dipolar self-induced bosonic Josephson junction. Europhys. Lett. 2011, 94, 10004. [Google Scholar] [CrossRef] [Green Version]
  43. Wang, C.; Kevrekidis, P.G.; Frantzeskakis, D.J.; Malomed, B.A. Effects of long-range nonlinear interactions in double-well potentials. Phys. D Nonlinear Phenom. 2011, 240, 805–813. [Google Scholar] [CrossRef] [Green Version]
  44. Adhikari, S.K. Self-trapping of a dipolar Bose-Einstein condensate in a double well. Phys. Rev. A 2014, 89, 043609. [Google Scholar] [CrossRef] [Green Version]
  45. Zhang, A.-X.; Xue, J.-K. Dipolar-induced interplay between inter-level physics and macroscopic phase transitions in triple-well potentials. J. Phys. B At. Mol. Opt. Phys. 2012, 45, 145305. [Google Scholar] [CrossRef]
  46. Fortanier, R.; Zajec, D.; Main, J.; Wunner, G. Dipolar Bose–Einstein condensates in triple-well potentials. J. Phys. B At. Mol. Opt. Phys. 2013, 46, 235301. [Google Scholar] [CrossRef] [Green Version]
  47. Bouchoule, I.; Molmer, K. Spin squeezing of atoms by the dipole interaction in virtually excited Rydberg states. Phys. Rev. A 2002, 65, 041803. [Google Scholar] [CrossRef]
  48. Henkel, N.; Nath, R.; Pohl, T. Three-dimensional roton excitations and supersolid formation in rydberg-excited bose-einstein condensates. Phys. Rev. Lett. 2010, 104, 195302. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Honer, J.; Weimer, H.; Pfau, T.; Büchler, H.P. Collective many-body interaction in rydberg dressed atoms. Phys. Rev. Lett. 2010, 105, 160404. [Google Scholar] [CrossRef] [Green Version]
  50. Pupillo, G.; Micheli, A.; Boninsegni, M.; Lesanovsky, I.; Zoller, P. Strongly correlated gases of rydberg-dressed atoms: Quantum and classical dynamics. Phys. Rev. Lett. 2010, 104, 223002. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  51. Johnson, J.E.; Rolston, S.L. Interactions between Rydberg-dressed atoms. Phys. Rev. A 2010, 82, 033412. [Google Scholar] [CrossRef] [Green Version]
  52. Li, W.; Hamadeh, L.; Lesanovsky, I. Probing the interaction between Rydberg-dressed atoms through interference. Phys. Rev. A 2012, 85, 053615. [Google Scholar] [CrossRef] [Green Version]
  53. DeSalvo, B.J.; Aman, J.A.; Gaul, C.; Pohl, T.; Yoshida, S.; Burgdörfer, J.; Hazzard, K.R.A.; Dunning, F.B.; Killian, T.C. Rydberg-blockade effects in autler-townes spectra of ultracold strontium. Phys. Rev. A 2016, 93, 022709. [Google Scholar] [CrossRef] [Green Version]
  54. Hsueh, C.-H.; Wang, C.-W.; Wu, W.-C. Vortex structures in a rotating rydberg-dressed bose-einstein condensate with the lee-huang-yang correction. Phys. Rev. A 2020, 102, 063307. [Google Scholar] [CrossRef]
  55. Maucher, F.; Henkel, N.; Saffman, M.; Królikowski, W.; Skupin, S.; Pohl, T. Rydberg-induced solitons: Three-dimensional self-trapping of matter waves. Phys. Rev. Lett. 2011, 106, 170401. [Google Scholar] [CrossRef] [PubMed]
  56. Cinti, F.; MacRì, T.; Lechner, W.; Pupillo, G.; Pohl, T. Defect-induced supersolidity with soft-core bosons. Nat. Commun. 2014, 5, 4235. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Hsueh, C.H.; Tsai, Y.C.; Wu, W.C. Excitations of one-dimensional supersolids with optical lattices. Phys. Rev. A 2016, 93, 063605. [Google Scholar] [CrossRef]
  58. McCormack, G.; Nath, R.; Li, W. Dynamical excitation of maxon and roton modes in a Rydberg-Dressed Bose-Einstein Condensate. Phys. Rev. A 2020, 102, 023319. [Google Scholar] [CrossRef]
  59. Lauer, A.; Muth, D.; Fleischhauer, M. Transport-induced melting of crystals of Rydberg dressed atoms in a one-dimensional lattice. New J. Phys. 2012, 14, 095009. [Google Scholar] [CrossRef] [Green Version]
  60. Lan, Z.; Minar, J.; Levi, E.; Li, W.; Lesanovsky, I. Emergent Devil’s Staircase without Particle-Hole Symmetry in Rydberg Quantum Gases with Competing Attractive and Repulsive Interactions. Phys. Rev. Lett. 2015, 115, 203001. [Google Scholar] [CrossRef] [Green Version]
  61. Angelone, A.; Mezzacapo, F.; Pupillo, G. Superglass Phase of Interaction-Blockaded Gases on a Triangular Lattice. Phys. Rev. Lett. 2016, 116, 135303. [Google Scholar] [CrossRef] [Green Version]
  62. Chougale, Y.; Nath, R. Ab initio calculation of Hubbard parameters for Rydberg-dressed atoms in a one-dimensional optical lattice. J. Phys. B At. Mol. Opt. Phys. 2016, 49, 144005. [Google Scholar] [CrossRef]
  63. Li, Y.; Geißler, A.; Hofstetter, W.; Li, W. Supersolidity of lattice bosons immersed in strongly correlated Rydberg dressed atoms. Phys. Rev. A 2018, 97, 023619. [Google Scholar] [CrossRef] [Green Version]
  64. Zhou, Y.; Li, Y.; Nath, R.; Li, W. Quench dynamics of Rydberg-dressed bosons on two-dimensional square lattices. Phys. Rev. A 2020, 101, 013427. [Google Scholar] [CrossRef] [Green Version]
  65. Barbier, M.; Geißler, A.; Hofstetter, W. Decay-dephasing-induced steady states in bosonic rydberg-excited quantum gases in an optical lattice. Phys. Rev. A 2019, 99, 033602. [Google Scholar] [CrossRef] [Green Version]
  66. Jau, Y.Y.; Hankin, A.M.; Keating, T.; Deutsch, I.H.; Biedermann, G.W. Entangling atomic spins with a Rydberg-dressed spin-flip blockade. Nat. Phys. 2016, 12, 3487. [Google Scholar] [CrossRef] [Green Version]
  67. Zeiher, J.; van Bijnen, R.; Schauß, P.; Hild, S.; Choi, J.Y.; Pohl, T.; Bloch, I.; Gross, C. Many-body interferometry of a Rydberg-dressed spin lattice. Nat. Phys. 2016, 12, 3835. [Google Scholar] [CrossRef]
  68. Zeiher, J.; Choi, J.Y.; Rubio-Abadal, A.; Pohl, T.; Van Bijnen, R.; Bloch, I.; Gross, C. Coherent many-body spin dynamics in a long-range interacting Ising chain. Phys. Rev. X 2017, 7, 041063. [Google Scholar] [CrossRef] [Green Version]
  69. Guardado-Sanchez, E.; Spar, B.M.; Schauss, P.; Belyansky, R.; Young, J.T.; Bienias, P.; Gorshkov, A.V.; Iadecola, T.; Bakr, W.S. Quench Dynamics of a Fermi Gas with Strong Long-Range Interactions. arXiv 2020, arXiv:2010.05871. [Google Scholar]
  70. Borish, V.; Marković, O.; Hines, J.A.; Rajagopal, S.V.; Schleier-Smith, M. Transverse-Field Ising Dynamics in a Rydberg-Dressed Atomic Gas. Phys. Rev. Lett. 2020, 124, 063601. [Google Scholar] [CrossRef] [Green Version]
  71. McCormack, G.; Nath, R.; Li, W. Nonlinear dynamics of Rydberg-dressed Bose-Einstein condensates in a triple-well potential. Phys. Rev. A 2020, 102, 063329. [Google Scholar] [CrossRef]
  72. Wolf, A.; Swift, J.B.; Swinney, H.L.; Vastano, J.A. Determining Lyapunov exponents from a time series. Phys. D Nonlinear Phenom. 1985, 16, 285–317. [Google Scholar] [CrossRef] [Green Version]
  73. Andreev, A.V.; Balanov, A.G.; Fromhold, T.M.; Greenaway, M.T.; Hramov, A.E.; Li, W.; Makarov, V.V.; Zagoskin, A.M. Emergence and control of complex behaviors in driven systems of interacting qubits with dissipation. NPJ Quantum Inf. 2021, 7, 1–7. [Google Scholar] [CrossRef]
  74. Eckel, S.; Lee, J.G.; Jendrzejewski, F.; Murray, N.; Clark, C.W.; Lobb, C.J.; Phillips, W.D.; Edwards, M.; Campbell, G.K. Hysteresis in a quantized superfluid ’atomtronic’ circuit. Nature 2014, 506, 200–203. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  75. Trenkwalder, A.; Spagnolli, G.; Semeghini, G.; Coop, S.; Landini, M.; Castilho, P.; Pezzè, L.; Modugno, G.; Inguscio, M.; Smerzi, A.; et al. Quantum phase transitions with parity-symmetry breaking and hysteresis. Nat. Phys. 2016, 12, 826–829. [Google Scholar] [CrossRef] [Green Version]
  76. Bürkle, R.; Vardi, A.; Cohen, D.; Anglin, J.R. Probabilistic hysteresis in integrable and chaotic isolated hamiltonian systems. Phys. Rev. Lett. 2019, 123, 114101. [Google Scholar] [CrossRef] [Green Version]
  77. Scully, M.O.; Zubairy, M.S. Quantum Optics, 1st ed.; Cambridge University Press: Cambridge, UK, 1997. [Google Scholar]
  78. Castro, E.R.; Chávez-Carlos, J.; Roditi, I.; Santos, L.F.; Hirsch, J.G. Quantum-classical correspondence of a system of interacting bosons in a triple-well potential. arXiv 2021, arXiv:2105.10515. [Google Scholar] [CrossRef]
  79. Dey, A.; Cohen, D.; Vardi, A. Adiabatic Passage through Chaos. Phys. Rev. Lett. 2018, 121, 250405. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  80. Dey, A.; Cohen, D.; Vardi, A. Many-body adiabatic passage: Quantum detours around chaos. Phys. Rev. A 2019, 99, 033623. [Google Scholar] [CrossRef] [Green Version]
  81. Datseris, G. Dynamicalsystems.jl: A julia software library for chaos and nonlinear dynamics. J. Open Source Softw. 2018, 3, 598. [Google Scholar] [CrossRef]
  82. Baier, G.; Klein, M. Maximum hyperchaos in generalized Hénon maps. Phys. Lett. A 1990, 151, 281–284. [Google Scholar] [CrossRef]
  83. Baier, G.; Sahle, S. Design of hyperchaotic flows. Phys. Rev. E 1995, 51, R2712–R2714. [Google Scholar] [CrossRef]
  84. Kapitaniak, T.; Thylwe, K.-E.; Cohen, I.; Wojewoda, J. Chaos-hyperchaos transition. Chaos Solitons Fractals 1995, 5, 2003–2011. [Google Scholar] [CrossRef]
  85. Tarkhov, A.E.; Wimberger, S.; Fine, B.V. Extracting lyapunov exponents from the echo dynamics of bose-einstein condensates on a lattice. Phys. Rev. A 2017, 96, 023624. [Google Scholar] [CrossRef] [Green Version]
  86. Lerose, A.; Pappalardi, S. Bridging entanglement dynamics and chaos in semiclassical systems. Phys. Rev. A 2020, 102, 032404. [Google Scholar] [CrossRef]
  87. Pausch, L.; Carnio, E.G.; Rodríguez, A.; Buchleitner, A. Chaos and Ergodicity across the Energy Spectrum of Interacting Bosons. Phys. Rev. Lett. 2021, 126, 150601. [Google Scholar] [CrossRef]
  88. Kollath, C.; Roux, G.; Biroli, G.; Läuchli, A.M. Statistical properties of the spectrum of the extended Bose-Hubbard model. J. Stat. Mech. 2010, 2010, P08011. [Google Scholar] [CrossRef]
  89. Chen, Y.; Cai, Z. Persistent oscillations versus thermalization in the quench dynamics of quantum gases with long-range interactions. Phys. Rev. A 2020, 101, 023611. [Google Scholar] [CrossRef] [Green Version]
Figure 1. (Color online) The extended Bose-Hubbard Chain and quenching schemes. (a) Nearest-neighbor (U) and next-nearest-neighbor (V) interactions between atoms in a one-dimensional optical lattice (lattice constant d). The tilting of the lattice is denoted by the parameter γ . We consider a linear quench in (b) γ and (c) U towards a non-zero value (solid). When γ (U) returns to the initial value (both solid and dashed), this is a hysteresis quench. The rate to quench γ (U) is α ( β ). See text for details of the soft-core interaction and quenching protocols.
Figure 1. (Color online) The extended Bose-Hubbard Chain and quenching schemes. (a) Nearest-neighbor (U) and next-nearest-neighbor (V) interactions between atoms in a one-dimensional optical lattice (lattice constant d). The tilting of the lattice is denoted by the parameter γ . We consider a linear quench in (b) γ and (c) U towards a non-zero value (solid). When γ (U) returns to the initial value (both solid and dashed), this is a hysteresis quench. The rate to quench γ (U) is α ( β ). See text for details of the soft-core interaction and quenching protocols.
Photonics 08 00554 g001
Figure 2. Eigenenergies, Bogoliubov spectra, and Lyapunov exponents when varying the tilt γ . We show (a) the nonlinear eigenenergy, the Bogoliubov spectra of (b) the ground state and (c) the first excited state, and (d) the Lyapunov exponents of the ground state. The nonlinearity dominates when | γ | is small, leading to loops in the eigenenergy. The Bogoliubov spectra are all real when the system is in the ground state (b). The Bogoliubov spectra have complex components (red region) when the system is in the first-excited state. Positive Lyapunov exponents indicate the system exhibits chaos dynamically, which appear mostly in the loop region of the eigenenergy. Parameters are L = 3 and U = 2 V = 5 .
Figure 2. Eigenenergies, Bogoliubov spectra, and Lyapunov exponents when varying the tilt γ . We show (a) the nonlinear eigenenergy, the Bogoliubov spectra of (b) the ground state and (c) the first excited state, and (d) the Lyapunov exponents of the ground state. The nonlinearity dominates when | γ | is small, leading to loops in the eigenenergy. The Bogoliubov spectra are all real when the system is in the ground state (b). The Bogoliubov spectra have complex components (red region) when the system is in the first-excited state. Positive Lyapunov exponents indicate the system exhibits chaos dynamically, which appear mostly in the loop region of the eigenenergy. Parameters are L = 3 and U = 2 V = 5 .
Photonics 08 00554 g002
Figure 3. Eigenenergies and Lyapunov exponents as a function of U. We show eigenenergy for (a) L = 3 and (b) L = 5 when the trap is balanced ( γ = 0 ). Level crossings are found when the interaction is strong. Starting from the ground state, we calculate Lyapunov exponents for (c) L = 3 and (d) L = 5 . For a given U, Lyapunov exponents of same value but opposite signs appear in pairs.
Figure 3. Eigenenergies and Lyapunov exponents as a function of U. We show eigenenergy for (a) L = 3 and (b) L = 5 when the trap is balanced ( γ = 0 ). Level crossings are found when the interaction is strong. Starting from the ground state, we calculate Lyapunov exponents for (c) L = 3 and (d) L = 5 . For a given U, Lyapunov exponents of same value but opposite signs appear in pairs.
Photonics 08 00554 g003
Figure 4. (Color online) Final population distribution of the ground state. The population by quenching γ with (a) scheme I and (b) scheme II is shown for L = 3 . In the numerical simulation, γ i = γ f = 10 and the interaction strength is U = 5 . The interaction U is quenched with (c) scheme III and (d) scheme IV , where U i = 0 , U f = 10 and γ = 0 , respectively. In all the figures, the quench rates ( α or β ) are 1 (blue), 0.1 (green), and 0.01 (red). The total number of trajectories is M = 100 . The target final state is shown as the large black circle.
Figure 4. (Color online) Final population distribution of the ground state. The population by quenching γ with (a) scheme I and (b) scheme II is shown for L = 3 . In the numerical simulation, γ i = γ f = 10 and the interaction strength is U = 5 . The interaction U is quenched with (c) scheme III and (d) scheme IV , where U i = 0 , U f = 10 and γ = 0 , respectively. In all the figures, the quench rates ( α or β ) are 1 (blue), 0.1 (green), and 0.01 (red). The total number of trajectories is M = 100 . The target final state is shown as the large black circle.
Photonics 08 00554 g004
Figure 5. (Color online) Bogoliubov spectra and Lyapunov exponents of the localized state. Dynamically unstable regions (dark red) for (a) L = 3 and (b) L = 5 are shown as a function of U and γ . Panels (c,d) give the Lyapunov exponents as a function of U. Random perturbation to the initial state are examined for (e) L = 3 and (f) L = 5 . The red lines show the maximal Lyapunov exponents in (c,d), correspondingly. Here, γ = 0 in panels (c,d).
Figure 5. (Color online) Bogoliubov spectra and Lyapunov exponents of the localized state. Dynamically unstable regions (dark red) for (a) L = 3 and (b) L = 5 are shown as a function of U and γ . Panels (c,d) give the Lyapunov exponents as a function of U. Random perturbation to the initial state are examined for (e) L = 3 and (f) L = 5 . The red lines show the maximal Lyapunov exponents in (c,d), correspondingly. Here, γ = 0 in panels (c,d).
Photonics 08 00554 g005
Figure 6. (Color online) Final population distribution of the localized state. The first and second row show the linear and hysteresis quench of γ . Here, U = 5 , γ i = γ f = 10 . The third and fourth row show the linear and hysteresis quench of U, with U i = 0 , U f = 10 , and γ = 0 . In (ad), we consider three sites and the initial thermal state is Ψ ¯ = [ 0.1 e i ϕ 1 , 0.98 e i ϕ 2 , 0.1 e ϕ 3 ] , with ϕ j ( j = 1 , 2 , 3 ) being a random number in [ 0 , 2 π ] . In (eh), L = 5 and the initial state is Ψ ¯ = [ 0.005 e i ϕ 1 , 0.98 e i ϕ 2 , 0.005 e i ϕ 3 , 0.005 e i ϕ 4 , 0.005 i ϕ 5 ] , with ϕ j ( j = 1 , , 5 ) being randomly distributed in [ 0 , 2 π ] . The small fraction in sites other than the localized state is used to trigger the hopping dynamics. The other parameters are the same as the ones in Figure 4.
Figure 6. (Color online) Final population distribution of the localized state. The first and second row show the linear and hysteresis quench of γ . Here, U = 5 , γ i = γ f = 10 . The third and fourth row show the linear and hysteresis quench of U, with U i = 0 , U f = 10 , and γ = 0 . In (ad), we consider three sites and the initial thermal state is Ψ ¯ = [ 0.1 e i ϕ 1 , 0.98 e i ϕ 2 , 0.1 e ϕ 3 ] , with ϕ j ( j = 1 , 2 , 3 ) being a random number in [ 0 , 2 π ] . In (eh), L = 5 and the initial state is Ψ ¯ = [ 0.005 e i ϕ 1 , 0.98 e i ϕ 2 , 0.005 e i ϕ 3 , 0.005 e i ϕ 4 , 0.005 i ϕ 5 ] , with ϕ j ( j = 1 , , 5 ) being randomly distributed in [ 0 , 2 π ] . The small fraction in sites other than the localized state is used to trigger the hopping dynamics. The other parameters are the same as the ones in Figure 4.
Photonics 08 00554 g006
Figure 7. (Color online) Lyapunov Exponents vs. System Size. The maximum Lyapunov exponent and total number of positive Lyapunov exponents are shown in (a,b) for the ground state configuration. Panels (c,d) show the same quantity for the localized state. The larger the value of U is, the larger the maximal Lyapunov exponent. The maximal Lyapunov exponent decreases with increasing L. The number of Lyapunov exponents increases and then decreases with increasing L. For the localized state, η increases almost linearly with increasing L. In each panel, U = 1 (square), 3 (circle), and 5 (triangle).
Figure 7. (Color online) Lyapunov Exponents vs. System Size. The maximum Lyapunov exponent and total number of positive Lyapunov exponents are shown in (a,b) for the ground state configuration. Panels (c,d) show the same quantity for the localized state. The larger the value of U is, the larger the maximal Lyapunov exponent. The maximal Lyapunov exponent decreases with increasing L. The number of Lyapunov exponents increases and then decreases with increasing L. For the localized state, η increases almost linearly with increasing L. In each panel, U = 1 (square), 3 (circle), and 5 (triangle).
Photonics 08 00554 g007
Figure 8. (Color online) Poincaré Sections of the ground state and localized state on the U 1 -plane. The Poincaré sections are shown for the ground state (ac) and the localized state (df). Each point represents a numerical realization. We consider L = 5 (a,d), L = 10 (b,e), and L = 20 (c,f). Other parameters are U = 3 and γ = 0 .
Figure 8. (Color online) Poincaré Sections of the ground state and localized state on the U 1 -plane. The Poincaré sections are shown for the ground state (ac) and the localized state (df). Each point represents a numerical realization. We consider L = 5 (a,d), L = 10 (b,e), and L = 20 (c,f). Other parameters are U = 3 and γ = 0 .
Photonics 08 00554 g008
Figure 9. (Color online) Areas of the Poincaré Sections. We compare the fitted area (open shapes) of the Poincaré section with λ m (solid) for both the ground state (a) and the localized state (b), respectively. The blue circles are for U = 3 , and red triangles for U = 5 . In both situations, γ = 0 .
Figure 9. (Color online) Areas of the Poincaré Sections. We compare the fitted area (open shapes) of the Poincaré section with λ m (solid) for both the ground state (a) and the localized state (b), respectively. The blue circles are for U = 3 , and red triangles for U = 5 . In both situations, γ = 0 .
Photonics 08 00554 g009
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

McCormack, G.; Nath, R.; Li, W. Hyperchaos in a Bose-Hubbard Chain with Rydberg-Dressed Interactions. Photonics 2021, 8, 554. https://doi.org/10.3390/photonics8120554

AMA Style

McCormack G, Nath R, Li W. Hyperchaos in a Bose-Hubbard Chain with Rydberg-Dressed Interactions. Photonics. 2021; 8(12):554. https://doi.org/10.3390/photonics8120554

Chicago/Turabian Style

McCormack, Gary, Rejish Nath, and Weibin Li. 2021. "Hyperchaos in a Bose-Hubbard Chain with Rydberg-Dressed Interactions" Photonics 8, no. 12: 554. https://doi.org/10.3390/photonics8120554

APA Style

McCormack, G., Nath, R., & Li, W. (2021). Hyperchaos in a Bose-Hubbard Chain with Rydberg-Dressed Interactions. Photonics, 8(12), 554. https://doi.org/10.3390/photonics8120554

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop