Hyperchaos in a Bose-Hubbard chain with Rydberg-dressed interactions

We study 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 groundstate and localized state configuration. 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\'e 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, and collective and emergent nonlinear dynamics of BECs with long-range interactions.


I. INTRODUCTION
Over the past two decades, Bose-Einstein condensates (BECs) of ultracold atomic gases have become an ideal system to study both quantum and nonlinear dynamics, due to the high controllability over 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, 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 in a single site.This is in contrast to the homogeneous, superfluid state, which form the groundstate 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 meanfield 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 bifur-< l a t e x i t s h a 1 _ b a s e 6 4 = " S K m q a w Q O 3 p E I E O 2 Z q Q 0 X j V i 3 N d E = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K q M e i F 4 8 t m F p o Q 9 l s J + 3 a z S b s b o R S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 M B V c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o p Z N M M f R Z I h L V D q l G w S X 6 h h u B 7 V Q h j U O B D + H o d u Y / P K H S P J H 3 Z p x i E N O B 5 B F n 1 F i p 6 f f K F b f q z k F W i Z e T C u R o 9 M p f 3 X 7 C s h i l Y Y J q 3 f H c 1 A Q T q g x n A q e l b q Y x p W x E B 9 i x V N I Y d T C Z H z o l Z 1 b p k y h R t q Q h c / X 3 x I T G W o / j 0 H b G 1 A z 1 s j c T / / M 6 m Y m u g w m X a W Z Q s s W i K B P E J G T 2 N e l z h c y I s S W U K W 5 v J W x I F W X G Z l O y I X j L L 6 + S V q 3 q X V Z r z Y t K / S a P o w g n c A r n 4 M E V 1 O E O G u A D A 4 R n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f s 3 e M 4 A = = < / l a t e x i t > U < l a t e x i t s h a 1 _ b a s e 6 4 = " l d E C R D l K J 7 k L g y 8 T 4 j g C w 2 M X W 6 8 = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m K q M e i F 4    cations of the eigenspectra and chaotic dynamics, when confined in harmonic trap [39,40].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 groundstate 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 Fig. 1(a).Various theoretical studies on the static and dynamical properties of Rydberg-dressed atoms con-arXiv:2108.09683v2[cond-mat.quant-gas]24 Aug 2021 fined 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].In Ref. [71] we have shown that self-trapping dynamics of Rydberg-dressed BECs can be controlled in a triple-well potential through mean-field and quantum mechanical analysis.
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 dynamical stability of the groundstate and localized state, where dependence of the largest, and total number of positive Lyapunov exponents [72,73] on the dressed interaction and system size is explored.We probe the chaotic dynamics by employing both linear and a hysteresis quench of the potential bias and dressed interaction [74][75][76].
The paper is organized as follows.In Sec.II 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 static (eigenenergies and Bogoliubov spectra) and dynamical properties (Lyapunov exponents) of the groundstate and localized state configurations in Sec.III, and Sec.IV, respectively.Dynamics driven by both the linear and hysteresis quenching parameters are explored with different initial states.In Sec.V 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 proportional almost linearly to the largest Lyapunov exponent.We conclude our work in Sec.VI.

II. MODEL AND METHOD
A. Extended Bose-Hubbard model in the semiclassical limit Our setting consists of N bosonic atoms confined in a one-dimensional lattice with lattice constant d, as depicted in Fig. 1(a).The Rydberg-dressing induces longrange interactions between atoms at different sites.Taking into account of hopping between nearest-neighbor sites, we obtain an extended Bose-Hubbard Hamiltonian of L sites [52] where âj (â † 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, nj = â † j â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 Ĥint = g 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 Λ 1(a)].Both the soft-core radius R and C 6 can be tuned by laser parameters [48].In this work, we will restrict to the onsite, nearest-neighbor (Λ j,j±1 ) and nextnearest-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 meanfield approximation where the bosonic operator is described by a classical field ψ j , i.e. âj ≈ ψ j √ N , and â † This yields the semiclassical Hamiltonian H ≈ Ĥ/N , The dynamics of the classical field ψ j is obtained via the canonical equation idψ j /dt = ∂H/∂ψ * j , yielding the coupled GP equations 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 strength.The onsite interaction W takes into account contributions from both the s-wave and soft-core interaction.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 interaction to be U = 2V 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 ( 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 â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].

B. Nonlinear eigenenergies and Bogoliubov spectra
Though the Hamiltonian ( 2) is Hermitian, the densitydependent 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 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 is obtained selfconsistently.
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], where L = H0 + 2U P − µ, and N = −U P. H0 and P are L × L block matrices.From Eq. (3), we obtain the matrix elements ψ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, whose rate is determined by the imaginary part of the spectra.

C. 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 Due to the normalization condition, we need to solve a 2L − 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 Uk > 0. These intersecting points form the 2D Poincaré section.To be specific, we will evaluate the Poincaré section of variable (R 2 , I 2 ) on the U 1 plane.
The strength of 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 will consider both the groundstate and a localized state initially.In a localized state, nearly all the condensate sits in a single site, which can be stable (i.e. the selftrapping state) when the nonlinear interaction is strong.In this work, the Lyapunov exponents λ j (j = 1, • • • , 2L) are calculated via DynamicalSystems.jl, a fast and reliable Julia library to determine the dynamics of nonlinear systems [81].We have checked that it gives consistent data with the method in Ref. [72].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].

D. Quenching schemes
In Sections.III and IV 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 function 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 Fig. 1(b).
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 The corresponding scheme is shown by the solid and dashed curve in Fig. 1(b).
Scheme III : In addition to quenching the level bias we also change the two-body interaction strength through a linear ramp, where U i is the initial interaction strength and β is the quench rate.This is shown by the solid curve in Fig. 1(c).Note that the next-nearest-neighbor interaction V depends on time as well due to the relation U = 2V .Scheme IV : The hysteresis counterpart of the interaction quench is given by where U f is the final interaction strength, with τ = U f /β.

III. STABILITY OF THE GROUNDSTATE A. 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 Fig. 2(a).We find 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 2L 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 groundstate, i.e., the Bogoliubov spectra B are real, as shown in Fig. 2(b).This is in contrast to excited eigenstates, whose Bogoliubov spectra have imaginary components.As an example, the Bogoliubov spectra of the first excited state is shown in Fig. 2(c).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 Fig. 2(d) Lyapunov exponents are shown for the groundstate 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 groundstate could gain exponential growth, and hence drives the system away from the groundstate.
To further understand roles played by the nonlinearity, we calculate eigenenergies as a function of the interaction strength U shown in Fig. 3(a) and (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 groundstate of the nonlinear system are shown in Figs.3(c) and 3(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.Importantly, the number of Lyapunov exponents now depends on L. For L = 3, one obtains single positive Lyapunov exponent when U 4. When L = 5, there are 3 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 is conserved in the Bose-Hubbard chain.

B. 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 quenching the tilt potential and interaction strength given by Eqs. ( 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 groundstate), 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 groundstate 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 [Fig.4(a)].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.Increasing U , the population is distributed into a larger region of phase space, i.e., it occupies a 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 Eq. ( 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 [Fig.4(b)].Here the hysteresis quench has allowed for a large level of reversibility in the dynamics [76], as the chaotic regions have not been triggered.Increasing the quench rate α, the population distributions cluster around much smaller regions in phase space, than the one shown in panel (a).
In Fig. 4(c) we quench the interaction according to Eq. ( 9).The initial states depend on the value of γ.For example the groundstate 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 Fig. 3 (a) and Lyapunov exponent in Fig. 3(c)].We have verified that by increasing γ the associated randomness with the final states decreases, as the number of crossings in the eigenenergy will decrease.
In case of the hysteresis quench of U , we find that the results [Fig.4(d)] are similar to the linear quench.When looking at γ = 0, the final states do not return to the initial value.As shown in Fig. 3(c), the Lyapunov exponent of the groundstate becomes positive when U 4, which causes the final state more random, i.e. a broader distribution of the densities.As the tilt γ increases, we have verified that chaos is gradually suppressed, as the population localizes in the trap corresponds to the lowest energy state.In order to trigger chaotic dynamics in the tilted case, stronger interactions are needed in general.

A. 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 groundstate 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 to be Ψ ] with 1 and φ j to be a random phase.This choice furthermore insures that the energy of different initial states are almost identical.
In Figs.5(a) and (b), dynamical unstable regions in the Bogoliubov spectra for L = 3 and L = 5 are shown (highlighted with dark red color).In the unstable region, B develops imaginary components, which depend on U , γ and L. In case of L = 3, the condensate is localized in the middle site initially, meaning the Bogoliubov spectra are symmetric with respect to γ. Fig. 5(a) shows that the system is dynamically unstable when U is small, in par- The small fraction in sites other than the localized state is used to trigger the hopping dynamics.Other parameters are same with the one in Fig. 4.
ticular when the lattice is balanced (|γ| is small).This is not surprising, as the localized state is not the groundstate, nor the system supports 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 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 Fig. 5(c) the Lyapunov exponents for L = 3 show an unusually symmetric shape when U 4. The exponents are a smooth function of U , and reaches maximal value around U = 5.Further increasing U , the positive Lyapunov exponents decrease.This indicates that the localized configuration could exhibit chaotic dynamics for 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 [Fig.5(d)], 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 wavefunction of the traps on either side of the localized site, with 0 < ε < 0.01.In Fig. 5(e) and (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.

B. Quench dynamics
For U = 5 and L = 3, a linear quench [Fig.6(a)] 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 Fig. 6(b) 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 Eqs. ( 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 manor 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 a complete localization of the condensate on 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 Fig. 6(e) 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 site, implying 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 lost completely, and is replaced with a chaotic distributions, largely due to the presence of hyperchaos [see Fig. 5(d)].

V. SCALING OF LYAPUNOV EXPONENTS WITH THE SYSTEM SIZE
In the following we will investigate 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 can not 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 groundstate and the localized state, will be examined in detail.
In Fig. 7(a) the largest Lyapunov exponent λ m for the groundstate configuration is shown.When 2 ≤ L ≤ 4, the values of λ m are small in general.This is due to the fact that chaos has not be triggered [see Fig. 5(b) and (c)].When L > 4, the situation changes as chaos is already found with the given U .We find λ 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 Fig. 7(b).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 .Recently it has been shown that the largest Lyapunov exponents in the BH model can be obtained from the echo dynamics of the condensate [85].Similar technique could be applied to extract the largest Lyapunov exponents studied here.
Figs. 7(c) and (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 groundstate, 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 deviates 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 state, η > 1 can be seen even with relatively weak interaction (e.g.U = 1), leading to more pronounced hyperchaotic dynamics.
The total number of nonlinear differential equations is 2L (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 Fig. 2(d ].As the extended Bose-Hubbard model is a Hamiltonian system, not only the sum of the Lyapunov exponents vanishes, but also conservative quantities, such as the energy and particle number, are found in the dynamics.This indicates that the maximal number of the Lyapunov exponents is L − 2 but not L.For sufficiently large L, the total number of positive Lyapunov exponents is smaller than L − 2, as the nonlinear interaction becomes smaller.For the groundstate, one can estimate the interaction energy for a given site to be 2(U + V )/L 2 approximately, 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 as 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 Fig. 8, showing that profiles of the Poincaré section depend on the system size and the initial state.When L = 5 the area is largest [Fig.8(a)] and decrease with increasing L in case of the groundstate [Fig.8(b) and (c)].For different L, the profile of the Poincaré section is largely symmetric with respect to R 2 = 0 and I 2 = 0.In case of the localized state, similar dependence on L is found, as depicted in Fig. 8(d)-(f).We note two differences compared to the groundstate ones.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 Fig.The area is largely determined by the largest Lyapunov exponent.To verify this, we find the area of the Poincaré section approximately through numerically fitting the Poincaré section, shown in Fig. 9.For the groundstate, 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 groundstate could be a boundary effect when L is small, as the localized state suffers less from the boundary effect.

VI. CONCLUSION AND OUTLOOK
We have investigated the chaotic and hyperchaotic dynamics of a one-dimensional Bose-Hubbard chain of Rydberg-dressed BECs in the semiclassical regime.We have shown that both the groundstate and localized state can have positive Lyapunov exponents, even though the corresponding Bogoliubov spectra are real valued.As a result, small perturbations to these states lead to large fluctuations, which are captured by the quench dynamics.We have found that hyperchaos emerges in both the groundstate 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 have shown that η grows with the system size L when U is large.So far our investigations are focusing 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 eigenspectra on the Bose-Hubbard model with onsite [87] and long-range interactions [88,89].It is therefore worthwhile to explore features of chaos and hyperchaos due to the Rydberg dressed interaction in the quantum regime.
8 t 2 F p o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W yx W M S q E 1 C N g k t s G W 4 E d h K F N A o E P g T j 2 5 n / 8 I R K 8 1 j e m 0 m C f k S H k o e c U W O l Z r t f r r h V d w 6 y S r y c V C B H o 1 / + 6 g 1 i l k Y o D R N U 6 6 7 n J s b P q D K c C Z y W e q n G h L I x H W L X U k k j 1 H 4 2 P 3 R K z q w y I G G s b E l D 5 u r v i Y x G W k + i w H Z G 1 I z 0 s j c T / / O 6 q Q m v / Y z L J D U o 2 W J R m A p i Y j L 7 m g y 4 Q m b E x B L K F L e 3 E j a i i j J j s y n Z E L z l l 1 d J u 1 b 1 L q u 1 5 k W l f p P H U Y Q T O I V z 8 O A K 6 n A H D W g B A 4R n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f t P u M 4 Q = = < / l a t e x i t > V < l a t e x i t s h a 1 _ b a s e 6 4 = " c 8 + X q F a D W 2 u s 5 W d j S H 5 R n + b 0 j n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f y j O M 7 w = = < / l a t e x i t > d < l a t e x i t s h a 1 _ b a s e 6 4 = " X j e C 9 a H / I u m e Z Y U H z h y x A S l F g 9 M = " > A A A B 6 3 i c b V B N S 8 N A E J 3 4 W e t X 1 a O X x S J 4 K k k R 9 V j 0 4 r G C / Y A 2 l M 1 2 0 y 7 d 3 Y T d i V B C / 4 I X D 4 p 4 9 Q 9 5 8 9 + Y t D l o 6 4 O B x 3 s z z M w L Y i k s u u 6 3 s 7 a + s b m 1 X d o p 7 + 7 t H x x W j o 7 b N k o M 4 y 0 W y c h 0 A 2 q 5 F J q 3 U K D k 3 d h w q g L J O 8 H k L v c 7 T 9 x Y E e l H n M b c V 3 S k R S g Y x V z q I 0 0 G l a p b c + c g q 8 Q r S B U K N A e V r / 4 w Y o n i G p m k 1 v Y 8 N 0 Y / p Q Y F k 3 x W 7 i e W x 5 R N 6 I j 3 M q q p 4 t Z P 5 7 f O y H m m D E k Y m a w 0 k r n 6 e y K l y t q p C r J O R X F s l 7 1 c / M / r J R j e + K n Q c Y J c s 8 W i M J E E I 5 I / T o b C c I Z y m h H K j M h u J W x M D W W Y x V P O Q v C W X 1 4 l 7 X r N u 6 r V H y 6 r j d s i j h K c w h l c g A f X 0 I B 7 a E I L G I z h G V 7 h z V H O i / P u f C x a 1 5 x i 5 g T + w P n 8 A S N 3 j k 8 = < / l a t e x i t > ⌧ < l a t e x i t s h a 1 _ b a s e 6 4 = " N 2 9 t S f 3 k w FIG. 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.

FIG. 2 .
FIG.2.Eigenenergies, Bogoliubov spectra and Lyapunov exponents when varying the tilt γ.We show (a) the nonlinear eigenenergy, Bogoliubov spectra of (b) the groundstate and (c) the first excited state , and (d) the Lyapunov exponents of the groundstate.The nonlinearity dominates when |γ| is small, leading to loops in the eigenenergy.The Bogoliubov spectra are all real when the system is in the groundstate (b).The Bogoliubov spectra have complex components (red region) when the system is in the firstexcited 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 = 2V = 5.

FIG. 3 .
FIG. 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 groundstate, 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.

FIG. 4 .
FIG. 4. (color online) Final population distribution of the groundstate.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 IVwhere Ui = 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.

FIG. 5 .
FIG. 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) and (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) and (d), correspondingly.Here γ = 0 in panels (c)-(d).

FIG. 7 .
FIG. 7. (color online) Lyapunov Exponents vs System Size.The maximum Lyapunov exponent and total number of positive Lyapunov exponents are shown in (a) and (b) for the groundstate configuration.Panels (c) and (d) show the same quantity for the localized state.The larger 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 nearly linearly with increasing L. In each panel, U = 1 (square), 3 (circle) and 5 (triangle).

FIG. 8 .
FIG. 8. (color online) Poincaré Sections of the groundstate and localized state on the U1-plane.The Poincaré sections are shown for the groundstate (a-c) and the localized state (d-f).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.

FIG. 9 .
FIG. 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 groundstate (a) and localized state (b), respectively.The blue circles are for U = 3, and red triangles for U = 5.In both situations γ = 0.