How to calculate quantum quench distributions with a weighted Wang-Landau Monte Carlo
HHow to calculate quantum quench distributionswith a weighted Wang-Landau Monte Carlo
Simone Ziraldo , and Giuseppe E. Santoro , , SISSA, Via Bonomea 265, I-34136 Trieste, Italy CNR-IOM Democritos National Simulation Center, Via Bonomea 265, I-34136Trieste, Italy International Centre for Theoretical Physics (ICTP), P.O.Box 586, I-34014Trieste, Italy
Abstract.
We present here an extension of the Wang-Landau Monte Carlomethod which allows us to get very accurate estimates of the full probabilitydistributions of several observables after a quantum quench for large systems,whenever the relevant matrix elements are calculable, but the full exponentialcomplexity of the Hilbert space would make an exhaustive enumeration impossiblebeyond very limited sizes. We apply this method to quenches of free-fermionmodels with disorder, further corroborating the fact that a Generalized GibbsEnsemble fails to capture the long-time average of many-body operators whendisorder is present.PACS numbers: 05.70.Ln, 75.10.Pq , 72.15.Rn, 02.30.Ik a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a y ONTENTS Contents1 Introduction 22 Weighted Wang-Landau algorithm 43 Quantum quenches 5
A sudden quench of the Hamiltonian parameters is perhaps the simplest form of out-of-equilibrium dynamics that a closed quantum system can experience. Experimentswith “virtually isolated” cold atomic species in optical lattices [1, 2] have transformedthis seemingly theoretical dream into a rich and lively stage. Several fundamentalissues of theoretical quantum statistical physics, like the onset of thermalization whichis generally expected to occur for a closed quantum system after a sudden quench[3, 4, 5], or the “breakdown of thermalization” [6, 7] expected when the system isintegrable or nearly integrable, are now of experimental relevance [8]. We refer thereader to a recent review [9] for an extensive introduction to such non-equilibriumquantum dynamics issues.The issue we want to tackle in this paper is the following. Suppose you performa quantum quench of the Hamiltonian parameters, abruptly changing, at t = 0, fromˆ H → ˆ H . If | Ψ (cid:105) denotes the initial quantum state at t = 0, and | α (cid:105) the eigenstatesof ˆ H with energy E α , the ensuing quantum dynamics would lead to averages for anygiven operator ˆ A given by: A ( t ) ≡ (cid:104) Ψ | e i ˆ Ht ˆ Ae − i ˆ Ht | Ψ (cid:105) = (cid:88) α | c α | A αα + (cid:88) α (cid:48) (cid:54) = α e i ( E α (cid:48) − E α ) t c ∗ α (cid:48) A α (cid:48) α c α , (1)where c α ≡ (cid:104) α | Ψ (cid:105) and A α (cid:48) α ≡ (cid:104) α (cid:48) | ˆ A | α (cid:105) . The first (time-independent) term in theprevious expression dominates the long-time average of A ( t ), and is usually known as diagonal average [5] (cid:104) ˆ A (cid:105) D ≡ (cid:88) α | c α | A αα . (2)To calculate it, in principle, we should take the sum over all the (many-body)eigenstates | α (cid:105) of ˆ H — an exponentially large number of states —, calculating foreach of them the overlap c α and the associated diagonal matrix element A αα . Luckily, (cid:104) ˆ A (cid:105) D can be calculated for many problems, notably those that can be reduced toquadratic fermionic problems, by circumventing in one way or another exponentiallylarge sums: for instance, through a detour to time-dependent single-particle Green’sfunction and the use of Wick’s theorem, see e.g. [10, 11]. But suppose that you wantto know more than just the diagonal average (cid:104) ˆ A (cid:105) D , and pretend to have informationon the whole distribution of the values of A αα accessed after the quench [12, 13], i.e., ρ D ( A ) ≡ (cid:88) α | c α | δ ( A − A αα ) , (3) ONTENTS (cid:104) ˆ A (cid:105) D is just the average: (cid:104) ˆ A (cid:105) D = (cid:82) d A ρ D ( A ) A . Here there is, evidently, aproblem: knowing the distribution of A requires exploring the full many-body Hilbertspace, summing over the eigenstates | α (cid:105) , and this exhaustive enumeration wouldrestrict our calculations to exceedingly small sample sizes, although all informationon c α and A αα might in principle be easy to calculate, or in any case accessible, forinstance by just solving a one-body problem (hence, for much larger sizes). A similarproblem occurs in considering, for instance, the corresponding ‡ generalized Gibbsensemble (GGE) average [14, 15, 16, 17, 18] (cid:104) ˆ A (cid:105) GGE ≡ (cid:88) α e − (cid:80) µ λ µ I αµ Z GGE A αα , (4)where λ µ are Lagrange multipliers which constrain the mean value of each of theconstants of motion ˆ I µ to their t = 0 value, (cid:104) Ψ | ˆ I µ | Ψ (cid:105) = Tr (cid:104) ˆ ρ GGE ˆ I µ (cid:105) , I αµ ≡ (cid:104) α | ˆ I µ | α (cid:105) , Z GGE is the GGE partition function, and ˆ ρ GGE ≡ e − (cid:80) µ λ µ ˆ I µ /Z GGE . Once again, for“quadratic problems” this average is rather simply calculated in terms of single-particlequantities, but the corresponding distribution ρ GGE ( A ) ≡ (cid:88) α e − (cid:80) µ λ µ I αµ Z GGE δ ( A − A αα ) , (5)requires a difficult sum over the Hilbert space. § Concerning the issue of thermalization after a quantum quench, we might indeedexpect that, if the system is well described by a GGE ensemble, not only the meanvalues of ρ D ( A ) and ρ GGE ( A ) are equal, i.e., (cid:104) ˆ A (cid:105) D ≡ (cid:82) d A ρ D ( A ) A = (cid:104) ˆ A (cid:105) GGE ≡ (cid:82) d A ρ
GGE ( A ) A , but also the two distributions should be closely related; at least thisis what a good statistical ensemble should do.Quite generally, we might formulate the problem as follows: how can we obtaininformation on weighted distributions (or density of states) ρ w ( A ) ≡ (cid:88) α w α δ ( A − A α ) , (6)with positive weights w α , when both w α and A α are “easily calculated”, but thesum over α runs over an exponentially large “configuration space”? As discussedbefore, examples of this are the diagonal distribution, where w α = | c α | , the GGEdistribution, where w α = e − (cid:80) µ λ µ I αµ /Z GGE , but also the microcanonical distribution,where w α is a window characteristic function for the microcanonical shells, etc. AMonte Carlo algorithm to perform such exponentially large sums in configuration spaceseems unavoidable. We stress that this is so even if one is considering quenches infree-fermion models, where the relevant many-particle states | α (cid:105) and matrix elementsare easy to write down and calculate. The alternative of using exact diagonalizationmethods would put a strong limit on the size of the problem which can be studied.In this paper we introduce a Monte Carlo method — obtained by a rather naturalextension of the Wang-Landau algorithm (WLA) [19, 20, 21] — which will allow us tocompute weighted distributions of the form of Eq. (6). The Wang-Landau algorithm, ‡ The generalized Gibbs ensemble is the relevant ensemble for quenches with quadratic fermionicmodels, but a similar situation occurs for the usual statistical ensembles, i.e., microcanonical,canonical and grand-canonical. § For the GGE ensemble things might be worked out by an appropriate representation of the Dirac’sdelta, or by computing the moment generating function of ρ GGE ( A ). These tricks however would notwork, for instance, for the microcanonical distribution. ONTENTS
2. Weighted Wang-Landau algorithm
Let us consider a system with a discrete configuration space, where configurations canbe labeled with an index α . Given a physical observable ˆ A , we define its weighted(coarse-grained) density of states: ρ w ( A ) ≡ (cid:88) α w α δ AA α , (7)with w α a positive weight. Here δ AA α is a Kronecker delta, or, if the possible valuesof A α are too dense to keep them all, a suitable histogram-window-function coarse-graining of the Dirac’s delta. When w α = 1, we recover the usual density of states ρ ( A ), and the WLA can be used to estimate it [19, 20, 21]. We will now show that,by properly modifying the WLA, we can compute ρ w ( A ) for generic w α s.To understand the gist of the approach, consider a generic positive function ˜ ρ w ( A )— which is our best guess for the desired ρ w ( A ) —, and set up a Markov chain randomwalk in which, given a state α , a new state α (cid:48) is generated with a trial probability T ( α (cid:48) | α ), which we will take to be symmetric, T ( α (cid:48) | α ) = T ( α | α (cid:48) ), and accepted withprobability: R ( α (cid:48) | α ) = Min (cid:20) , w α (cid:48) w α ˜ ρ w ( A α )˜ ρ w ( A α (cid:48) ) T ( α | α (cid:48) ) T ( α (cid:48) | α ) (cid:21) . (8)With this standard Metropolis Monte Carlo prescription, we know that, after aninitial transient, we will visit the configurations α with an equilibrium distribution P eq α fulfilling the detailed balance condition and given by: P eq α = C w α ˜ ρ w ( A α ) , where C is a normalization constant. As in the WLA [19], while the random walkgoes on, we collect a histogram h ( A ), updating h ( A α ) → h ( A α ) + 1 at each visitedstate α . At equilibrium, after N s steps, the “mean” histogram will then be given by: h ( A ) = N s (cid:88) α P eq α δ AA α = N s C (cid:88) α w α ˜ ρ w ( A α ) δ AA α = N s C ρ w ( A )˜ ρ w ( A ) . (9)Exactly as for the WLA [19], if our guess for ˜ ρ w ( A ) is a good approximation to ρ w ( A ),the histogram h ( A ) will be “almost flat” (see below). Obviously, during the randomwalk, together with the histogram h ( A ) we also update our guessed ˜ ρ w ( A ). Therefore,closely inspired by the WLA [19], we propose the following algorithm: ONTENTS f >
1, and set ln ˜ ρ w ( A ) = 0 and h ( A ) = 0 for all valuesof A ;(1) Start the Monte Carlo procedure using Eq. (8) and update at each step thehistogram and the weighted density of states with the rules h ( A α ) → h ( A α ) + 1and ln ˜ ρ w ( A α ) → ln ˜ ρ w ( A α ) + ln f ;(2) Stop the random walk when h ( A ) is “almost flat” (for instance [19], when h ( A ) > . h for all values of A , where h is the mean histogram value). For theprevious observations, at the end of this step ln ˜ ρ w ( A ) is a good approximation toln ρ w ( A ) with a discrepancy of order ln f ;(3) Reduce the value of f → √ f , reset h ( A ) = 0 and restart the procedure from step(1) using the ˜ ρ w ( A ) just obtained. Stop this loop when ln f is smaller than thedesired discrepancy (cid:15) .A similar extension of the WLA has been already been introduced for the particularcase in which w α is the Boltzmann distribution, with the aim of computing the freeenergy profile as a function of a reaction coordinate [24, 25]. In the present paper, wewill use this algorithm to compute distributions where the weights are not Boltzmann-like, but rather associated to quantum quenches.Let us return for a moment to the original WLA. A first trivial observation isthat, as it should be, the weighted-WLA with w α = 1 /N coincides with the WLA. Inmany situations, when the size of the configuration space is too big and the density ofstates ranges over too many orders of magnitude, it is convenient, in computing ρ ( A ),to run many WLA over small domains ∆ ( i ) A = [ A ( i )min , A ( i )max ]. But then the updaterule of the standard WLA has to be changed to avoid that, during the random walk, A α leaves the domain ∆ ( i ) A . This trick was already used in the first papers by Wangand Landau, when dealing with the largest sizes [19]. To avoid leaks from ∆ ( i ) A , theempirical solution was to reject any proposal to states α (cid:48) with A α (cid:48) / ∈ ∆ ( i ) A , withoutany update of ˜ ρ ( A ) and h ( A ). With this prescription, however, there are “boundaryeffects”, actually a systematic underestimation of the density of states at the bordersof the intervals [26]. Schulz et al. [26] showed phenomenologically that such boundaryeffects are eliminated by using the rather obvious update rule: given a proposal α (cid:48) , if A α (cid:48) is outside the interval we remain in α and we update h ( A ) and ln ˜ ρ ( A ) using thestate α , otherwise we accept α (cid:48) with the usual rule. This update rule is just what isobtained, rigorously, by using our weighted-WLA. Indeed, the density of states in arestricted range ∆ ( i ) A is proportional to a weighted density of states in which w α = 1when A α ∈ ∆ ( i ) A , and zero otherwise. With these weights, the update rule of ourweighted-WLA is exactly the one obtained phenomenologically by Schulz et al. [26].
3. Quantum quenches
In this section we come back to the initial problem of computing the distributions ρ D ( A ) and ρ GGE ( A ) related to quantum quenches. We will show that with theweighted-WLA we can compute these distributions for sizes inaccessible with anexhaustive enumeration.We concentrate on quantum quenches in two models possessing a free fermionicdescription. The first model we considered is the fermionic Anderson model with ONTENTS H A ≡ − t L (cid:88) j =1 (cid:16) ˆ c † j ˆ c j +1 + h . c . (cid:17) + L (cid:88) j =1 h j ˆ c † j ˆ c j , (10)where ˆ c † j (ˆ c j ) creates (destroys) a fermion at site j , t is the nearest-neighbor hoppingintegral and h j is an uncorrelated on-site random potential uniformly distributedin the range [ − W/ , W/ H A and in presence of any W >
0, all the single-particle eigenstates of ˆ H A are exponentially localized. Thesecond model we considered still describes spinless fermions hopping on a chain, butnow the hopping is long-ranged [28]:ˆ H lrh = (cid:88) j j t j j (ˆ c † j ˆ c j + h . c . ) , (11)where t j j is a (real) hopping integral between sites j and j . We will take the t j j ’sto be random and long-ranged, with a Gaussian distribution of zero mean, (cid:104) t j j (cid:105) = 0,and variance given by: (cid:104) t j j (cid:105) = 11 + (cid:16) | j − j | β (cid:17) γ . (12)Here γ is a real positive parameter setting how fast the hoppings’ variance decayswith distance. Notice that, for j = j , we have (cid:104) t j j (cid:105) = 1 for any γ , hence themodel has also on-site Gaussian disorder; by increasing the distance between the twosites | j − j | , the variance of the hopping integral decreases with a power law. Thepeculiarity of this long-range-hopping model is that, although being one-dimensionaland regardless of the value of β (which hereafter is fixed to 1), it has an Andersontransition from (metallic) extended eigenstates, for γ <
1, to (insulating) power-lawlocalized eigenstates for γ > γ , long-range hoppings are capable of overcoming the localization due todisorder. Having access, in the same model, to physical situations in which the finaleigenstates are extended ( γ <
1) or localized ( γ >
1) will clearly show the role thatspatial localization plays in disrupting the ability of the GGE to describe the after-quench dynamics. Physically, spatial localization prevents the different “modes” ofthe system from having an infinite reservoir.For the considered quenches, we use as initial Hamiltonian ˆ H the clean chainwith W = 0 and the same boundary conditions of the final Hamiltonian, i.e., periodicboundary conditions when quenching to ˆ H A and open boundary conditions whenquenching to ˆ H lrh . The corresponding initial state | Ψ (cid:105) will be the filled Fermi sea,i.e., the ground state of ˆ H with N F = L/
2, where N F is the number of fermions. Thereason behind this simple choice for ˆ H is that the “stationary state” reached does notdepend, qualitatively, on the initial Hamiltonian being ordered or not, see Ref. [10].The final Hamiltonian will be the Anderson model ˆ H A with W = 2, or the long-rangehopping chain ˆ H lrh with γ = 0 . N F = L/ t >
0. To get a smoother size dependenceof the computed quantities, the smaller size realizations are obtained by cutting anequal amount of sites at the two edges of the largest realization.
ONTENTS L in terms of new fermionic operatorsˆ d † µ = L (cid:88) j =1 u jµ ˆ c † j , (13)where u jµ are the wave functions of the eigenmodes of energy (cid:15) µ : ˆ H A / lrh = (cid:80) µ (cid:15) µ ˆ d † µ ˆ d µ .The energies (cid:15) µ and the associated wave functions u jµ are obtained, for any givendisorder realization of a chain of size L by numerically diagonalizing the L × L one-body hopping matrix.Given an observable ˆ A , consider the two distributions introduced before: ρ D ( A ) ≡ (cid:88) α | c α | δ ( A − A αα ) (14) ρ GGE ( A ) ≡ (cid:88) α e − (cid:80) µ λ µ n αµ Z GGE δ ( A − A αα ) , (15)where δ ( x ) is the Dirac’s delta, {| α (cid:105)} are the many-body eigenstates of ˆ H , A αβ ≡(cid:104) α | ˆ A | β (cid:105) , c α ≡ (cid:104) α | Ψ (cid:105) , and n αµ = 0 , µ in the many-body eigenstate | α (cid:105) . These functions give the weighted distributions of A αα in the diagonal and GGE ensembles.Let us discuss a few technical details on the implementation we made, beforediscussing the physics emerging from our calculations. Notice that the sum over α is effectively restricted to the canonical Hilbert space H N with a fixed number ofparticles N = N F in the diagonal ensemble, since c α ≡ (cid:104) α | Ψ (cid:105) = 0 if N α (cid:54) = N F .No such restriction is in principle present in the GGE case, where the sum over α runs over the grand-canonical Hilbert space. By definition, the distributions are suchthat (cid:104) ˆ A (cid:105) D = (cid:82) A ρ D ( A ) d A and (cid:104) ˆ A (cid:105) GGE = (cid:82) d A ρ
GGE ( A ) A , where the integration isover the domain of A αα . As customary in any numerical finite-size study, one reallyneeds to consider a coarse-grained version of these distributions, obtained by splittingthe domain of A into small intervals ∆ ( i ) of amplitude ∆. Such a coarse-graineddistribution has exactly the form of a weighted density of states, see Eq. (7), with w α = | c α | / ∆ in the diagonal case, and w α = e − (cid:80) µ λ µ n αµ / (∆ Z GGE ) in the GGE case.The configuration space {| α (cid:105)} (i.e., the canonical Hilbert space H N for the diagonaldistribution and the full Hilbert space for the GGE) over which the two weighteddistributions are defined is discrete and grows exponentially with the system size.The weighted-WLA is therefore the appropriate tool for the numerical computationof ρ D ( A ) and ρ GGE ( A ). The eigenstates | α (cid:105) which appear in the definition of ρ D ( A )have a fixed number of fermions N F (the same of the initial state), while in ρ GGE ( A )the number of particles can change. In the weighted-WLA, for the diagonal ensemble,we use therefore a “particle conserving” proposal scheme: given a state | α (cid:105) , the state | α (cid:48) (cid:105) is given by moving at random a fermion in one of the unoccupied single-particleeigenstates. In this case, the ratio w α (cid:48) /w α which appears in Eq. (8), is equal to: w α (cid:48) w α = | c α (cid:48) | | c α | , where the coefficient | c α | is the square of the determinant of a N F × N F matrix (see[31, App. D] for the explicit expression of | c α | ). For the GGE case, instead, we donot have restrictions on the number of fermions and, given a state | α (cid:105) , we generate a ONTENTS | α (cid:48) (cid:105) by changing the occupation of a randomly selected single-particle eigenstate µ . In this case: w α (cid:48) w α = e ± λ µ , where the + ( − ) sign appears when the mode µ is initially occupied (empty). Let usrecall that the Lagrange’s multipliers λ µ are obtained by requiring (cid:104) Ψ | ˆ d † µ ˆ d µ | Ψ (cid:105) = (cid:104) ˆ d † µ ˆ d µ (cid:105) GGE . This condition, written explicitly, reads:e λ µ = 1 (cid:80) ν n ν (cid:12)(cid:12)(cid:12) [ u † u ] νµ (cid:12)(cid:12)(cid:12) − , (16)where u and u are L × L matrices whose elements u jν and u jµ are the single-particlewavefunctions of the initial Hamiltonian ˆ H and the final one ˆ H , and n ν = 0 , ν th eigenstate of ˆ H in the initial state. The difference in thecomputational effort on computing the ratio w α (cid:48) /w α in the two ensembles is evident:in the diagonal case at each step we have to compute the determinant of a N F × N F matrix, while in the GGE we have just to recover the value of e λ µ (they can becomputed and stored before the Monte Carlo calculation because their number is L ).Here we will show results for sizes up to L = 256, where both ρ D ( A ) and ρ GGE ( A )can be computed and compared. (For the GGE ensemble, we could reach L = 1024without problem.) In the numerical computations we used a minimum value of theWL parameter (cid:15) = ln f min = 10 − , and we split the domain of A in L bins. Noticethat the domain of A in ρ GGE ( A ) is always larger than the domain of ρ D ( A ) because,in the GGE, the many-body eigenstates do not have a restriction on the number N F of fermions.In the next two subsections we show the results obtained with the weighted-WLA for the calculation of ρ D ( A ) and ρ GGE ( A ) for two observables, the total energyand the local density. The physical picture emerging from the calculation of the fulldistribution function of the after-quench energy and local-density confirms and extendsthe results discussed in Ref. [10, 11]. In particular, we find clear differences betweenthe diagonal and GGE distributions, even at the level of the variances, whenever adisorder-induced spatial localization is at play in the after-quench Hamiltonian. The first observable we consider is the total energy: Here A αα → E α = (cid:80) µ (cid:15) µ n αµ ,where n αµ = (cid:104) α | ˆ d † µ ˆ d µ | α (cid:105) = 0 , | α (cid:105) . In Fig. 1 we show the distributions ln[ ρ D ( E )] /L and ln[ ρ GGE ( E )] /L , computedfor L = 128 and L = 256, for the three cases we have studied, i.e., quenches from aninitially ordered half-filled chain ˆ H towards: 1) a long-range hopping Hamiltonianˆ H lrh with extended eigenstates ( γ = 0 .
5, top), 2) ˆ H lrh with localized eigenstates( γ = 2, center), and 3) an Anderson Hamiltonian ˆ H A with a disorder width W = 2(bottom). Observe, first, that the distributions ρ D ( E ) and ρ GGE ( E ) shown in Fig. 1have identical average (denoted by a solid vertical line) (cid:104) ˆ H (cid:105) D = (cid:90) d E ρ D ( E ) E = (cid:90) d E ρ
GGE ( E ) E = (cid:104) ˆ H (cid:105) GGE . This result comes directly from the fact that the energy does not fluctuate in time(i.e., the diagonal energy coincides with the average energy (cid:104) Ψ | ˆ H | Ψ (cid:105) ) and GGE ONTENTS Figure 1.
Value of ln[ ρ D ( E )] /L and ln[ ρ GGE ( E )] /L computed with the weighted-WLA. The gray curves are obtained with L = 128, while the black ones with L = 256. The solid vertical lines are the average energy after the quench, i.e. (cid:104) Ψ | ˆ H | Ψ (cid:105) , for L = 256. The three panels are obtained starting from the groundstates of clean chains and quenching to different disordered Hamiltonians: panel(a) long-range hopping with γ = 0 . γ = 2 (localized eigenstates) and panel (c) Anderson model with W = 2. For the computations we used a single disorder realization and, to get asmoother size dependence, the smaller size realization is obtained by cutting anequal amount of sites at the two edges of the larger realization. These distributionsare obtained for a single realization of the couplings, but we verified that, for largesizes, the results are self-averaging. fixes the occupation of the fermionic eigenstates in such a way as to exactly reproduce (cid:104) Ψ | ˆ H | Ψ (cid:105) . The form of the two distributions, however, differs considerably, mostnotably at the extremes of the spectrum, and for the Anderson model case. Let usnow consider the fluctuations of the energies in both distributions. In the diagonalensemble the variance is: σ E, D = (cid:90) d E ρ D ( E ) E − (cid:104) ˆ H (cid:105) D = (cid:104) ˆ H (cid:105) D − (cid:104) ˆ H (cid:105) D , (17)where the expression on the right-hand side holds only for the Hamiltonian (it wouldnot apply to arbitrary operators, because ( A αα ) (cid:54) = (cid:104) α | ˆ A | α (cid:105) ). An entirely similarexpression applies to the GGE case. Since the energy is an extensive operator, it is ONTENTS e = E/L ,which are simply given by σ e, D = σ E, D /L , and σ e, GGE = σ E, GGE /L . On pretty generalgrounds, for quenches of local non-integrable Hamiltonians, it is known [12, 13] that σ e, D → L → ∞ . Indeed, as shown in Fig. 2 both σ e, D and σ e, GGE decrease to 0 for L → ∞ for the three considered cases. For our quadraticproblems, however, we can say a bit more. First of all, from the explicit expression inEq. (17) after very simple algebra (mainly using Wick’s theorem), we arrive at: σ e, GGE = 1 L (cid:88) µ (cid:15) µ n µ (cid:0) − n µ (cid:1) , (18) σ e, D = σ e, GGE − L (cid:88) µ (cid:54) = µ (cid:15) µ (cid:15) µ (cid:12)(cid:12) G µ µ (cid:12)(cid:12) , (19)where G µ µ ≡ (cid:104) Ψ | ˆ d † µ ˆ d µ | Ψ (cid:105) is the t = 0 one-body Green’s function. The off-diagonal elements of G µ µ play here an important role, and the second term in σ e, D originates from the fact that, by definition, GGE does not include correlations betweendifferent eigen-modes, i.e., (cid:104) ˆ d † µ ˆ d µ (cid:105) GGE = 0, when µ (cid:54) = µ . Figure 2.
Variances σ e, D = σ E, D /L (empty circles) and σ e, GGE = σ E, GGE /L (solid triangles) as a function of the size L . The data are obtained using the sameset of quenches used in Fig. 1 and the values are computed using Eqs. (19). Errorbars are calculated by averaging over 20 different realizations of the disorder.The dashed lines are power law fits σ e ∼ L − s , where s ≈ H lrh , s ≈ .
82 when γ = 0 .
5, and s ≈ .
95 when γ = 2. Notice the observable difference between σ e, D and σ e, GGE when the finaleigenstates are localized.
Let us first consider the Anderson model case. Assuming, as done so far, abounded distribution of disorder, we are guaranteed that a finite bound (cid:15) max exists
ONTENTS | (cid:15) µ | ≤ (cid:15) max for any L . With this assumption, it is easy show that σ e, GGE hasto go to zero at least as 1 /L for L → ∞ . Indeed, the occupation factors appearing in σ e, GGE are such that 0 ≤ n µ (cid:0) − n µ (cid:1) ≤ /
4. Hence: σ e, GGE ≤ (cid:15) L (cid:88) µ n µ (cid:0) − n µ (cid:1) ≤ (cid:15) L . (20)The same statement can be made for σ e, D , because the difference between the twovariances has a similar upper bound: | σ e, D − σ e, GGE | ≤ L (cid:88) µ (cid:54) = µ | (cid:15) µ || (cid:15) µ | (cid:12)(cid:12) G µ µ (cid:12)(cid:12) ≤ (cid:15) L (cid:88) µ (cid:54) = µ (cid:12)(cid:12) G µ µ (cid:12)(cid:12) ≤ (cid:15) L , (21)where we used that (cid:80) µ µ (cid:12)(cid:12) G µ µ (cid:12)(cid:12) = N F ≤ L . Nevertheless, although both σ e, D and σ e, GGE go to 0 as 1 /L for the Anderson model, they do so with a different pre-factor,see Fig. 2 and comments below.For the quenches to ˆ H lrh , a bound (cid:15) max for the single-particle spectrum is inprinciple not defined: one can think of rare realizations in which the hopping is large atarbitrarily large distances, which would give an unbounded distribution of eigenvalues (cid:15) µ . Indeed, the behavior of both σ e, D and σ e, GGE suggest, see Fig. 2, that the power-law approach to 0 might be slower than 1 /L , i.e., as L − s with s < s ≈ . γ = 0 . s ≈ .
95 for γ = 2). While this might be a finite-size artifact,we find it intriguing that such deviations are quite clearly seen for quenches to ˆ H lrh :they might be due to the power-law nature of the hopping integral variance.Concerning the similarity between σ e, D and σ e, GGE , we observe that the twoessentially coincide for the case of a quench to ˆ H lrh with extended eigenstates, whilethere is a small discrepancy for the quench to ˆ H lrh with localized eigenstates, anda quite clear different pre-factor in the Anderson model case, σ e, D ∼ C D /L and σ e, GGE ∼ C GGE /L with C GGE < C D . This different pre-factor can be understood byanalyzing the term (cid:80) µ (cid:54) = µ (cid:15) µ (cid:15) µ | G µ µ | which appears in Eq. (19). In Fig. 3, panel(b), we show the structure of the matrix | G µ µ | for the three quench cases. We dividethis matrix into four sectors, one for each sign of the single-particle energies (cid:15) µ and (cid:15) µ : in two of these quadrants the product (cid:15) µ (cid:15) µ is positive (top-right and bottom-left), and in the others is negative. For quenches to ˆ H lrh this matrix is almost equallydistributed in all the four sectors: the sum (cid:80) µ (cid:54) = µ (cid:15) µ (cid:15) µ | G µ µ | has cancellations,leading to σ e, GGE ≈ σ e, D for large sizes. For quenches to ˆ H A , on the contrary, thematrix | G µ µ | is mainly concentrated in the sectors in which (cid:15) µ (cid:15) µ <
0, leading to σ e, GGE < σ e, D .Finally, let us comment on one aspect of the distributions shown in Fig. 1 whichcan be easily understood from the single-particle occupations shown in Fig. 3. We seethat, when the after-quench Hamiltonian is the Anderson model, ρ D ( E ) has both mode(i.e., maximum value) and average very close to the ground state energy: the quenchexcites mostly the low-energy part of the many-body spectrum. On the contrary, forboth the quenches towards ˆ H lrh , mode and average are almost in the middle of themany-body spectrum; there, indeed, the quench is more dramatic: we are going fromthe ground state of a chain with nearest-neighbor hopping to a disordered chain withlong-range hopping. This is evident by looking at the occupations n µ ≡ (cid:104) Ψ | ˆ d † µ ˆ d µ | Ψ (cid:105) as a function of the single-particle energy (cid:15) µ , shown in Fig. 3, panel (a). By definition,only the eigenstates of ˆ H with (cid:15) ν < | Ψ (cid:105) . The quench to ˆ H A only ONTENTS
Figure 3.
Panel (a): occupations n µ = (cid:104) Ψ | ˆ d † µ ˆ d µ | Ψ (cid:105) as a function of thesingle-particle energy (cid:15) µ . Panel (b): representation of the matrix | G µ µ | . Forthe diagonal and off-diagonal elements we add a black pixel when the value exceedstheir mean value. For the diagonal elements the mean value is x ≡ (cid:80) µ ( n µ ) /L ,while for the off-diagonal elements the mean value is ( N F − xL ) /L ( L − N F is the number of fermions in the initial state, and we used the relation (cid:80) µ µ | G µ µ | = N F (see [31, App. D]). The vertical and horizontal linesindicate the indexes at which the single-particle energies (cid:15) µ and (cid:15) µ change sign,and the signs shown in the four quadrants are those of the product (cid:15) µ (cid:15) µ . Forthe two panels we used L = 256 and the same quenches used in Fig. 1 and Fig. 2. slightly modifies the initial occupations: n µ , apart for fluctuations due to disorder, goessmoothly from 1, in the lower part of the single-particle spectrum, to 0, in the highestpart of the spectrum. On the contrary, for the quenches towards ˆ H lrh , the single-particle spectrum is entirely excited, both the positive and the negative energy part.This explains why, for these quenches, the after-quench energy (cid:104) Ψ | ˆ H | Ψ (cid:105) = (cid:80) µ (cid:15) µ n µ is near the center of the many-body spectrum. Let us now consider the local density ˆ n j ≡ ˆ c † j ˆ c j , perhaps the simplest one-bodyobservable. For definiteness, we concentrate on j = L/
2, the center of the chain.It is important to stress that we are going to consider the fluctuations of ˆ n j before anypossible average over the sites j : averaging over the sites j an intensive local operatorwould effectively send to zero the fluctuations in the thermodynamic limit [13], whilewe will show that, for a fixed j , finite fluctuations survive in the thermodynamic limitwhen the eigenstates are localized, due to disorder.The diagonal and GGE distributions ρ D ( n ) and ρ GGE ( n ) are now constructed usingthe matrix elements n αα ≡ (cid:104) α | ˆ n j | α (cid:105) = (cid:80) µ | u jµ | n αµ , where n αµ = 0 , | α (cid:105) . In Fig. 4 we plot ln[ ρ D ( n )] and ONTENTS ρ GGE ( n )], computed for the three quenches discussed before. The case of a quenchto ˆ H lrh with γ = 2 . n canassume is actually split in two separated domains, one just above n = 0 and one justbelow n = 1, and the mean value is exactly in the middle, where no values of n αα happen to fall. This is due to the strong spatial localization of the eigenstates. As weshow in Fig. 5, at fixed j , the value of | u jµ | is strongly localized in a single eigenstate˜ µ . This implies that the value n αα = (cid:80) µ | u jµ | n αµ has a strong jump when we movefrom a state | α (cid:105) in which n α ˜ µ = 0, to the state | α (cid:105) in which n α ˜ µ = 1. For the quenchto ˆ H A , with W = 2, the localization is not strong enough to produce such a gap: wehowever expect this to happen for larger values of the disorder amplitude W . Figure 4.
Distributions of the local density ˆ n j at the center of the chain, j = L/
2, for the diagonal ensemble and the GGE, and for L = 256 (black curves)or L = 128 (gray curves). In panel (a), we plot ln[ ρ D ( n )] /L and ln[ ρ GGE ( n )] /L ,while in panels (b) and (c) ln[ ρ D ( n )] and ln[ ρ GGE ( n )]. The vertical lines are thediagonal and GGE average of ˆ n j , which coincide for the local density. The threepanels are obtained using the same quenches of Fig. 1. Since ˆ n j is a one-body operator, the diagonal and GGE averages coincide [11],and therefore, the mean value of the two distributions is the same: (cid:90) d n ρ D ( n ) n = (cid:90) d n ρ GGE ( n ) n . ONTENTS Figure 5.
Squared single-particle wavefunction | u jµ | as a function of theeigenstates index µ , at fixed site j = L/
2. We have taken L = 256 and thethree panels are obtained using the same quenches of Fig. 1. We also note that, for this observable, ˆ n mj = ˆ n j for any positive integer m , andtherefore (cid:104) ˆ n mj (cid:105) D = (cid:104) ˆ n mj (cid:105) GGE . However, this does not allow us to conclude that the twodistributions ρ D ( n ) and ρ GGE ( n ) coincide, since, unlike the case of the total energy, wehave that, for instance: (cid:90) d n ρ D / GGE ( n ) n (cid:54) = (cid:104) ˆ n j (cid:105) D / GGE . The variance of the two distributions can be computed by exploiting again Wick’stheorem. We find that: σ n, GGE = (cid:88) µ | u jµ | n µ (cid:0) − n µ (cid:1) (22) σ n, D = σ n, GGE − (cid:88) µ (cid:54) = µ | u jµ | | u jµ | (cid:12)(cid:12) G µ µ (cid:12)(cid:12) . (23)In Fig. 6 we plot σ n, GGE and σ n, D as a function of size. We see that, in bothensembles, the variances vanish as 1 /L when quenching to ˆ H lrh with γ = 0 . finite when quenching to ˆ H A and to ˆ H lrh with γ = 2, i.e.,when the final Hamiltonian has localized eigenstates. These results agree with thefindings of Ref. [32], who show that, for large L , the variance of few-body intensive(but not site-averaged) observables remains finite both in the microcanonical ensembleand in the diagonal ensemble for the Aubry-Andr´e model. ONTENTS σ n, GGE , we see that it is related to an inverse participationratio (IPR): the sum is over the eigenstates µ , each µ weighted with the correspondingoccupation factor 0 ≤ n µ (cid:0) − n µ (cid:1) ≤ / σ n, GGE ≤ (cid:88) µ | u jµ | = IPR j , (24)where the last equality defines the IPR at fixed site j . This shows that, whenever theIPR goes to zero, i.e., when the final Hamiltonian has delocalized eigenstates, σ n, GGE goes to zero as well. For a final Hamiltonian with localized eigenstates we have insteadthe opposite: there is at least one eigenstate ˜ µ localized around j , and therefore thereis a single-particle wavefunction u j ˜ µ which does not vanish in the thermodynamiclimit; if the initial occupation n µ of this eigenstate is such that 0 < n µ <
1, then σ n, GGE remain finite in the thermodynamic limit.Concerning σ n, D , Eq. (23) can be rewritten as: σ n, D = σ n, GGE − δ jj , (25)where δ jj denotes the mean squared time-fluctuations of the single-particle Green’sfunction G j j ( t ) [10]: δ j j ≡ lim t →∞ t (cid:90) t d t (cid:48) | δG j j ( t (cid:48) ) | , (26) δG j j ( t ) ≡ G j j ( t ) − G j j being the time fluctuation with respect to the long-timeaverage G j j . Physically, δ jj is the averaged long-time fluctuation of the local densityˆ n j = ˆ c † j ˆ c j . In Refs. [10, 11] we have shown that if the final Hamiltonian has extendedeigenstates, then δ jj ≈ /L for large sizes, while δ jj remains finite when the finalHamiltonian has localized eigenstates. This explains all the features shown in Fig. 6,in particular the clear difference between σ n, D and σ n, GGE in all cases.
4. Summary and conclusions
In this paper we have introduced a Monte Carlo method — obtained by a rathernatural extension of the Wang-Landau algorithm [19, 20, 21] — which allows tocompute quite general weighted distribution functions of the form relevant to quantumquenches, see Eq. (6). We have used this approach to analyze quantum quenchesfor free-fermion Hamiltonians in presence of disorder. For these systems, thanks toWick’s theorem, after-quench expectation values and time averages require a modestcomputational effort, proportional to a power-law of the size L [11]. However, thecalculation of full probability distributions — like the diagonal ensemble distribution ρ D ( A ), Eq. (3), or the GGE one ρ GGE ( A ), Eq. (5) — would still require a sum overan exponential number of terms, hence unfeasible beyond very small sizes.Although quadratic, hence with an extensive number of conserved quantities,these free-fermion problems are not described by the GGE ensemble whenever thedisorder is such that the after-quench eigenstates are localized. More precisely, whilethe GGE ensemble is known to correctly capture the long-time average of any one-body operator, almost “by construction” [11], it does not capture correlations inducedby the spatial localization of the eigenstates. Our study further explored this issueby explicitly calculating and comparing the full probability distributions of both theenergy and the local density in the two relevant ensembles. ONTENTS Figure 6.
Plot of the variances σ n, D (circles) and σ n, GGE (triangles) as a functionof size. The data are obtained using the same set of quenches used in Fig. 2 andthe values are computed using Eq. (23). Error bars are calculated by averagingover 20 different realizations of the disorder. The dashed lines are power-law fits σ n ∼ L − s , where s ≈ Concerning the energy, we have explicitly verified that the form of the twodistributions for the diagonal and GGE ensembles differs considerably, most notablyat the extremes of the spectrum, and for the Anderson model case. More in detail,we have verified that, regardless of the final Hamiltonian, the averaged fluctuations ofthe energy-per-site, (cid:2) σ e, D (cid:3) av and (cid:2) σ e, GGE (cid:3) av , go to zero in the thermodynamic limit,see Fig. 2, in agreement with the general analysis of Refs. [12, 13]. Nevertheless, wefind that there is a clearly detectable difference in the two variances when the finalHamiltonian has localized eigenstates.In addition to the energy, we studied the local density distributions. For thisobservable, it was already known that, even in presence of disorder and localization,the GGE expectation value coincides with the diagonal average [11], a property true,more generally, for any one-body operator [11]. Our numerical results confirm thateven if the averages of ρ D ( n ) and ρ GGE ( n ) coincide, the two distributions are differentwhen localization is present, with clearly detectable differences already at the levelof the variance, see Fig. 6: σ n, GGE and σ n, D differ by a quantity which representsthe averaged long-time fluctuations of the local density [10, 11], which remain finitewhenever the final Hamiltonian has localized eigenstates.Other many-body operators, like density-density correlations, might be analyzedin a similar way. Here, even the average values are not in general well described by theGGE distribution whenever localization is at play [11]: we expect, once again, clearlyvisible discrepancies between the diagonal and GGE distributions in such cases.In conclusion, as explained above, the weighted-WLA we have presented ONTENTS
Acknowledgments
We acknowledge discussions with A. Laio, A. Russomanno and A. Silva. Research wassupported by MIUR, through PRIN-2010LLKJBX-001, by SNSF, through SINERGIAProject CRSII2 136287 1, by the EU-Japan Project LEMSUPER, and by the EU FP7under grant agreement n. 280555.
References [1] Immanuel Bloch, Jean Dalibard, and Wilhelm Zwerger. Many-body physics with ultracoldlattices.
Rev. Mod. Phys. , 80:885–964, 2008.[2] Maciej Lewenstein, Anna Sanpera, Veronica Ahufinger, Bodgan Damski, Aditi Sende, and UjjwalSen. Ultracold atoms in optical lattices: mimicking condensed matter physics and beyond.
Advances in Physics , 56:243–379, 2007.[3] J. M. Deutsch. Quantum statistical mechanics in a closed system.
Phys. Rev. A , 43:2046–2049,1991.[4] Mark Srednicki. Chaos and quantum thermalization.
Phys. Rev. E , 50:888–901, 1994.[5] M. Rigol, V. Dunjko, and M. Olshanii. Thermalization and its mechanism for generic isolatedquantum systems.
Nature , 452:854–858, 2008.[6] Marcos Rigol. Breakdown of thermalization in finite one-dimensional systems.
Phys. Rev. Lett. ,103:100403, 2009.[7] Marcos Rigol. Quantum quenches and thermalization in one-dimensional fermionic systems.
Phys. Rev. A , 80:053607, 2009.[8] Toshiya Kinoshita and Trevor Wenger and David S. Weiss. A quantum newton’s cradle.
Nature(London) , 440:900, 2006.[9] Anatoli Polkovnikov, Krishnendu Sengupta, Alessandro Silva, and Mukund Vengalattore.Nonequilibrium dynamics of closed interacting quantum systems.
Rev. Mod. Phys. , 83:863,2011.[10] Simone Ziraldo, Alessandro Silva, and Giuseppe E. Santoro. Relaxation dynamics of disorderedspin chains: Localization and the existence of a stationary state.
Phys. Rev. Lett. , 109:247205,Dec 2012.[11] Simone Ziraldo and Giuseppe E. Santoro. Relaxation and thermalization after a quantumquench: Why localization is important.
Phys. Rev. B , 87:064201, Feb 2013.[12] M. Rigol, V. Dunjko, and M. Olshanii. Thermalization and its mechanism for generic isolatedquantum systems.
Nature , 452(7189):854–858, 2008.[13] Giulio Biroli, Corinna Kollath, and Andreas M. L¨auchli. Effect of rare fluctuations on thethermalization of isolated quantum systems.
Phys. Rev. Lett. , 105:250401, Dec 2010.[14] Marcos Rigol, Alejandro Muramatsu, and Maxim Olshanii. Hard-core bosons on opticalsuperlattices: Dynamics and relaxation in the superfluid and insulating regimes.
Phys. Rev.A , 74:053616, 2006.[15] Marcos Rigol, Vanja Dunjko, Vladimir Yurovsky, and Maxim Olshanii. Relaxation in acompletely integrable many-body quantum system: An ab initio study of the dynamics ofthe highly excited states of 1d lattice hard-core bosons.
Phys. Rev. Lett. , 98:050405, 2007.[16] T. Barthel and U. Schollw¨ock. Dephasing and the steady state in quantum many-particlesystems.
Phys. Rev. Lett. , 100:100601, Mar 2008.[17] Pasquale Calabrese, Fabian H. L. Essler, and Maurizio Fagotti. Quantum quench in thetransverse-field ising chain.
Phys. Rev. Lett. , 106:227203, Jun 2011.
ONTENTS [18] Miguel A. Cazalilla, Anibal Iucci, and Ming-Chiang Chung. Thermalization and quantumcorrelations in exactly solvable models. Phys. Rev. E , 85:011133, Jan 2012.[19] Fugao Wang and D. P. Landau. Efficient, multiple-range random walk algorithm to calculatethe density of states.
Phys. Rev. Lett. , 86:2050–2053, Mar 2001.[20] Fugao Wang and D. P. Landau. Determining the density of states for classical statistical models:A random walk algorithm to produce a flat histogram.
Phys. Rev. E , 64:056101, Oct 2001.[21] D.P. Landau and F. Wang. Determining the density of states for classical statistical models bya flat-histogram random walk.
Computer Physics Communications , 147:674 – 677, 2002.[22] R. E. Belardinelli, S. Manzi, and V. D. Pereyra. Analysis of the convergence of the 1 t and wang-landau algorithms in the calculation of multidimensional integrals. Phys. Rev. E , 78:067701,Dec 2008.[23] Pedro Ojeda, Martin E. Garcia, Aurora Londo˜no, and N. Y. Chen. Monte carlo simulations ofproteins in cages: Influence of confinement on the stability of intermediate states.
BiophysicalJournal , 96(3):1076 – 1082, 2009.[24] Evelina B. Kim, Roland Faller, Qiliang Yan, Nicholas L. Abbott, and Juan J. de Pablo. Potentialof mean force between a spherical particle suspended in a nematic liquid crystal and asubstrate.
The Journal of Chemical Physics , 117(16):7781–7787, 2002.[25] M. M¨uller and J.J. de Pablo. Simulation techniques for calculating free energies. In
ComputerSimulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 1 ,volume 703 of
Lecture Notes in Physics , pages 67–126. Springer Berlin Heidelberg, 2006.[26] B. J. Schulz, K. Binder, M. M¨uller, and D. P. Landau. Avoiding boundary effects in wang-landausampling.
Phys. Rev. E , 67:067102, Jun 2003.[27] M. Gertsenshtein and V. Vasilev. Waveguides with random inhomogeneities and brownianmotion in the lobachevsky plane.
Theory of Probability & Its Applications , 4(4):391–398,1959.[28] Alexander D. Mirlin, Yan V. Fyodorov, Frank-Michael Dittes, Javier Quezada, and Thomas H.Seligman. Transition from localized to extended eigenstates in the ensemble of power-lawrandom banded matrices.
Phys. Rev. E , 54:3221–3230, Oct 1996.[29] E. Cuevas, M. Ortu˜no, V. Gasparian, and A. P´erez-Garrido. Fluctuations of the correlationdimension at metal-insulator transitions.
Phys. Rev. Lett. , 88:016401, Dec 2001.[30] Imre Varga. Fluctuation of correlation dimension and inverse participation number at theanderson transition.
Phys. Rev. B , 66:094201, Sep 2002.[31] S. Ziraldo.
Thermalization and relaxation after a quantum quench in dis-ordered Hamiltonians . PhD thesis, SISSA, Trieste, 2013. Available at: .[32] Kai He, Lea F. Santos, Tod M. Wright, and Marcos Rigol. Single-particle and many-bodyanalyses of a quasiperiodic integrable system after a quench.