Bosonic quantum dynamics following a linear interaction quench in finite optical lattices of unit filling
aa r X i v : . [ qu a n t - ph ] N ov Bosonic quantum dynamics following a linear interaction quenchin finite optical lattices of unit filling
S.I. Mistakidis, G.M. Koutentakis,
1, 2 and P. Schmelcher
1, 2 Zentrum f¨ur Optische Quantentechnologien, Universit¨at Hamburg,Luruper Chaussee 149, 22761 Hamburg, Germany The Hamburg Centre for Ultrafast Imaging, Universit¨at Hamburg,Luruper Chaussee 149, 22761 Hamburg, Germany (Dated: November 20, 2017)The nonequilibrium ultracold bosonic quantum dynamics in finite optical lattices of unit fillingfollowing a linear interaction quench from a superfluid to a Mott insulator state and vice versa isinvestigated. The resulting dynamical response consists of various inter and intraband tunnelingmodes. We find that the competition between the quench rate and the interparticle repulsion leadsto a resonant dynamical response, at moderate ramp times, being related to avoided crossings inthe many-body eigenspectrum with varying interaction strength. Crossing the regime of weak tostrong interactions several transport pathways are excited. The higher-band excitation dynamicsis shown to obey an exponential decay possessing two distinct time scales with varying ramptime. Studying the crossover from shallow to deep lattices we find that for a diabatic quenchthe excited band fraction decreases, while approaching the adiabatic limit it exhibits a non-linearbehavior for increasing height of the potential barrier. The inverse ramping process from strongto weak interactions leads to a melting of the Mott insulator and possesses negligible higher-bandexcitations which follow an exponential decay for decreasing quench rate. Finally, independentlyof the direction that the phase boundary is crossed, we observe a significant enhancement of theexcited to higher-band fraction for increasing system size.Keywords: nonequilibrium dynamics; linear interaction quench; superfluid to Mott insulator phasetransition; interband tunneling; intraband tunneling.
I. INTRODUCTION
During the past two decades, ultracold atoms in opti-cal lattices emerged as a versatile system to investigatemany-body (MB) phenomena [1–4]. A prominent exam-ple is the experimental observation of the superfluid (SF)to Mott insulator (MI) quantum phase transition [5, 6]which, among others, demonstrated a pure realization ofthe Bose-Hubbard model [7, 8]. Moreover, lattice sys-tems constitute ideal candidates for studying nonequilib-rium quantum phase transitions [9–14], where a numberof defects, induced by time-dependent quenches [15, 16],appear in the time evolving state. The Kibble-Zurekmechanism of such defect formation [15, 17–20], origi-nally addressed in the context of classical phase transi-tions [21, 22], has been tested in different ultracold MBsettings [23–28] and refers to the rate of topological defectformation induced by quenches across phase transitions.Quench dynamics of ultracold bosons confined in anoptical lattice covering the MI-SF transition in both di-rections has been vastly used to examine both the Kibble-Zurek mechanism [19, 25, 26, 29–32] and the approach tothe adiabatic response limit [29, 33–43, 45]. Referringto the Kibble-Zurek mechanism, recent studies [30, 31]of a linear quench across the MI-SF and back in theone-dimensional Bose-Hubbard model demonstrated thatthe excitation density and the correlation length satisfythe Kibble-Zurek scaling for limited ranges of quenchrates. On the other hand and focussing on the slowquench dynamics across the MI-SF transition many re- cent works evinced the formation and melting of Mottdomains [34, 37–40], the growth of interparticle correla-tions [32, 41–44] and the consequent equilibration process[16, 35, 45].Despite the enormous theoretical and experimental ef-forts in this field, the response of such systems subjectedto abrupt or quasi-adiabatic quenches has not been com-pletely understood and deserves further investigation. Inparticular, a nonadiabatic quench inevitably excites thesystem, and a number of defects including higher-bandexcitations can be formed during the dynamics. The lat-ter implies the necessity to consider a multiband treat-ment [46] of the nonequilibrium correlated dynamics andto obtain information about the higher-band excitationspectrum [47–52] being inaccessible by the lowest Bose-Hubbard model or mean-field (MF) methods. Promisingcandidates for such investigations constitute few-bodysystems [53, 54] being accessible by current state of theart experiments [55–58]. In this context, it is possibleto track the microscopic quantum mechanisms [59–62]consisting of intraband and interband tunneling, namelytunneling within the same or between energetically dif-ferent single-particle bands respectively, and to avoid fi-nite temperature effects. Such few-body systems do notserve as a platform to confirm the Kibble-Zurek scalinghypothesis due to their finite size [43]. However, theyprovide useful insights into the largely unexplored scal-ing of few-body defect density including the formationand melting of Mott domains and the excited to higher-band fraction participating in the dynamics.In the present work we consider few bosons confinedin an optical lattice of unit filling. Thereby, the groundstate for increasing interaction strength experiences thefew-body analogue of the SF to MI transition. Wefirst analyze the MB eigenspectrum for varying inter-particle repulsion, revealing the existence of narrow andwide avoided crossings between states of the zeroth andfirst excited band. Then, we apply a linear interactionquench (LIQ) protocol either from weak to strong in-teractions (positive LIQ) or inverserly (negative LIQ)covering in both cases the diabatic to nearly adiabaticcrossing regimes. As a consequence we observe a dy-namical response consisting of the lowest band tunnel-ing and higher-band excitations. Overall, we find anenhanced dynamical response at moderate quench ratesrather than in the abrupt or almost adiabatic regimes.The lowest band dynamics consists of first and secondorder tunneling [63–65]. These modes can be furthermanipulated by tuning either the interaction strengthafter the quench (postquench interaction) or the heightof the potential barriers in the optical lattice. Further-more, we show that following a positive LIQ the excitedto higher-band fraction obeys a bi-exponential decay forvarying ramp time. The latter decay law possesses twotime scales being related to the width of the existingavoided crossings in the eigenspectrum. However, theinterband tunneling [66, 67], with varying height of thepotential barrier exhibits a more complex behavior. Fordiabatic quenches it decreases, while for smaller quenchrates it scales non-linearly possessing a maximum at acertain height of the potential barrier. The latter be-havior manifests the strong dependence of the excitedto higher-band fraction on the quench rate. Moreover,the excited fraction for a varying postquench interactionstrength features different scaling laws. Approaching theregion of the corresponding avoided crossing it exhibitsa non-linear growth, while for stronger interactions it in-creases almost linearly. On the contrary, for a negativeLIQ we observe the melting of the MI. Here, the lowestband transport (intraband tunneling) is reduced whencompared to the inverse scenario, while the excited frac-tion is negligible obeying an exponential decay both withvarying ramp time and potential height. Finally, for bothpositive and negative LIQs the higher-band fraction issignificantly enhanced for increasing system size.This work is organized as follows. In Sec. II we in-troduce our setup and outline the multiband expansionbeing used for a microscopic analysis of the dynamics.Sec. III presents the resulting dynamics induced by a LIQconnecting the weakly to strongly correlated regimes andback in a triple well of unit filling. To extend our find-ings in Sec. IV we discuss the LIQ dynamics for largerlattice systems of unit filling. We summarize and discussfuture perspectives in Sec. V. Appendix A describes ourcomputational methodology.
II. SETUP AND ANALYSIS TOOLS
The Hamiltonian of N identical bosons each of mass M confined in an one-dimensional m -well optical latticeemploying a LIQ protocol reads H = N X i =1 (cid:18) p i M + V sin ( kx i ) (cid:19) + g ( t, τ ) X i To analyze the LIQ induced dynamics of our system, itis instructive first to demonstrate the dependence of theeigenstates on the interparticle repulsion. Thus, we firstinvestigate the eigenspectra of the system with varyinginteraction strength [Sec. III A], which are subsequentlyrelated to the LIQ dynamics [Secs. III B and III C]. A. Eigenspectra In contrast to the discrete Bose-Hubbard model, herewe employ a continuum space Hamiltonian [see Eq. (1)],which allows us to resolve quench induced higher-bandexcitations [75]. For completeness we note that the Bose-Hubbard model is adequate for the theoretical descrip-tion of the quench dynamics in deep lattices (i.e. large V ) and for relatively small quench amplitudes when com-pared to the band gap. To expose the underlying physicalprocesses that lead to the emergence of such MB excitedstates [49–52] we examine, below, the eigenenergies ofthree particles confined in a triple well potential as afunction of the interaction strength g . FIG. 1. Dependence of the lowest 25 eigenenergies E i on theinteraction strength g . The system consists of three bosonsconfined in a triple well with potential depth ( a ) V = 4 and( b ) V = 10. Solid (dashed) lines represent parity even (odd) E i ’s. Wide (narrow) avoided crossings possessing a width δE > . 01 ( δE < . 01) are denoted by solid (dashed) circles.The eigenenergies of the eigenstates that do not contributeto any wide avoided crossing are shown in grey. The energyregions E g , E ∗ and the subbands S , SP , T , SE and HE aremarked by the respective bars. The grey scaled interactioninterval g C in ( a ) denotes the region of avoided crossings be-tween the SE and T as well as the T and HE classes for thecase of a shallow lattice V = 4. The solid boxes indicate theSF to MI transition. For better visibility of the avoided cross-ings in the vicinity of g = 0 (see solid boxes), ( a ) and ( b )provide the lowest 10 eigenenergies of ( a ) and ( b ) respectively. First we focus on the case of a relatively shallow triplewell, namely V = 4 . 0, see Fig. 1 ( a ). For very smallinteractions, g ≃ 0, the MB eigenstates are energeti-cally categorized according to their corresponding par-ticle configuration in terms of single-particle bands. Forinstance, Fig. 1 ( a ) shows that the eigenstates of thesystem are predominantly bunched onto two energy re-gions denoted by E g and E ∗ respectively. The eigen-states lying within E g possess zero higher-band excita-tions, while those bunched onto E ∗ refer to states withone single-particle excitation to the first excited band.The width of the aforementioned energy regions (band-width) depends on the tunneling coupling between thedifferent sites. Note here that the term tunneling cou-pling refers to the corresponding inverse tunneling rate[7, 8]. The distance between E ∗ and E g (band gap) ischaracterized by the band gap between the ground andthe first excited band of the non-interacting system. Re-garding the decomposition of each eigenstate in terms of(localized) Wannier number states [see Eq. (3)] it turnsout that it is an admixture of all the energetic classes S , SP and T . The latter indicates the spatial delocaliza-tion of the bosons within the triple well, and thereforemanifests the few-body analogue of the SF phase for lowinteractions.For increasing interaction strength the energy expec-tation value of the number states belonging to the SP and T classes increases strongly. The same holds for theeigenenergies of the eigenstates to which the aforemen-tioned number states are contributing. For 0 < g / . a ).These avoided crossings manisfest the tunneling couplingbetween the S , SP and T number states of the same par-ity [66, 75]. The aforementioned interaction regime cor-responds in our few-body system to the region where thetransition from the SF to the MI phase takes place. For g ≥ . a ). Theground state of the system is dominated by the S classmanifesting the few-body analogue of the MI phase. The SP and T class dominated eigenstates are also bunchedtogether forming the SP and T subbands. Moreover, theeigenstates of the T subband (being the most sensitive tointerparticle repulsion) experience wide (see solid circles)and narrow (see dashed circles) [76] avoided crossingswith the eigenstates possessing a higher-band excitation.The wide avoided crossings are related to the onset ofthe cradle mode [59, 60] and are a consequence of theinteraction induced decay of an SP or T state causedby the scattering of one of the bosons that reside in thesame well to the first excited state of an adjacent site[77]. This latter process gives rise to the so-called cradlemode which represents a dipole-like intrawell oscillationin the outer wells of the finite lattice (for a detailed de-scription on the generation and properties of this modesee [59, 60]).The same overall behavior can be observed for deeperlattices, see for instance in Fig. 1 ( b ) the case of V = 10.As shown, the differences between shallow and deep lat-tices are mainly quantitative [78]. The bandwidth of eachsubband decreases as a consequence of the reduced tun-neling coupling, while the corresponding subband energygap increases. Most importantly, the transition from theSF to the MI state is realized for smaller interactionswhen compared to the case of V = 4 [see in particularFigs. 1 ( a ) and ( b )]. Moreover, as a consequence of the increased band gap, the positions of the avoided cross-ings related to the cradle mode (see solid circles in Fig.1) occur at larger interparticle repulsions g . Concluding,all differences in the eigenspectrum caused by V can betraced back to the increased intrawell localization of theparticles for deeper lattices.In the following sections, we study the dynamical re-sponse induced by a LIQ from g i = 0 to g f = 2 and back.The choice of g i = 0 ensures that the initial state withina positive LIQ consists of an admixture of all energeticclasses, while due to the postquench interaction, g f = 2,the system remains adequately detuned from the cradlemode even in the fully diabatic (i.e. abrupt) limit. More-over, following a negative LIQ from g i = 2 to g f = 0 thesystem initially resides in a MI state and finally in an ad-mixture of several lowest band states. Taking advantageof the above presented eigenspectra, we expect that thehigher-band dynamics strongly depends on the rampingrate of the SF to MI transition resulting in a differentpopulation of the T class (eigen)states. B. Quench dynamics from SF to MI phase Let us first focus on the positive LIQ dynamics tostrong interactions, namely from g i = 0 to g f = 2,of a system consisting of three initially non-interactingbosons confined in a triple well. The initial MB state isan admixture of the available lowest band number states[see also Fig. 1 ( a )], where the main contribution stemsfrom the SP category due to the hard wall boundaryconditions. To gain an overview of the system’s dynam-ical response induced by the LIQ we employ the fidelity F ( t ; τ ) = |h Ψ τ (0) | Ψ τ ( t ) i| , being the overlap betweenthe time evolved and the initial (ground) state [79–81].Fig. 2 ( a ) presents F ( t ; τ ) for a shallow lattice ( V = 4)with varying ramp time τ . As it can be seen, F ( t ; τ ) per-forms oscillations in time with multiple frequencies, whilefor increasing ramp time the system significantly departsfrom its initial state [i.e. F ( t ; τ ) overall decreases]. Inter-estingly enough, for intermediate τ ’s the fidelity fluctu-ates more prominently, see e.g. 15 < τ < 35, indicatingan enhanced dynamical response. The same overall re-sponse is observed for deeper lattices, see for instanceFig. 2 ( d ), but the corresponding signatures of enhancedresponse for larger τ ’s become more faint due to the in-creased energy gap in the respective MB eigenspectra [seealso Fig. 1]. To further elaborate on the overall responseof the system induced by the LIQ we show in Fig. 2( b ) the time averaged fidelity of Fig. 2 ( a ) for differ-ent ramp times, namely ¯ F ( τ ) = T − τ R Tτ dtF ( t ), where T = 500. We observe that for increasing τ , ¯ F ( τ ) mainlydecreases possessing also some small amplitude deforma-tions displayed as dips in ¯ F ( τ ). This overall decrease of¯ F ( τ ) implies that the system departs more prominentlyfrom its initial state for increasing τ and it is a man-ifestation of the Landau-Zener mechanism [19, 82, 83].Namely, the more adiabatically we drive the system itdeparts stronger from its initial state. Finally, it is ob-served that for large enough τ , namely τ > 70, ¯ F tendsto a constant value which indicates a tendency to ap-proach the adiabatic ramping rate. The small amplitudedeformations in ¯ F ( τ ) refer to different resonant responseregions at specific intervals of τ . These strong responseregions indicate that the specific combination of τ ’s andquench amplitude are more efficient to cross the exist-ing avoided crossings [see Fig. 1] and as a consequenceto depart from the initial state. This latter behavior isalready imprinted in F ( t ; τ ) at the final instant of theramping, i.e. t = τ . Indeed, within the τ intervals that¯ F ( τ ) exhibits dips (humps) the system departs signifi-cantly (negligibly) from its initial state at time t = τ , seealso Fig. 2 ( a ). A careful inspection of the eigenspectra,shown in Fig. 1 ( a ), indicates that for smaller V ’s thegap between the different energetic subbands reduces andtherefore the corresponding tunneling processes can bemore pronounced. The latter suggests that for shallower[deeper] lattices the system is perturbed more efficiently[inefficiently] resulting in an enhanced [reduced] dynam-ical response. To examine in more detail the dynamicalresponse, induced by the LIQ, we employ the normalizedvariance of the fidelity K ( τ ) = p σ F ( τ )¯ F ( τ ) , (5)where σ F ( t ) = T − τ R Tτ dt [ F ( t ; τ ) − ¯ F ( τ )] denotes thetemporally integrated variance of the fidelity. K ( τ )serves as a measure for the mean fidelity variation fromits mean value and it is bounded to take values withinthe interval [0 , K ( τ ) → K ( τ ) → τ ∈ (5 , , (15 , 35) the system can be drivenaway from its initial state in the most prominent manner[ K ( τ ) increases] and therefore the dynamical response ismaximized (see also below).To identify the participating tunneling modes inducedby the positive LIQ, we inspect the fidelity spectrum[61, 62], see Fig. 2 ( e ). As shown the LIQ triggers fivedistinct tunneling modes onto the system. The first fourmodes located at ω , , , ≈ . , , . , . V = 10, g f = 2 refer to intraband tunneling within the SP class, SP → S , SP → T , and T → S categories respectively.The mode at ω ≈ . 75 indicates interband interactionassisted tunneling between the S and HE states. Notehere that the SP → S , and T → S lowest band modesrefer to second order tunneling [63–65], while the othersdenote single-particle transport. The positions of thesefrequency branches remain insensitive to varying τ , seethe inset in Fig. 2 ( e ). It is also worth mentioning atthis point that the dominant dynamical modes SP → S and SP → T are responsible for the maximized dynami-cal response [see Fig. 2 ( c )] within the ramping intervals τ ∈ (5 , 10) and (15 , 35) respectively. To highlight the cor-related MB character of the above mentioned tunneling ω [ ω R ] A F [ a . u ] 20 40 60 80 1001234 τ [ ω − R ] ω [ ω R ] t [ ω − R ] τ [ ω − R ] ¯ F ( τ ) 20 40 60 80 1000.30.40.50.6 τ [ ω − R ] K ( τ ) τ [ ω − R ] ( a )( d ) ( b )( c )( e ) V = 10, g i = 0, g f = 2, MB V = 10, g i = 0, g f = 2, MF V = 4, g i = 0, g f = 2, MB V = 10, g i = 0, g f = 1, MB V = 10, g i = 0, g f = 3, MB FIG. 2. Fidelity evolution for varying ramp time τ after a LIQfrom g i = 0 to g f = 2. Shown are the cases of ( a ) a shallow V = 4 and ( d ) a deep V = 10 triple well. ( b ) Mean fidelity¯ F ( τ ) and ( c ) normalized variance of the fidelity K ( τ ) in thecase of a shallow lattice for varying τ . ( e ) Spectrum of thefidelity F ( ω ), following a LIQ with τ = 8 for different barrierheights, initial and final interactions within the MB approachor the MF approximation (see legend). Inset: F ( ω ) followinga LIQ from g i = 0 to g f = 2 for varying ramp time τ . Thesystem consists of three initially non or weakly interactingbosons (see legend) confined in a triple well. processes we also show the corresponding fidelity spec-trum calculated within the MF approximation. Here,we observe that the system features only the first twotunneling modes, namely transport within the SP and SP → S categories. These mode frequencies are also pos-itively shifted namely ω ′ , ≈ . , . 35 when comparedto the MB approach. As a next step, let us demonstratepossible control protocols of the tunneling dynamics. Anefficient way to manipulate the transport frequency isto tune the height of the potential barriers. Then, thetunneling frequency can be enhanced by using a shallowlattice as the corresponding energy gaps between the dif-ferent subbands become smaller [compare Figs. 1 ( a ),( b )], thereby making each tunneling process more favor-able. Indeed, as shown in Fig. 2 ( e ) the single-particletunneling modes for V = 4 are located at larger frequen-cies ω ′′ , ≈ . , . V = 10 with theformer being the most energetically favorable as it pos-sesses the highest amplitude. On the other hand, the twoparticle transport modes, namely SP → S and S → T ,are negatively shifted with respect to their correspond-ing frequencies for V = 10. The latter is a consequenceof the decreasing subband energy gap occuring for shal-lower lattices, see also Fig. 1. Alternatively, a similar FIG. 3. Mean excitation probability of at least one bosonto reside in an excited band for varying ( a ) ramp time τ ofthe LIQ, ( b ) barrier height V and ( c ) postquench interactionstrength g f . In all cases we consider a LIQ from an initiallynon to a final strongly correlated state unless stated otherwise(see legends). Different curves correspond to different param-eter values (see legend), while the lines belong to a numericalfitting and provide a guide to the eye. The time-scale τ a de-notes the border between diabatic and moderate ramping for V = 10. The grey scaled interaction interval in ( c ) denotedby g C refers to the region of avoided crossings when V = 4between the SE and T as well as the T and HE categories,see also Fig. 1 ( a ). manipulation of the tunneling frequency can be achievedby tuning the postquenched state, by e.g. the value of g f . We observe that for smaller (larger) g f ’s the corre-sponding tunneling branches are negatively (positively)shifted when compared to g f = 2 and V = 10 [see Fig.2 ( e )]. It is important to note here that the frequencypeaks located at ω ′′′ ≈ . , . SP states and the S class. Sum-marizing, it has been shown that by tuning either theheight of the potential barrier or the postquenched statevia g f we can manipulate the location and the intensityof the tunneling frequency branches.Next we investigate the excitation dynamics, and inparticular the time averaged probability of finding atleast one boson within a single-particle state of a higher-band [see also Eq. (4)]. Fig. 3 ( a ) presents ¯ P exc ( τ ) forvarying τ . As it is evident and further confirmed by a di-rect numerical fit, the mean excitation probability obeysthe bi-exponential law¯ P exc ( τ ) = Ae − τ/τ + Be − τ/τ , (6)where A , B are positive constants. Here, the two timescales introduced by the positive constants τ , τ of the bi-exponential are necessary in order to describe accu-rately the mean excitation dynamics covering the suddento adiabatic interaction quench regimes. The existenceof these two time scales can be explained by the behaviorof the number states that contribute to the MB state ateach ramping interval. For convenience let us refer to theramping time at the border between diabatic and moder-ate to nearly adiabatic ramping as τ a [see also Fig. 3 ( a )].Then, following the LIQ for τ < τ a where ¯ P exc ( τ ) exhibitsa rapid decay, the avoided crossing located at g = 0 inthe eigenspectra is crossed diabatically and the MB stateafter the quench consists of a superposition of all lowestband number state classes. However for τ > τ a , where¯ P exc ( τ ) features a slow decay, the crossing at g = 0 is tra-versed in a more adiabatic manner and the system afterthe quench is again in a superposition of the above men-tioned states but now the state | , , i possesses the maincontibution. Alternatively, the above description can beexplained by the fact that the SP and T categories forvarying τ exhibit different scalings. Indeed, for τ << τ a the population of the SP increases while the T slowly de-creases for increasing τ . However, for τ ≥ τ a both cate-gories decay in favor of the S state. Having described themechanism behind this bi-exponential decay of ¯ P exc ( τ )let us demonstrate whether this law is robust also withinthe MF approximation or upon varying the initial stateof the system. As it can be seen in Fig. 3 ( a ) ¯ P exc ( τ )exhibits the same bi-exponential law also within the MFapproximation but the predicted amount of excitationsis lesser when compared to the MB approach. There-fore, we can infer that the MF approximation predictsthe qualitative decay of excitations but fails to capturethe quantitative amount of excitations. Turning again tothe MB approach, we observe that by initializing the sys-tem in a weakly interacting state ¯ P exc ( τ ) shows the samedecay law but fewer excitations are produced than start-ing from g i = 0. This can be explained by the fact thatthe population of the T category in the initial state isstrongly suppressed for increasing g i . Additionally, Fig.3 ( a ) shows ¯ P exc ( τ ) for a shallower lattice. ¯ P exc ( τ ) againexhibits a bi-exponential decay but most importantly itis observed that at τ ≈ 19 the curves for V = 10 and V = 4 cross each other. As a consequence for τ < τ > 19] ¯ P exc ( τ ) is larger [smaller] for the shallower lat-tice. The above-mentioned crossing between the ¯ P exc ( τ )for different heights of the potential barriers serves as astarting-point for investigating in the following the be-havior of the mean excitation probability for varying V and fixed τ , see Fig. 3 ( b ). A compatible double Gaussianfit ¯ P exc ( V ) = A e − ( V − C C ) + B e − ( V − D D ) (where A , B , C , C , D , D refer to positive constants) is pro-vided as a guide to the eye. As it can be seen, the scalingof ¯ P exc ( V ) depends strongly on the considered rampingtime τ . Indeed, proximally to the diabatic limit, e.g. see τ =1, ¯ P exc ( V ) decreases for increasing V . Interestinglyenough, for larger τ , e.g. see τ =25, ¯ P exc ( V ) exhibits acompletely different behavior. Initially it increases forincreasing V until it reaches a maximum value and thendecreases for larger V . For instance, at τ = 25 ¯ P exc ( V )increases until V ≈ V > 6. We also remark here thatthe maximum of ¯ P exc ( V ) is displaced for varying τ , e.g.for τ = 10 is located at V ≈ . τ = 50 isat V ≈ 8. The above described alternating behavior of¯ P exc ( V ) for the different τ regions can be traced back tothe previously observed crossing of ¯ P exc ( τ ) for different V , shown in Fig. 3 ( a ) at τ = 19. Indeed, for τ < τ > 19) ¯ P exc ( τ ) decreases (increases) for increasing V ,i.e. ¯ P exc ( τ < , V = 4) > ¯ P exc ( τ < , V = 10).This is due to the fact that the fast (slow) decay pro-cess described by τ ( τ ) being related to the T → SP ( SP → S ) mode is stronger (weaker) for increasing V ,see also Fig. 3 ( a ). Comparing ¯ P exc ( V ) for a fixed V and varying τ , we observe that ¯ P exc ( V ) decreases in auniform manner for increasing τ , a result that is in ac-cordance with our previous observations for the ¯ P exc ( τ )bi-exponential decay, see Fig. 3 ( a ). However, for deeperlattices ( V > 12) the decrease of ¯ P exc ( V ) almost haltsafter a certain τ (mainly τ > c ). ¯ P exc ( g f ) exhibits two distinct response re-gions with respect to g f . These regions are classified bythe location g C of the existing avoided crossings betweenthe SE and T as well as T and HE categories in the cor-responding MB eigenspectrum, see the grey scaled areawhich incorporate the solid elipses in Fig. 1 ( a ). Indeed,when g f approaches g C ¯ P exc ( g f ) exhibits a non-lineargrowth as shown in the grey scaled interaction intervalof Fig. 3 ( c ), while for g f > g C it increases in a linearmanner. Overall, ¯ P exc ( g f ) grows for larger g f ’s as weimport more energy to the system and thus excite morehigher-band states. On the contrary, ¯ P exc ( g f ) reducesfor larger τ ’s because for more adiabatic LIQs the effectof the existing avoided crossings is smeared out.In the next subsection we examine the negative LIQdynamics, within the triple well system, from the MI tothe SF correlated regimes. In particular, we shall elabo-rate in detail how such a negative LIQ alters the overallresponse of the system consisting of inter and intrabandtunneling dynamics. C. Quench dynamics from MI to SF phase The system is initialized within the strongly correlatedregime, g i = 2, and therefore the ground state corre-sponds to the spatial contribution | , , i . To inducethe dynamics we perform a negative LIQ to a weaklyor non-interacting state and examine the dynamical re-sponse of the system. Fig. 4 ( a ) presents F ( t ; τ ) in ashallow triple well, V = 4, for varying ramp time τ and g f = 0. As in the positive LIQ scenario F ( t ) exhibitsoscillations during the evolution possessing here, how-ever, only a few small value frequencies. The fact that F ( t ; τ ) < S , SP and T classes, see also Fig. 1 ( a ). Tojustify the latter, we resort to the probability of specificnumber states that belong to the above classes. Indeed,Figs. 4 ( c ), ( d ) and ( e ) present P ( t ) = | h , , | Ψ( t ) i | , P ( t ) = | h , , | Ψ( t ) i | = | h , , | Ψ( t ) i | and P ( t ) = | h , , | Ψ( t ) i | respectively. We observe that for dia-batic LIQs the three modes are almost equally populated,while for large τ ’s the SP mainly contributes to the dy-namics and the other states possess a decaying amplitudein time. The reduced amplitude of | , , i for larger τ indicates that the system is significantly perturbed forproximally adiabatic LIQs. To further elaborate on thesystem’s dynamical response, Fig. 4 ( f ) illustrates thecorresponding normalized variance of the fidelity K ( τ )for varying τ . We observe an almost monotonic decreaseof K ( τ ) for increasing τ , especially for τ > 10, suggestingthat the dynamical response is enhanced in the abruptlimit and reduces when the adiabatic limit is approached.The same overall response is also observed for a deepertriple well, see Fig. 4 ( b ). However, the correspondingdecaying response for increasing τ is more faint due tothe increased energy gap between the MB eigenstates,see also Fig. 1.To investigate the tunneling dynamics we rely on thefidelity spectrum, see Fig. 4 ( g ). It is observed that,in the case of a deep lattice ( V = 10), following a LIQto g = 0 only the S → SP tunneling mode is excitedwhich is located at ω ≈ . 1. Employing the correspond-ing MF approximation this tunneling mode hardly sur-vives with frequency ω ′ = 0 . 001 [hardly visible in Fig.4 ( g )]. Turning again to the MB correlated approach,we investigate the influence on the tunneling dynamicsof the postquench state and the dependence on the po-tential barrier. Examining first the situation of a weaklycorrelated postquench state, namely g f = 0 . 05, we ob-serve that two lowest band tunneling modes exist. Thefirst mode located at ω ≈ . 02 refers to the transport S → SP and the second at ω ≈ . 11 to tunneling S → T . On the other hand, following a LIQ in a shallowlattice to a non-interacting final state mainly three tun-neling modes are triggered. The first located at ω ′ = 0 . S → SP , while the remaining twopossess frequencies ω ′ = 0 . 35 and ω ′′ = 0 . 45 which cor-respond to tunneling between the S and different statesbelonging to the T class. Note that higher frequencypeaks, e.g. at ω ≈ . 62, refer to higher order lowestband transitions, for instance SP → T , and possess re-duced amplitudes.To complement our study, let us explore the meanamount of higher-band excitations induced by a negativeLIQ. As in the previous section we examine the impact ofthe mean excitation dynamics as a function of the ramp-ing time τ [see Fig. 5 ( a )] or the height of the potentialbarrier V [see Fig. 5 ( b )]. Here, we observe that themean excitation dynamics for varying ramp time obeysan exponential decay¯ P exc ( τ ) = A e − τ/τ , (7) ω [ ω R ] A F [ a . u ] V = 10, g i = 2, g f = 0, MB V = 10, g i = 2, g f = 0 . 05, MB V = 4, g i = 2, g f = 0, MB V = 10, g i = 2, g f = 0, MF τ [ ω − R ] τ [ ω − R ] τ [ ω − R ] 20 40 60 80 100 K ( τ ) t [ ω − R ] τ [ ω − R ] t [ ω − R ] τ [ ω − R ] t [ ω − R ] τ [ ω − R ] ( a ) ( b )( d )( f )( g )( c ) ( e ) FIG. 4. Fidelity evolution for varying ramp time τ aftera LIQ from g i = 2 to g f = 0. Shown are the cases of( a ) a shallow V = 4 and ( b ) a deep V = 10 triple well.Probabilities of specific number state configurations duringthe evolution, namely ( c ) P ( t ) = | < , , | Ψ( t ) > | , ( d ) P ( t ) = | < , , | Ψ( t ) > | = | < , , | Ψ( t ) > | and ( e ) P ( t ) = | < , , | Ψ( t ) > | . ( f ) Normalized variance ofthe fidelity K ( τ ) in the case of a shallow lattice for varyingramp time τ . ( g ) Spectrum of the fidelity following a LIQwith τ = 8 for different barrier heights, initial and final in-teractions within the MB approach or the MF approximation(see legend). The system consists of three initially stronglyinteracting bosons (see legend) confined in a triple well. where A , τ correspond to positive constants, and τ characterizes the rate of the decay. As shown the amountof produced excitations is negligible for all ramp times,see Fig. 5 ( a ). It reduces for a LIQ to a weakly correlatedinstead of the non-interacting state [see the inset of Fig. 5( a )] and it is enhanced for shallower lattices. Finally, thecorresponding MF approximation shows the same qual-itative behavior but quantitatively fails to predict thecorrect amount of excitations. In addition, we find thatthe mean excitation dynamics follows again an exponen-tial decay with respect to the barrier height, namely¯ P exc ( V ) = A e − V /V , (8)where A , V refer to positive constants and V is theinverse rate of the decay. As it can be deduced from Fig.5 ( b ) again the produced amount of excitations is negli-gible and reduces even further when a larger ramp timeis considered. Finally, let us examine the dependence ofthe excited to higher-band fraction on the postquenchstate for fixed ramping rate, shown in the inset of Fig. 5 τ [ ω − R ] ¯ P e x c ( τ ) V [ E R ] ¯ P e x c ( V ) V = 10, g i = 2, g f = 0, MB V = 10, g i = 2, g f = 0, MF V = 10, g i = 2, g f = 0 . 05, MB V = 4, g i = 2, g f = 0, MB τ = 0 . τ = 1 τ = 10 τ = 50 τ = 100 g f [ E R k − ] ¯ P e x c ( g f ) × -3 V = 4, τ = 10 V = 4, τ = 100 V = 10, τ = 10 × -5 τ [ ω − R ] ¯ P e x c ( τ ) V = 10, g i = 2, g f = 0, MB V = 10, g i = 2, g f = 0 . 05, MB ( b )( a ) FIG. 5. Mean excitation probability of at least one boson toreside in an excited band for varying ( a ) ramp time τ of theLIQ (see a magnified version in the inset), ( b ) barrier height V and in the inset of ( b ) postquench interaction strength g f with g i = 3. In all cases we consider a LIQ from an initiallystrongly to a final non correlated state unless stated otherwise(see legends). Different curves correspond to different param-eter values (see legend), while the lines stem from a numericalfitting and provide a guide to the eye. ( b ). We observe that the population of excited states isoverall negligible, and in particular it is greatly reducedwhen we enter the SF regime namely g f < . IV. LIQ DYNAMICS IN EXTENSIVE UNITFILLING LATTICE SYSTEMS Let us now investigate the dynamical response forlarger unit filling setups. Here, we mainly focus on afive well optical lattice and consider a LIQ from a SFto a MI-like state and vice versa. Concerning the ini-tial state of the system in the SF phase it consists ofan admixture of the Wannier number states | , , , , i , | , , , , i , | , , , , i , | , , , , i , | , , , , i , whilein the MI phase the dominant contribution stems fromthe | , , , , i state. Note also that due to the underly-ing spatial symmetry of the system all the correspondingparity symmetric states contribute as well. FIG. 6. Fidelity evolution for varying ramp time τ after apositive LIQ from g i = 0 to g f = 2 within ( a ) the MB ap-proach and ( b ) the MF approximation. ( c ), ( d ) The same as( a ), ( b ) but following a negative LIQ from g i = 2 to g f = 0.The system consists of five initially ( a ), ( b ) non and ( c ), ( d )strongly interacting bosons confined in five wells. ¯ P exc ( τ ) forvarying ramp time τ of the ( e ) positive and ( f ) negative LIQ.( g ), ( h ) The same as ( e ), ( f ) but for varying barrier height V . Different curves correspond to different parameter values(see legends), while the lines stem from a numerical fittingand provide a guide to the eye. The legend shown in ( f ) [( h )]is the same also for ( e ) [( g )]. To trigger the dynamics we employ either a positiveor a negative LIQ protocol. Fig. 6 ( a ) shows F ( t ; τ ) fol-lowing a positive LIQ for varying τ . As in the triplewell case, we observe that F ( t ; τ ) exhibits an oscilla-tory behavior in time possessing also multiple frequen-cies which mainly correspond to lowest band transportand a few interband tunneling modes. The amplitudeof the F ( t ; τ ) oscillation becomes larger for increasing τ when referring to a fixed time instant t . This indicatesthat the system maximally departs from its initial statein the proximity of an adiabatic LIQ. Then by exploiting τ we can adjust the sweeping rate of the avoided cross-ings in the corresponding MB eigespectrum (not shownhere for brevity) and thus control the population of theresulting excitations. Interestingly enough and in con-trast to the triple well system here ¯ F ( t > τ ; τ ) < . 3, suggesting that larger unit filling setups can be drivenout-of-equilibrium in a more efficient manner which is amanifestation of the Anderson orthogonality catastrophe[81, 84, 85]. As a consequence more modes can be trig-gered pointing also to the opportunity for enhanced in-terband tunneling (see below). To demonstrate the MBcharacter of the dynamics we next present F ( t ; τ ) em-ploying the corresponding MF approximation, see Fig. 6( b ). Despite the overall tendency for enhanced dynamicalresponse for increasing τ , F ( t ; τ ) differs significantly fromthe MB approach. Indeed, a reduced number of modes isparticipating most of which refer to single-particle trans-port while 0 . < ¯ F ( t > τ ; τ ) < . a ). The existing modes possess positive shiftedfrequencies, when compared to the MB approach.Next we examine the dynamical response induced bya negative LIQ, namely the MI melting, employing again F ( t ; τ ), see Fig. 6 ( c ). We observe that F ( t ; τ ) per-forms oscillations of both large amplitude and period,while for larger τ ’s the system becomes more perturbed[hardly visible in Fig. 6 ( c )]. After the LIQ the MB stateconsists of a superposition of several lowest band num-ber states among which the states of double occupancy,e.g. | , , , , i , possess the main contribution and be-come dominant as we tend to an adiabatic ramping. Theabove result is in agreement with the negative LIQ dy-namics of the triple well, however here F ( t ; τ ) performsoscillations of both larger amplitude and period. Withinthe corresponding MF approximation, see Fig. 6 ( d ), F ( t ; τ ) shows the same overall qualitative decrease forlarger τ ’s but quantitatively the dynamical response isaltered significantly. For instance, F ( t ; τ ) and its oscil-lation frequencies are larger when compared to the MBapproach.Turning to the investigation of the induced interbandtunneling we employ the mean excitation probability¯ P exc ( t ; τ ) for varying τ . Following a positive [negative]LIQ ¯ P exc ( τ ) obeys a bi-exponential [an exponential] de-cay law, see Fig. 6 ( e ) [Fig. 6 ( f )] similar to the triplewell case. Note also that, as before, following a nega-tive LIQ ¯ P exc ( τ ) is negligible. Additionally, ¯ P exc ( τ ) isenhanced for larger unit filling setups, while the MF ap-proximation predicts a smaller (enhanced) amount of ex-citations for positive (negative) LIQs. Following a posi-tive LIQ there is a crossing point (here at τ = 8) between¯ P exc ( τ )’s, as in the triple well case, that refer to differ-ent V ’s. This behavior of ¯ P exc ( τ ) seems to be quite ro-bust in different setups and suggests the existence of twoexcitation time scales concerning the ramping rate and V . Namely, within a positive diabatic (adiabatic) LIQthe excited to higher-band fraction is larger (smaller) forshallower lattices. To complement our study on the ex-citation dynamics we investigate ¯ P exc ( V ) for varying V both for a positive LIQ, see Fig. 6 ( g ), as well as a neg-ative LIQ, see Fig. 6 ( h ). We observe that within apositive LIQ scenario ¯ P exc ( V ) strongly depends on τ . Inagreement to the triple well case, ¯ P exc ( V ) decreases forsmall τ [e.g. see τ = 1 in Fig. 6 ( g )], while for larger0rampings ¯ P exc ( V ) shows a maximum at a specific regionof V [e.g. see τ = 10 in Fig. 6 ( g )]. However, for anegative LIQ ¯ P exc ( V ) follows an exponential decay withrespect to V . Overall, in both positive and negative LIQsthe excited to higher-band mean fraction increases for amore diabatic ramping and also for larger systems of unitfilling. Finally, let us note that for the five well system¯ P exc ( V ) is not negligible suggesting that for more ex-tensive unit filling lattices the occupation of higher-bandstates might be unavoidable even in the course of thenegative LIQ dynamics. V. CONCLUSIONS AND OUTLOOK We have explored the nonequilibrium quantum dynam-ics following a linear interaction quench protocol in repul-sively interacting few boson ensembles confined in finiteoptical lattices. The focus has been on unit filling suchthat the ground state of the system for increasing inter-action strength exhibits the transition from a SF to aMI phase. To realize this transition and obtain the in-teraction dependence of the occupation of the numberstates, we first calculate the many-body eigenspectrumfor varying interparticle repulsion. Here, the existence ofmultiple avoided crossings before and after the transitionis elucidated. To induce the dynamics we perform a LIQand cross dynamically, with a finite ramp rate, the afore-mentioned transition from both directions. Subsequently,we explore the dynamical response caused by the LIQ andin particular examine its dependence on several systemparameters, such as the height of the potential barrier.It is important to note here that within our multibandtreatment the dynamical response consists of both thelowest band tunneling and the excited to higher-bandsfraction.Crossing the weak to strong interaction regimes yieldsthe excitation of several lowest band tunneling pathwaysconsisting of single and two particle transport. Further-more, a rich interband tunneling dynamics is observedpossessing mainly a single excitation to the first or sec-ond excited band. Analyzing in more detail the excitedto higher-band fraction we examine its dependence on thequench ramp rate and barrier height. We find that it ex-hibits a bi-exponential decay for decreasing quench rate.This decay law introduces two different time scales in theexcitation dynamics, which are directly related to the di-abatic or adiabatic crossing of the transition respectivelyand can be further explained by the behavior of the par-ticipating number states. Furthermore, the excited frac-tion follows a more complex scaling for varying height ofthe potential barrier. For diabatic quenches it reduces,while for larger ramp times it exhibits a non-linear be-havior. The latter can be explained by exploiting thedependence on the ramp time of the excited fraction forshallow and deep lattices. Finally, the higher-band dy-namics strongly depends on the postquench state, namelywhen we tend to the region of an existing avoided cross- ing it is described by a non-linear growth, while for largerquench amplitudes it increases in an almost linear man-ner.In constrast to the above, the dynamical response fol-lowing a LIQ from strong to weak interactions is reducedand mainly comprises of the lowest band tunneling dy-namics. Indeed, in this case following the quench wecan excite only a few tunneling modes, while the excitedto higher-band fraction is negligible and obeys an expo-nential decay for varying ramp time. Here, the lowestband approximation seems to describe the induced dy-namics accurately. Finally, we made an attempt to gen-eralize our results by considering larger systems, e.g. afive well setup, and showing the robustness of the abovementioned scalings as well as the enhancement of the ex-cited to higher-band fraction.Finally, let us comment on possible future extensions ofour work. An interesting alternative of the present workwould be to investigate the dynamical response inducedby a LIQ in repulsively interacting dipolar bosons [86]upon crossing the corresponding SF to supersolid transi-tion point. Certainly, the study of bosonic or fermionicspinor ensembles confined in an optical lattice is an in-triguing perspective. Here, the inclusion of the spin de-gree of freedom enriches the phase diagram [1, 87, 88]and as a consequence it might alter significantly the dy-namical response. Appendix A: The Computational Approach:MCTDHB To solve the many-body (MB) Schr¨odinger equation( i ~ ∂ t − H ) | Ψ( t ) i = 0 of the interacting bosons as aninitial value problem | Ψ(0) i = | Ψ i , we rely on theMulti-Configuration Time-Dependent Hartree methodfor Bosons (MCTDHB) [73, 74, 89]. The latter hasalready been applied for a wide set of nonequilibriumbosonic settings, e.g. see [59–62, 75, 89–93]. Thismethod allows for a variationally optimal truncation ofthe Hilbert space as we employ a time-dependent movingbasis where the system can be instantaneously optimallyrepresented by time-dependent permanents. The MBwavefunction is expanded in terms of the bosonic numberstates | n , n , ..., n M ; t i , that built upon time-dependentsingle-particle functions (SPFs) | φ i ( t ) i , i = 1 , , ..., M ,and time-dependent weights C ~n ( t ) | Ψ( t ) i = X ~n C ~n ( t ) | n , n , ..., n M ; t i . (A1)Here M is the number of SPFs and the summation ~n isover all the possible combinations n i such that the to-tal number of bosons N is conserved. Note that in thelimit in which M approaches the number of grid pointsthe above expansion is equivalent to a full configura-tion interaction approach. Furthermore, in the case of M = 1 the MB wavefunction is given by a single per-1manent | n = N ; t i and the method reduces to the time-dependent Gross Pitaevskii mean-field approximation.To determine the time-dependent wave function | Ψ( t ) i we calculate the equations of motion for the coefficients C ~n ( t ) and the SPFs | φ i ( t ) i . Following the Dirac-Frenkel[94, 95] variational principle, h δ Ψ | i∂ t − ˆ H | Ψ i = 0, weobtain the well-known MCTDHB equations of motion[73, 74, 89]. These equations consist of a set of M non-linear integrodifferential equations of motion for theSPFs being coupled to the ( N + M − N !( M − linear equations ofmotion for the coefficients C ~n ( t ). Finally, let us remarkthat in terms of our implementation we use an extendedversion of MCTDHB being referred to in the literatureas the Multi-Layer Multi-Configuration Time-DependentHartree method for bosonic and fermionic Mixtures (ML-MCTDHX) [96–98]. This computational package is par-ticularly suitable for treating systems consisting of dif-ferent bosonic, fermionic species, while for the case of asingle bosonic species it reduces to MCTDHB.For the numerical implementation, the SPFs are ex-panded within a so-called primitive basis {| k i} of dimen-sion M p . As a primitive basis for the SPFs we have useda sine discrete variable representation, which intrinsicallyintroduces hard-wall boundaries at both ends of the po-tential. To obtain the n -th MB eigenstate we rely on theso-called improved relaxation scheme, being summarizedas follows. First, we initialize the system with an ansatzset of SPFs {| φ (0) i i} , diagonalize the Hamiltonian withina basis spanned by the SPFs and set the n -th obtainedeigenvector as the C (0) ~n -vector. Then, we propagate theSPFs in imaginary time within a finite time interval dτ , update the SPFs to {| φ (1) i i} and repeat the above stepsuntil the energy of the state converges within the pre-scribed accuracy. In turn, we perform a time-dependentquench on the strength of the interparticle repulsion andstudy the evolution of | Ψ( t ) i in the m -well potential byutilizing the appropriate Hamiltonian within the MCT-DHB equations of motion.To track the numerical error and guarantee the ac-curate performance of the numerical integration for theMCTDHB equations of motion we impose the followingoverlap criteria |h Ψ | Ψ i − | < − and |h ϕ i | ϕ j i − δ ij | < − for the total wavefunction and the SPFs respec-tively. The dimension of the used primitive basis con-sists of 300 spatial grid points in the case of a triple welland 500 spatial grid points for the five well potential.Furthermore, to ensure the convergence of our simula-tions we have used up to 9 (10) optimized single particlefunctions for the triple (five) well, thereby observing asystematic convergence of our results. An auxilliary in-dicator for convergence is provided by the population ofthe lowest occupied natural orbital kept always below0 . ACKNOWLEDGMENTS The authors gratefully acknowledge funding by theDeutsche Forschungsgemeinschaft (DFG) in the frame-work of the SFB 925 ”Light induced dynamics andcontrol of correlated quantum systems” and by theexcellence cluster ”The Hamburg Centre for UltrafastImaging-Structure, Dynamics and Control of Matter atthe Atomic Scale”. [1] M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A.Sen, and U. Sen, Adv. in Phys. , 243 (2007).[2] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008).[3] V.I. Yukalov, Laser Phys. , 110 (2009).[4] M. Lewenstein, A. Sanpera, and V. Ahufinger, UltracoldAtoms in Optical Lattices: Simulating Quantum Many-body Systems (Oxford University Press, Oxford, 2012).[5] M. Greiner, O. Mandel, T. Esslinger, T. W. H¨ansch, andI. Bloch, Nature (London), , 39 (2002).[6] M. Greiner, O. Mandel, T. W. H¨ansch, and I. Bloch,Nature (London), , 51 (2002).[7] M. P. Fisher, P. B. Weichman, G. Grinstein, and D. S.Fisher, Phys. Rev. B , 546 (1989).[8] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, andP. Zoller, Phys. Rev. Lett. , 3108 (1998).[9] W.H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. , 105701 (2005).[10] A. Polkovnikov, Phys. Rev. B , 161201 (2005).[11] D. Chowdhury, U. Divakaran, and A. Dutta, Phys. Rev.E , 012101 (2010).[12] K. Sengupta, D. Sen, and S. Mondal, Phys. Rev. Lett. , 077204 (2008).[13] V. Mukherjee, A. Dutta, and D. Sen, Phys. Rev. B ,214427 (2008). [14] A. Dutta, R.R.P. Singh, and U. Divakaran, EPL ,67001 (2010).[15] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalat-tore, Rev. Mod. Phys. , 863 (2011).[16] N. Andrei, arXiv: (2016).[17] T. W. B. Kibble, Physics Today, , 47 (2007).[18] A. Del Campo, and W. H. Zurek, Int. J. Mod. Phys. A, , 1430018 (2014).[19] J. Dziarmaga, Adv. Phys. , 1063 (2010).[20] J. Beugnon, and N. Navon, J. Phys. B: At. Mol. Opt.Phys. , 022002 (2017).[21] T. W. Kibble, J. Phys. A , 1387 (1976).[22] W. H. Zurek, Nature (London) , 505 (1985).[23] G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo, andG. Ferrari, Nat. Phys. , 656 (2013).[24] M. Anquez, B.A. Robbins, H.M. Bharath, M. Bogus-lawski, T.M. Hoang, and M.S. Chapman, Phys. Rev.Lett. , 155301 (2016).[25] L.W. Clark, L. Feng, and C. Chin, Science, , 606(2016).[26] D. Chen, M. White, C. Borries, and B. DeMarco, Phys.Rev. Lett. , 235304 (2011).[27] C. Meldgin, U. Ray, P. Russ, D. Chen, D.M. Ceperley,and B. DeMarco, Nature Phys. , 646 (2016). [28] M. Aidelsburger, J.L. Ville, R. Saint-Jalm, S.Nascimb´ene, J. Dalibard, and J. Beugnon, Phys. Rev.Lett. , 190403 (2017).[29] R. Sch¨utzhold, M. Uhlmann, Y. Xu, and U. R. Fischer,Phys. Rev. Lett. , 200601 (2006).[30] B. Gardas, J. Dziarmaga, and W. H. Zurek, Phys. Rev.B, , 104306 (2017).[31] J. Dziarmaga, and W. H. Zurek, Sci. Rep. 4, 5950 (2014).[32] S. Braun, M. Friesdorf, S.S. Hodgman, M. Schreiber, J.P.Ronzheimer, A. Riera, M. del Rey, I. Bloch, J. Eisert, andU. Schneider, Proc. Nat. Acad. Sci. , 3641 (2015).[33] S. R. Clark, and D. Jaksch, Phys. Rev. A , 043612(2004).[34] J.S. Bernier, D. Poletti, P. Barmettler, G. Roux, and C.Kollath, Phys. Rev. A , 033641 (2012).[35] C.L. Hung, X. Zhang, N. Gemelke, and C. Chin, Phys.Rev. Lett. , 160403 (2010).[36] S.D. Huber, E. Altman, H.P. B¨uchler, and G. Blatter,Phys. Rev. B , 085106 (2007).[37] L. Vidmar, S. Langer, I.P. McCulloch, U. Schneider, U.Schollw¨ock, and F. Heidrich-Meisner, Phys. Rev. B ,235117 (2013).[38] N. Gemelke, X. Zhang, C.L. Hung, and C. Chin, Nature , 995 (2009).[39] M. Rigol, R.T. Scalettar, P. Sengupta, and G.G. Ba-trouni, Phys. Rev. B , 121103 (2006).[40] P. Sengupta, M. Rigol, G.G. Batrouni, P.J.H. Denteneer,and R.T. Scalettar, Phys. Rev. Lett. , 220402 (2005).[41] Y. Yanay, and E.J. Mueller, Phys. Rev. A , 013622(2016).[42] C. Kollath, A.M. L¨auchli, and E. Altman, Phys. Rev.Lett. , 180601 (2007).[43] F. M. Cucchietti, B. Damski, J. Dziarmaga, and W. H.Zurek, Phys. Rev. A , 023603 (2007).[44] P. Navez, and R. Sch¨utzhold, Phys. Rev. A, , 063603(2010).[45] S.S. Natu, K.R. Hazzard, and E.J. Mueller, Phys. Rev.Lett. , 125301 (2011).[46] M. Lacki, D. Delande, and J. Zakrzewski, New J. Phys. , 013062 (2013).[47] J. Kinnunen, and P. T¨orm¨a, Phys. Rev. Lett. , 070402(2006).[48] N. Proukakis, S. Gardiner, M. Davis, and M. Szy-maska, Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics, Vol. , World Scientific (2013).[49] J. Grond, A. I. Streltsov, L. S. Cederbaum, and O. E.Alon, Phys. Rev. A , 063607 (2012).[50] J. Grond, A. I. Streltsov, A. U. Lode, K. Sakmann, L. S.Cederbaum, and O. E. Alon, Phys. Rev. A , 023606(2013).[51] M. Theisen, and A.I. Streltsov, Phys. Rev. A , 053622(2016).[52] R. Beinke, S. Klaiman, L.S. Cederbaum, A.I. Streltsov,and O.E. Alon, Phys. Rev. A , 063602 (2017).[53] D. Blume, Physics , 74 (2010).[54] D. Blume, Rep. Progr. Phys. , 046401 (2012).[55] G. Z¨urn, F. Serwane, T. Lompe, A.N. Wenz, M.G. Ries,J.E. Bohn, and S. Jochim, Phys. Rev. Lett. , 075303(2012).[56] F. Serwane, G. Z¨urn, T. Lompe, T.B. Ottenstein, A.N.Wenz, and S. Jochim, Science , 336 (2011).[57] A. N. Wenz, G. Z¨urn, S. Murmann, I. Brouzos, T. Lompe,and S. Jochim, Science, , 457 (2013). [58] A. M. Kaufman, B. J. Lester, M. Foss-Feig, M. L. Wall,A. M. Rey, and C. A. Regal, Nature, , 208 (2015).[59] S. I. Mistakidis, L. Cao, and P. Schmelcher, J. Phys. B:At. Mol. Opt. Phys. , 225303 (2014).[60] S. I. Mistakidis, L. Cao, and P. Schmelcher, Phys. Rev.A , 033611 (2015).[61] J. Neuhaus-Steinmetz, S. I. Mistakidis, and P.Schmelcher, Phys. Rev. A , 053610 (2017).[62] S. I. Mistakidis, and P. Schmelcher, Phys. Rev. A ,013625 (2017).[63] S. F¨olling, S. Trotzky, P. Cheinet, M. Feld, R. Saers,A. Widera, T. M¨uller, and I. Bloch, Nature , 1029(2007).[64] F. Meinert, M.J. Mark, E. Kirilov, K. Lauber, P. Wein-mann, M. Gr¨obner, A.J. Daley, and H.C. N¨agerl, Science , 1259 (2014).[65] Y.A. Chen, S. Nascimbne, M. Aidelsburger, M. Atala,S. Trotzky, and I. Bloch, Phys. Rev. Lett. , 210405(2011).[66] L. Cao, I. Brouzos, S. Z¨ollner, and P. Schmelcher, NewJ. Phys. , 033032 (2011).[67] L. Cao, I. Brouzos, B. Chatterjee, and P. Schmelcher,New J. Phys. , 093011 (2012).[68] M. Olshanii, Phys. Rev. Lett. , 938 (1998).[69] T. K¨ohler, K. Goral, and P. S. Julienne, Rev. Mod. Phys. , 1311 (2006).[70] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev.Mod. Phys. , 1225 (2010).[71] J. I. Kim, V. S. Melezhik, and P. Schmelcher, Phys. Rev.Lett. , 193203 (2006).[72] P. Giannakeas, F. K. Diakonos, and P. Schmelcher, Phys.Rev. A , 042703 (2012).[73] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, J.Chem. Phys. , 154103 (2007).[74] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys.Rev. A , 033613 (2008).[75] G. M. Koutentakis, S. I. Mistakidis, and P. Schmelcher,Phys. Rev. A , 013617 (2017).[76] The selection rule leading to the formation of exact cross-ings is the intrawell parity symmetry. The latter pro-hibits transitions from the ground state of a particularsite (parity even) to the first excited state of the samesite (parity odd). This selection rule is weakly violatedfor the edge sites due to the hard wall boundary condi-tions leading to very narrow avoided crossings possessinga width δE / . E R k − [75].[77] The described process is a result of the couplings U , , , s,s,s,s ± = R dx | w (0) s ( x ) | w ∗ (0) s ( x ) w (1) s ± ( x ) = 0 which donot fall under the restriction described in [76].[78] The only striking qualitative difference is that the wideavoided crossings between the T and SE eigenstates ex-hibited for shallow lattices [see Fig. 1 ( a )] become narrowin the case of deep lattices [see Fig. 1 ( b )]. The latter re-sults from the fact that a transition from a T to a SE number state requires a tunneling mediated virtual tran-sition to an SP state, which is suppressed for deep lat-tices.[79] L.C. Venuti, and P. Zanardi, Phys. Rev. A , 022113(2010).[80] T. Gorin, T. Prosen, T. H. Seligman, and M. ˇZnidariˇc,Phys. Rep. , 33 (2006).[81] S. Campbell, M.. Garc´ıa-March, T. Fogarty, and T.Busch, Phys. Rev. A , 013617 (2014). [82] L. D. Landau, Phys. Sov. Un. , 2 (1932).[83] C. Zener, Proc. Royal Society of London A , 6 (1932).[84] P.W. Anderson, Phys. Rev. Lett. , 1049 (1967).[85] M. Knap, A. Shashi, Y. Nishida, A. Imambekov, D.A.Abanin, and E. Demler, Phys. Rev. X , 041020 (2012).[86] B. Chatterjee, I. Brouzos, L. Cao, and P. Schmelcher, J.Phys. B: At. Mol. Opt. Phys. , 085304 (2013).[87] S. Tsuchiya, S. Kurihara, and T. Kimura, Phys. Rev. A, , 043628 (2004).[88] M. Rizzi, D. Rossini, G. De Chiara, S. Montangero, andR. Fazio, Phys. Rev. Lett. , 240404 (2005).[89] A. I. Streltsov, O. E. Alon, and L. S. Cederbaum, Phys.Rev. Lett. , 030402 (2007).[90] A. I. Streltsov, K. Sakmann, O. E. Alon, and L. S. Ceder-baum, Phys. Rev. A , 043604 (2011). [91] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys.Rev. A , 013611 (2007).[92] O. E. Alon, A. I. Streltsov, and L. S. Cederbaum, Phys.Rev. A , 022503 (2009).[93] G.C. Katsimiga, S.I. Mistakidis, G.M. Koutentakis, P.G.Kevrekidis, and P. Schmelcher, New J. Phys. (in press)(2017).[94] J. Frenkel, in Wave Mechanics 1st ed. (Clarendon Press,Oxford, 1934), pp. 423-428.[95] P. A. Dirac, (1930, July). Proc. Camb. Phil. Soc. , 376.Cambridge University Press.[96] L. Cao, S. Kr¨onke, O. Vendrell, and P. Schmelcher, J.Chem. Phys. , 134103 (2013).[97] S. Kr¨onke, L. Cao, O. Vendrell, and P. Schmelcher, NewJ. Phys. , 063018 (2013).[98] L. Cao, V. Bolsinger, S.I. Mistakidis, G.M. Koutentakis,S. Kr¨onke, J.M. Schurer, and P. Schmelcher, J. Chem.Phys.147