You are an expert at summarizing long articles. Proceed to summarize the following text: 
out - of - equilibrium dynamics following a quantum quench is a topic of intense study at present . the theme is pursued primarily along two lines . the first one is about the equilibration and thermalization  mechanism of a quantum system @xcite , a fundamental yet still open issue in statistical physics . the second one is about the the real - time dynamical behavior of a many - body system @xcite , which is highly non - trivial in the regime where the quasi - particle picture breaks down .    among all the models investigated so far , the bose - hubbard model takes a special position . as a paradigmatic strongly - correlated model , it can be realized accurately with cold atoms in optical lattices , and especially , the parameters can be controlled ( e.g. changed suddenly ) to a high degree @xcite . this nice property makes it an ideal candidate to investigate quantum quench dynamics both theoretically and experimentally . up to now , in the few theoretical works on the quench dynamics of the bose - hubbard model @xcite , the state of the system before the quench is always assumed to be the ground state of the initial hamiltonian . that is , the system is assumed to be at zero temperature initially . however , in this paper we shall start from a thermal equilibrium state . one should note that this scenario is actually more experimentally relevant . because in current experiments , one generally gets not a single tube of cold atoms , but instead a two - dimensional array of one - dimensional lattices for the cold atoms @xcite . in other words , an _ ensemble _ of one - dimensional bose - hubbard models is obtained in one shot . moreover , in view of the fact that the cold atoms are at finite temperatures necessarily @xcite , it is reasonable to start from a thermal state described by a canonical ensemble density matrix [ see eq . ( [ rhoi ] ) below ] .    as emphasized by linden _ et al_.  @xcite , in the pursuit of thermalization , it is important to distinguish the two closely related but inequivalent concepts of equilibration and thermalization . the latter is much stronger and has the trademark feature of the boltzmann distribution , whereas the former refers only to the stationary property of the density matrix of a ( sub)system or some physical observables . it is highly possible that a system equilibrates but without thermalization . this is actually the case for the bose - hubbard model . as revealed both in previous works ( zero temperature case ) @xcite and in the present paper ( finite temperature case ) , the bose - hubbard model thermalizes only if the quench amplitude is not so large , at least at the finite sizes currently accessible . however , it will be shown below that in a much wider range of parameters , some generic physical observables equilibrate very well . among them are the populations on the bloch states , which are ready to measure by the typical time - of - flight experiment @xcite . remarkably , this is actually guaranteed for these quantities in the thermodynamic limit , i.e. , when the size of the system gets large enough . the equilibration behavior of the physical observables imposes a question . it is ready to recognize that the equilibration of the physical observables is largely an effect of interference cancelation . it never means that the density matrix has suffered any dephasing or decoherence . actually , the density matrix evolves unitarily and in the diagonal representation of the hamiltonian , its elements simply rotate at constant angular velocities . a natural question is then , does the time - dependence of the density matrix has any chance to exhibit it , given that it is almost absent in the average values of the physical observables ? this leads us to consider giving the system a second quench . the concern is , would the system yield different long - time behaviors if the second quench comes at different times ? it turns out that the answer depends on whether the intermediate hamiltonian is integrable or non - integrable . in the former case , the density matrix shows repeated appreciable recurrences and thus the dependence on the second quench time is apparent . in the latter case , on the contrary , the density matrix shows no sign of recurrence and quantitatively similar long - time dynamics is observed for quenches at different times . this paper is organized as follows . in sec . [ sec2 ] , the setting of the problem and the basic approaches are given . in sec . [ sec3 ] , the dynamics after a single quench is studied . the time - averaged density matrix and the real - time evolution of some physical observables are investigated in detail . based on the observation in this section , we proceed to study the scenario of a second quench in sec . [ sec4 ] . finally , we summarize the results in sec . the time - dependent hamiltonian of the bose - hubbard model is ( @xmath1 throughout this paper ) @xmath2 here @xmath3 is the number of sites ( the total atom number will be denoted as @xmath4 ) and @xmath5 ( @xmath6 ) is the creation ( annihilation ) operator for an atom at site @xmath7 . note that here periodic boundary condition is assumed . the parameters @xmath8 and @xmath0 are the nearest - neighbor hopping strength and the on - site atom - atom interaction strength , respectively . note that the dynamics of the system depends only on the ratio @xmath9 , thus we will set @xmath10 throughout . we say the system is quenched if @xmath0 is changed suddenly at some time from one value to another value . experimentally , for cold atoms in an optical lattice , this can be realized by using the feshbach resonance . assume that initially the parameter @xmath0 is of value @xmath11 ( the corresponding hamiltonian is denoted as @xmath12 ) , and the system is at thermal equilibrium and of inverse temperature @xmath13 . denote the @xmath14-th eigenvalue and eigenstate of @xmath12 as @xmath15 and @xmath16 , respectively . the initial density matrix of the system is then @xmath17 is the partition function and @xmath18 is the probability of occupying the eigenstate @xmath19 . note that @xmath20 is the dimension of the hilbert space @xmath21 . the density matrix at time @xmath22 is given formally as @xmath23 , with @xmath24 $ ] . here @xmath25 means time ordering . the hamiltonian @xmath26 is invariant under the translation @xmath27 . this indicates that the total quasi - momentum of the system @xmath28 , where @xmath29 is the creation operator for an atom in the @xmath30-th bloch state , is conserved . this property implies that if the full hilbert space is decomposed into @xmath3 subspaces according to the values of @xmath31 , i.e. , @xmath32 , the hamiltonian and the density matrix are always block - diagonal with respect to the @xmath31-subspaces , i.e.,@xmath33 and @xmath34 @xcite . it is then possible to study the dynamics in each subspace individually ( which saves a lot of computational resource ) and then gather the information together ( note that for the expectation values of quantities like @xmath35 , there are contributions from each subspace ) . here it is necessary to mention that though we should have done the gathering or averaging process for many quantities studied below , we would rather not do so , because it is observed that the system behaves quantitatively similar in all the @xmath31-subspaces @xcite . a single @xmath31-subspace captures the overall behavior very well . therefore , our strategy is to focus on some specific @xmath31-subspace ( @xmath36 actually ) and take the normalization @xmath37 . _ it is understood that in the following all hamiltonians , density matrices , eigenvalues , and eigenstates refer to those belonging to this specific @xmath31-subspace . we will drop the superscript @xmath31 for notational simplicity . _ suppose at time @xmath38 the system is quenched by changing the value of @xmath0 from @xmath11 to @xmath39 , which is then held on forever . the hamiltonian later will be denoted as @xmath40 , and the eigenvalues and eigenstates associated will be denoted as @xmath41 and @xmath42 , respectivley . in the representation of @xmath43 , the density matrix at time @xmath22 is then simply ( in this paper @xmath44 means quantum state averaging while @xmath45 means time averaging ) @xmath46 where @xmath47 is the dimension of the specific @xmath31-subspace . it will prove useful to define the time - averaged density matrix @xmath48 the time - averaged density matrix is of great relevance for our purposes . first , it is both time - independent and variable - independent . second , the time - averaged value of an arbitrary operator @xmath49 is given simply by @xmath50 . that is , the time - averaged density matrix contains the overall information of the dynamics of the system . actually , as we will see later , for some quantities which fluctuate little in time , the time - averaged density matrix tells almost a complete story . third , the process of averaging over time is a process of relaxation in the sense that the entropy associated with @xmath51 is definitely no less than that with the density matrix at an arbitrary time , i.e. , @xmath52 . this is a corollary of the klein inequality @xcite and is reasonable since @xmath53 contains all the information of @xmath51 while the inverse is invalid . the equality also means that @xmath54 will never be damped , and time - averaging is essential .    note that when @xmath55 , generally there is no degeneracy between the eigenvalues of @xmath40 . therefore the time - averaged density matrix is simply diagonal in the basis of @xmath56 , i.e. , @xmath57 with @xmath58 being the population on the eigenstate @xmath59 . in the special case of @xmath60 , the hamiltonian reduces to @xmath61 , with @xmath62 . in this case , each eigenvalue is of the form @xmath63 , under the constraints @xmath64 and @xmath65 , and there can be level degeneracy . however , we can always make some unitary transforms in each degenerate subspace to make sure that @xmath51 is in the form of ( [ diag ] ) . since the time - averaged density matrix provides an overall information of the dynamics of the system , we look into it first . in fig . [ fig1 ] , we consider the scenario of starting from the same initial condition ( @xmath66 , @xmath67 ) but quenching to six different values of @xmath39 @xcite . in each panel , the logarithms of @xmath68 are plotted against the eigenvalues @xmath69 ( red dots ) . we have compared @xmath51 with a canonical ensemble density matrix @xmath70 , which is defined as @xmath71 under the condition @xmath72 . here @xmath73 , the final inverse temperature , is the only fitting parameter . in fig . [ fig1 ] , the green dots which form a straight line correspond to @xmath70 . we see that @xmath51 exhibits many interesting features . in the case of @xmath60 , @xmath51 agrees well with @xmath70 throughout the spectrum . in the case of @xmath74 , @xmath51 agrees well with @xmath70 in the lower part of the spectrum , while deviates from it significantly in the higher part of the spectrum . but overall the two are in good agreement since the weight of the higher part is small . the case of @xmath75 is somewhat the reverse of the @xmath74 case . it is in the lower part of the spectrum that @xmath76 fluctuates wildly . in the higher part @xmath76 goes almost linearly . since the weight is dominated by the lower part , @xmath70 is not a good approximation of @xmath51 . in the strong interaction limits of @xmath77 , another feature takes the place . as a whole the red dots do not fall close to a single straight line , but they do form some stripes , and the stripes are almost parallel with a common slope close to @xmath78 . it is easy to recognize that each stripe corresponds to a bump in the density of states of @xmath40 .    in order to understand the various features in fig . [ fig1 ] , we rewrite @xmath68 as @xmath79 where @xmath80 is a probability distribution @xcite associated with @xmath81 . note that @xmath82 is an intrinsic property of @xmath83 independent of @xmath78 . we have tried to characterize the distribution @xmath82 by its mean @xmath84 , its second central moment @xmath85 , and its third central moment @xmath86 , which are defined as follows ,    [ mean2 ] @xmath87    these quantities are presented in fig . these data enable us to understand fig . suppose for a distribution @xmath82 with @xmath88 , we define a gaussian distribution @xmath89 which shares the same mean and variance with @xmath90 but has vanishing third central moment . replacing @xmath90 in eq . ( [ pm2 ] ) by @xmath91 , we get an approximation of @xmath68 , @xmath92 in fig . [ fig1 ] , @xmath93 are represented by the blue dots . we see that as a whole @xmath93 is a good approximation of @xmath68 , except at the lower part of the spectrum in fig . the reason is clear the @xmath94 s there are the largest throughout all the figures , which indicates that the corresponding distributions @xmath90 are wide and asymmetric and thus can not be well approximated with a gaussian distribution . now we can understand the good fittings in figs . [ fig1]c and [ fig1]d . in these two cases , @xmath84 is almost a linear function of @xmath69 , and @xmath85 does not vary so much , therefore the exponent in eq . ( [ pmprime ] ) goes almost linearly with @xmath69 . the situation is similar in the higher part of the spectrum in fig . [ fig2]b , and therefore we have a good linear fitting for the higher spectrum part in fig .  [ fig1]b . in contrast , in fig . [ fig2]e , @xmath84 varies wildly for adjacent @xmath69 , therefore we see in fig . [ fig1]e large fluctuations about the straight line . as for the parallel stripes in figs . [ fig1]a and [ fig1]f , they are also understandable in terms of figs . [ fig2]a and [ fig2]f , where @xmath84 form parallel stripes . it is numerically checked and can be argued that the slopes of the stripes are almost unity . actually we have @xmath95 where @xmath96 . note that in the limit of large @xmath97 , the kinetic term in the hamiltonian ( [ h ] ) can be viewed as a perturbation to the second interaction term . the spectrum of the latter is highly degenerate and consists of integral multipliers of @xmath39 . the effect of the perturbation is to mix up the eigenstates of the interaction hamiltonian with different eigenvalues and smooth the spectrum . that is why there are bumps in the density of states in figs . [ fig1]a and [ fig1]f and two adjacent bumps are placed roughly @xmath39 apart . by perturbation theory , it is easy to show that the second term in eq . ( [ slope ] ) varies on the order of @xmath98 among eigenstates belonging to the same bump . therefore , approximately we have @xmath99 for each bump and this explains why the stripes in figs . [ fig2]a and [ fig2]f are of slope unity . in turn it explains [ with the help of eq . ( [ pmprime ] ) ] why we have the parallel stripes in figs . [ fig1]a and [ fig1]f , and especially the slopes are approximately @xmath78 . it seems in fig . [ fig1 ] that @xmath70 is a good approximation of @xmath51 only when @xmath100 is small . in fig . [ fig3 ] , we employ the tools of distance @xmath101 , fidelity @xmath102 , and relative entropy @xmath103 ( for the definitions see @xcite ) between two density matrices to quantify the difference or resemblance between @xmath70 and @xmath51 . there it is clear that only in the range of @xmath104 , we have @xmath105 , which means @xmath51 is close to @xmath70 . in the subsequent subsection we will see that only in this range the expectation values of some generic physical observables according to @xmath51 and @xmath70 agree well .     the distance @xmath101 and fidelity @xmath102 between @xmath70 and @xmath51 , and the relative entropy of @xmath70 with respect to @xmath51 , as functions of @xmath39 . the initial state is the same as in fig . [ fig1].,scaledwidth=35.0% ]      we now proceed to study the time evolution of the system after the quench . in fig . [ fig2 ] , we show the time evolution of the populations on the bloch states @xmath106 . the six sub - figures correspond to those in fig . [ fig1 ] respectively . for all the @xmath39 s and all the @xmath30 s , @xmath107 equilibrate to their average values after a transient time , which is relatively longer in the cases of @xmath75 and @xmath108 . in the special case of @xmath60 , there is no fluctuation at all . the reason is simply that in this case , @xmath109 are conserved . we see that the time - averaged values of @xmath110 predicted by @xmath51 ( @xmath111 ) and @xmath70 ( @xmath112 ) agree relatively well in the cases of @xmath60 and @xmath113 . this is consistent with the closeness between @xmath51 and @xmath70 for these two values of @xmath39 , as revealed in fig . [ fig1 ] and fig . [ fig3 ] . here we would say the system thermalizes well in the @xmath74 case , however , we would refrain making the same statement for the @xmath60 case . the reason will be clear in the next section . figure [ fig4 ] is about a finite - sized system with some specific initial condition . however , here we have some general statements . we argue that in the thermodynamic limit ( @xmath114 with @xmath115 fixed ) , as long as initially the system is at finite - temperature thermal equilibrium and described by a canonical ensemble density matrix as ( [ rhoi ] ) , we should see steady behaviors of the physical variables like @xmath109 . let @xmath116 and let @xmath117 . the ensemble - averaged value of @xmath118 at time @xmath22 is @xmath119,\ ] ] where @xmath120 . its time - averaged value is @xmath121 here note that for a generic hamiltonian @xmath40 , there is no level degeneracy . the time - averaged value of @xmath122 is @xcite @xmath123 note that here it is assumed that there is no degeneracy of energy gaps . thus we have for the variance of @xmath124 in time , @xmath125 , @xmath126 since @xmath118 is semi - positive definite and bounded , we have @xmath127 . thus we have @xmath128 here we note that the summation is the square of the frobenius norm of @xmath53 in the representation of @xmath129 , which is invariant in all representations and is preserved by an arbitrary unitary evolution @xcite . explicitly , we have @xmath130    we argue that this quantity , which depends only on the initial state , decays exponentially with the size @xmath3 . let @xmath131 increase with @xmath14 . we have @xmath132 as @xmath133 . here in the @xmath134 relation we used the fact the ground state energy @xmath135 of @xmath12 scales linearly with @xmath3 and so does the free energy @xmath136 of the initial state @xcite . the coefficients @xmath137 and @xmath138 are independent of @xmath3 . moreover , it is easy to see that @xmath139 for any @xmath78 , with the equality taken only in the limit of @xmath140 or @xmath141 , and @xmath142 increases monotonically with @xmath143 . this makes sure that @xmath144 would not grow exponentially with @xmath3 and transcend unity .    with ( [ bound1 ] ) and ( [ bound2 ] ) , we get an upper bound for @xmath145 , @xmath146 where @xmath147 is some constant . the upper bound of @xmath145 helps us determine an upper bound for the probability of finding @xmath124 deviating away from the mean @xmath148 by a distance larger than @xmath149 . actually , following reimann @xcite , using the chebyshev inequality @xcite , we have @xmath150 for a fixed value of @xmath149 , the upper bound decreases exponentially with the size of the system according to ( [ bound3 ] ) . it then follows the statement above . here some comments are worthy . though in the derivation above we have in mind a sudden quench , it is easy to see that the conclusion actually applies to any type of quench ( e.g. , the hamiltonian can be changed continuously over some period , as in @xcite , or quenched multiple times as in sec . [ sec4 ] below ) , as long as after some point the hamiltonian is never changed again . the reason lies in that the frobenius norm of the density matrix @xmath54 is conserved under unitary evolutions , and thus is independent of the historical or the final values of @xmath26 , but is determined entirely by the initial state . as for the operator @xmath118 , only the properties of semi - positive - definiteness and boundedness are used . thus similar conclusions can apply to other operators such as @xmath151 and @xmath152 , or operators in other models . finally , it should be mentioned that the conclusion relies on the fact that the quantity in eq . ( [ fnorm ] ) is bounded by some exponentially decreasing function , which is the case only at finite temperatures ( @xmath153 ) . at zero temperature , the quantity in eq . ( [ fnorm ] ) is always equal to unity and thus the problem is still open . it is shown in fig . [ fig4 ] that after a finite transient time , the physical variables equilibrate to their average values exhibiting minimal fluctuations . moreover , it has been proven that the amplitudes of the fluctuations will decrease exponentially with the size of the system . therefore , the observation is that the system , described by the density matrix @xmath54 , is almost indistinguishable from a system described by the time - averaged density matrix @xmath51 , as far as the simple realistic physical variables are concerned . this is remarkable . because though @xmath54 evolves unitarily and suffers no loss of information of @xmath53 , it behaves as if it were fully decoherenced . the question is then , to what extent can we hold onto this proposition ? is it possible to distinguish @xmath54 and @xmath51 , or @xmath154 and @xmath155 ( @xmath156 ) , by some means ? motivated by this problem , we have considered the scenario of giving the quenched system a second quench . that is , after the first quench at @xmath38 which changes @xmath0 from @xmath11 to @xmath39 , at time @xmath157 , the system is quenched again by changing the value of @xmath0 from @xmath39 to @xmath158 , which is then held on forever . the concern is , would the long - time dynamics of the system depends on the specific time @xmath159 ?    denote the hamiltonian associated with @xmath158 as @xmath160 . the density matrix of the system later is given by @xmath161 ( for @xmath162 ) . as before , we are interested in the long - time averaged value of @xmath54 , @xmath163 since it has been shown and proven above that the dynamics of the system is to a large extent captured by the time - averaged density matrix . here the subscript indicates the dependence on the time @xmath159 . it is also useful to define the average of @xmath164 with respect to @xmath159 , @xmath165 the second equality means that @xmath166 is actually the time - averaged density matrix associated with an initial state @xmath51 [ see eqs . ( [ def_bar_rho ] ) and ( [ diag ] ) ] and a hamiltonian @xmath160 . one purpose of defining @xmath166 is to set a reference state independent of @xmath159 .    to gain an overall idea of the dependence on @xmath159 of the long - time dynamics , we have studied the distance between @xmath164 and @xmath166 @xcite , and the time - averaged value of @xmath167 , @xmath168 as functions of @xmath159 . note that the average value of @xmath169 with respect to @xmath159 is given by @xmath166 , @xmath170 this is another reason for defining @xmath166 . the quantities @xmath171 and @xmath172 are shown in fig . eight pairs of @xmath173 are examined with the same initial condition as in fig . we see that for all cases with @xmath55 , both @xmath101 and @xmath174 set down to their average values quickly . however , for the special case of @xmath60 , both @xmath101 and @xmath175 display repeated recurrences , without any sign of equilibration . the situation is the reverse of that in fig . [ fig4 ] , where @xmath176 does not show any fluctuations in the case of @xmath60 . this phenomenon is due to the recurrence of the density matrix @xmath54 to @xmath53 @xcite . from eq . ( [ rhot ] ) , we see that in the representation of @xmath129 , the @xmath177-th off - diagonal element of @xmath54 rotates at an angular frequency of @xmath178 . in the generic case of @xmath179 , the energy gaps @xmath178 are quite random and incommensurate , and thus recurrence of the density matrix is rare . more precisely , the span between two times when all the matrix elements of @xmath54 get ( nearly ) in phase again is extraordinarily large . on the contrary , in the special case of @xmath60 , all eigenvalues and hence all the energy gaps @xmath178 are integral combinations of the few basic frequencies @xmath180 , and thus the probability of recurrence is much higher to demonstrate that the sharp peaks in figs . [ fig5]c and [ fig5]d are due to recurrences of the density matrix @xmath54 to @xmath53 , we define the figure - of - merit of recurrence , @xmath181 where the prime means the summation is over @xmath182 such that @xmath183 . it is clear that @xmath184 and @xmath185 when and only when all the off - diagonal elements get in phase .    in figs . [ fig6]a and [ fig6]b , which share the same parameters as figs . [ fig5]c and [ fig5]e respectively , we have shown @xmath186 together with @xmath171 and @xmath174 . in fig . [ fig6]a , we see that every time @xmath171 and @xmath172 get close to their values at @xmath187 , @xmath186 shows a peak . in other words , there is a strong positive correlation between @xmath186 and @xmath171 and @xmath172 . in comparison , in fig . [ fig6]b , @xmath186 drops quickly from unity to less than @xmath188 and remains low all the time , and in turn @xmath171 and @xmath189 do not show any recurrence . to further consolidate the connection between the recurrence of @xmath54 and that of @xmath171 and @xmath174 , we have considered the case of @xmath190 . in this case , if @xmath60 , all the basic frequencies @xmath180 are commensurate , and thus there exist perfect recurrences , as shown in fig . there we see clearly that @xmath171 and @xmath172 return to their original values at @xmath187 periodically , and this happens when and only when @xmath191 returns to unity . however , once @xmath192 is set nonzero ( see fig .  [ fig6]d ) and thus the commensurability of the energy gaps is destroyed , the situation returns to that in fig . [ fig6]b . the fact revealed in fig . [ fig5 ] and fig . [ fig6 ] is quite interesting . the long - time dynamics of the system is sensitive or insensitive to the exact time when the second quench is applied , depending on whether the intermediate hamiltonian @xmath40 is integrable ( @xmath60 ) or non - integrable ( @xmath55 ) . in the integrable case , @xmath172 exhibits large fluctuations and repeated recurrences . the system retains the memory of the initial state under the control of the hamiltonian @xmath40 . by contrast , in the non - integrable case , @xmath193 go over to their average values ( predicted by @xmath166 ) after a transitory period , showing little dependence on @xmath159 afterwards . combined with fig . [ fig4 ] , the picture is that @xmath54 evolving under the control of a non - integrable hamiltonian , not only yields the expectation values of @xmath194 as if it were @xmath51 , but even responds to the second quench as if it were @xmath51 .    in fig . [ fig7 ] , we have checked this picture by studying the real time evolution of @xmath195 with @xmath196 under the double - quench scenario . the eight figures shown correspond to those in fig .  [ fig5 ] respectively . for each pair of @xmath173 , we have studied the evolution of @xmath197 for several different values of @xmath159 . we see that in all the cases with @xmath55 , as long as @xmath159 is larger than the transient time , which can be roughly read from fig . [ fig5 ] , the later evolution of @xmath198 is quantitatively independent of @xmath159 . on the contrary , in the case with @xmath60 , the later values of @xmath199 vary wildly for different values of @xmath159 . here it is instructive to combine fig . [ fig4 ] and fig . [ fig7 ] and compare . in the @xmath55 cases , there is a sense of typicality @xcite . the density matrix @xmath54 governed by @xmath40 is surely non - stationary . however , for @xmath54 at different times , they yield almost the same expectation values for the observables , and moreover , they share almost the same response to the same quench . in the case of @xmath60 , what fig . [ fig7 ] reveals is a good complement to that in fig . it demonstrates that it is inappropriate to say that the system thermalizes in this case , even though the density matrices and expectation values of the observables agree  since according to one s everyday experience , a system in thermal equilibrium should not show any time dependence . we have studied the quench dynamics of the bose - hubbard model both analytically and numerically . the issues of thermalization  and equilibration are investigated thoroughly .    on the thermalization  side , which concerns whether the quenched system behaves like a canonical ensemble , it is found that this is the case only for small - amplitude quenches ( at least for the finite - sized system investigated ) . however , the time - averaged density matrix does manifest many interesting features in different regimes . these features are self - consistently understood after a study of the overlaps between the eigenstates of @xmath12 and @xmath40 . here we would like to say that it is urgent and would be very helpful to develop some analytical tools so that some general relations between the eigen - systems of @xmath12 and @xmath40 can be established . these tools and relations would also be useful to determine whether the non - thermalizationphenomenon observed is just a finite - size effect .    on the equilibration side , which is about whether physical observables relax to stationary values without appreciable fluctuations , the result is that this is indeed the case for quantities as @xmath200 which are of most interest . moreover , it is proven analytically that for these quantities the fluctuations in time will decay exponentially with the size of the system . therefore , the overall picture is that generally the system equilibrates but without thermalization . the second quench reveals something more subtle . first , the subsequent dynamics depends or not on the second quench time @xmath159 according to @xmath60 or not . the underline reason is the recurrence or not of the initial density matrix , which in turn has its root in the eigenvalue statistics of the hamiltonian @xmath40 . this effect leaves us the impression that a non - integrable hamiltonian has more `` dephasing power '' than an integrable one . possibly it can be a tool to check the integrability of a hamiltonian . second , in the case of @xmath55 , it is found that the system described by @xmath154 responds to the second quench as if it were @xmath51 for @xmath159 larger than the transient time . this means that we can take the equilibration more serious@xmath154 and @xmath51 not only yield almost the same expectation values for the generic physical variables but also yield almost the same dynamics after a quench . moreover , the fact that the transient time is short indicates that the intermediate hamiltonian @xmath40 , which is non - integrable , is effective in `` dephasing '' the initial density matrix . in another perspective , the dynamics of the system is sensitive to the fluctuations of @xmath0 . this has the implication that in future experiments , accurate control of @xmath0 would be a necessity to interpret the results correctly . we are grateful to h.  t. yang , z.  x. gong , y .- h . chan , l.  m. duan , and d.  l. zhou for stimulating discussions and valuable suggestions . j. m. z. is supported by nsf of china under grant no .
we study quench dynamics of the bose - hubbard model by exact diagonalization . 
 initially the system is at thermal equilibrium and of a finite temperature . 
 the system is then quenched by changing the on - site interaction strength @xmath0 suddenly . both the single - quench and double - quench scenarios are considered . in the former case , 
 the time - averaged density matrix and the real - time evolution are investigated . 
 it is found that though the system thermalizes only in a very narrow range of the quenched value of @xmath0 , it does equilibrate or relax well in a much larger range . 
 most importantly , it is proven that this is guaranteed for some typical observables in the thermodynamic limit . in order to test whether it is possible to distinguish the unitarily evolving density matrix from the time - averaged ( thus time - independent ) , fully decoherenced density matrix , 
 a second quench is considered . 
 it turns out that the answer is affirmative or negative according to the intermediate value of @xmath0 is zero or not .