Loss of ergodicity in a quantum hopping model of a dense many body system with repulsive interactions
aa r X i v : . [ c ond - m a t . d i s - nn ] S e p Loss of ergodicity in a quantum hopping model ofa dense many body system with repulsive interactions
Kazue Matsuyama
Physics and Astronomy DepartmentSan Francisco State University San Francisco, CA 94132, USA (Dated: September 15, 2020)In this work we report on a loss of ergodicity in a simple hopping model, motivated by the Hubbard Hamil-tonian, of a many body quantum system at zero temperature, quantized in Euclidean time. We show that thisquantum system may lose ergodicity at high densities on a large lattice, as a result of both Pauli exclusion andstrong Coulomb repulsion. In particular we study particle hopping susceptibilities and the tendency towards par-ticle localization. It is found that the appearance and existence of quantum phase transitions in this model, in thecase of high density and strong Coulomb repulsion, depends on the starting configuration of particle trajectoriesin the numerical simulation. This is presumably the Euclidean time version of a breakdown of the eigenstatethermalization hypothesis in real time quantization.
I. INTRODUCTION
In recent years there has been some interest in the phe-nomena of non-ergodicity and very slow thermalization inquantum many-body systems [1] associated with many-bodylocalization in random potentials [2–4], and with certaintranslation-invariant quantum hopping models [5, 6]. Therehas also been a great deal of effort devoted to the sign problemin the Hubbard model at finite density, although this is still awork in progress. Motivated largely by the sign problem in theHubbard model, I investigate here the Euclidean-time quanti-zation, in discretized time, of a many-body hopping modelincorporating features which are reminiscent of the Hubbardmodel. What I will show here is that this quantized hoppingmodel, which can also be regarded as the statistical mechanicsof interacting particle trajectories, exhibits clear non-ergodicbehavior, i.e. a strong and qualitative dependence of expec-tation values on the initial configuration, which includes theappearance, or non-appearance, of quantum phase transitionsdepending on the initial state. This seems to be another exam-ple, in addition to the cases cited above, of non-ergodicity ina translation-invariant quantum system, in which there is pre-sumably a violation of the eigenstate thermalization hypothe-sis [7–10].
II. THE MODEL
The initial motivation was to approximate some features ofthe Hubbard Hamiltonian H = − t ∑ h j , i i , σ ( c † j σ c i σ + c † i σ c j σ )+ U ∑ j n j ↑ n j ↓ − µ ∑ j ( n j ↑ + n j ↓ ) (1)in a simpler model which is hopefully more tractable to nu-merical simulation at high densities. It is a long-standing con-jecture, in the condensed matter community, that this simpleHamiltonian system describes the behavior of strongly corre-lated electrons in solids. The roadblock, at least so far as stan-dard numerical algorithms are concerned, is the sign problem.There are a number of approaches, e.g. the complex Langevinequation [11], deformation of path integration into the com- plex plane (the “thimble” approach) [12], and the density ofstates method [13], which have been investigated, largely inthe high energy physics community, in an effort to deal withthe sign problem associated with the QCD phase diagram. Al-though successes with these methods are so far limited, therehave been recent efforts to import them to deal with the signproblem in many body systems [14–17]. Whether one or moreof these approaches will be useful in the case of the Hubbardmodel (or the QCD phase diagram, for that matter) is not yetclear; what is known is that the problem of finding a general solution of the sign problem lies in the NP-hard complexityclass [18].In this article we do not attempt the direct simulation of theHubbard model. Instead we construct a simplified hoppingmodel which is free of the sign problem, but which might in-corporate at least some of the same physics. The model con-tains two types of particles, which we refer to as “spin up”and “spin down,” with a property which we will refer to, a lit-tle loosely, as the “Exclusion Principle,” meaning that no twoparticles of the same type can occupy the same lattice site.Euclidean time is also discretized, and a particle may hop inone time step to either a nearest or next-nearest lattice site if(i) this transition is allowed by the Exclusion Principle; and(ii) the spatial separation between a particle’s position at timestep t , and the position at time steps t ±
1, does not exceedthe nearest and next-nearest criterion. The Euclidean action isgiven by S = N t ∑ t = ( κ n p ∑ n = j ( n , t ) + ∑ x , y V ( x , y , t ) ) = n p ∑ n = K ( n ) + ∑ x , y , t V ( x , y , t ) (2)where n p is the total number of particles on the finite lattice, N t is the extension in the time direction, with j ( n , t ) = n is at the same lattice site attime t +
11 if particle n is at a nearest or next-nearestsite at time t + n at time t , and V ( x , y , t ) = ( x , yU if two electrons of opposite spins at site x , y (4)at time t . We have K ( n ) = κ × the number of hops along thetrajectory of the n -th particle, and we count it as one “hop”whenever j ( n , t ) =
1. The nearest and next-nearest constraint,and the Exclusion Principle constraint, are understood.The quantized Hamiltonian of this system would followfrom the transfer matrix in the continuous time limit, althoughwe will not be concerned with that limit here. Periodic bound-ary conditions in the space and time dimensions are imposed,and we are able to vary the inverse temperature in lattice unitsby varying the extension in the time direction. For the mostpart we use N t = n p , and for an L × L lat-tice n p = L corresponds to half-filling, as in a true fermionicsystem.For numerical simulation of the model via importance sam-pling, employing the usual Metropolis algorithm, we calculatethe change in this action ∆ S = κ∆ j + ∆ V (5)resulting from a trial update in the trajectory of one of the par-ticles in the system. We will be interested in studying thebehavior of the system at high density (i.e. half-filling andabove) as the U and κ parameters are varied.It should be obvious that despite some similarities, themodel of eq. (2) is not the Hubbard model. It is really ahopping model describing a dense set of two types of distin-guishable particles (“spin up” and “spin down”) , with the con-straint that no more than one particle of either type can occupyany given site of the lattice. There is at least one historicalprecedent for treating fermions in that way, namely the origi-nal Hartree formulation of the Hartree-Fock approximation inatomic physics. The Hartree-Fock approximation, of course,consists of writing down the ground state of a set of electronsmoving in a central potential determined by self-consistency.In the original Hartree version, the minimal energy state of N electrons is chosen from many body states of the form Ψ ( , , ..., N ) = ψ s ( ) ψ s ( ) .... ψ s N ( N ) (6)where ψ s ( j ) indicates that particle j is in the energy eigen-state ψ s of the one-body Schrodinger equation for an electronmoving in a central potential. The restriction was the PauliExclusion Principle: no two particles could be in quantumstates labeled by the same set of integers s . Of course it was soon noted that this application of the Exclusion Principle isinsufficient, and that a many-body wave function of this kindis appropriate to distinguishable particles rather than identicalfermions. The remedy was to impose antisymmetrization viathe Slater determinant, and this improved method is knownas “Hartree-Fock.” However, it is worth noting that the origi-nal Hartree approximation was not so terrible at the quantita-tive level. Agreement with experiment was certainly improvedwith the appropriate antisymmetrization, but the Hartree ver-sion already gives reasonable results for atomic structure, dif-fering from the more accurate Hartree-Fock method, in esti-mates of atomic energy levels, at the 10-20% level [19]. Thehopping model I have just introduced imposes an exclusionprinciple on double occupancy, rather than on energy eigen-states. But one might hope that if the particles are reasonablylocalized at the quantum level, then there might still be someresemblance to the physics of the Hubbard model, despite theclear violation of Fermi-Dirac statistics. In any case, giventhe absence of any robust computational solution of the Hub-bard model away from half-filling, we believe that the inves-tigation of this simplified (and, as regards identical particlestatistics, evidently wrong) version of that model may still beworth pursuing. Of particular interest would be the occurrenceof quantum phase transitions in this Euclidean-time quantizedhopping model.
A. Observables of the hopping model
We numerically simulate the system we have described viathe standard Metropolis algorithm. The procedure is to gotime slice by time slice, at each time updating the location ofeach particle on the two-dimensional lattice. In the Metropolisupdate for each particle we choose at random a trial hoppingdirection (to a nearest or next-nearest neighbor site), and mea-sure the change in the number of hops ∆ j along the trajectoryof the given particle (this must be in the range − ≤ ∆ j ≤ ∆ S in(5), and the change is then accepted or rejected in the usualway. The trial move is constrained by the Exclusion Princi-ple, and the restriction that a change in particle position attime t should not result in a hop from time t − t , or time t to t +
1, such that the distance between sites at the earlierand later times exceeds the nearest or next-nearest neighborhopping limit.The two dimensional lattice size of the simulations has beenvaried from 10 ×
10 to 50 ×
50 sites in order to study the vol-ume dependance of the observables. Because there are equalnumbers of each type of particle, the maximum number ofparticles allowed on an L × L lattice is 2 L .Note that the Euclidean path integral represents the classi-cal statistical mechanics of a discretized system of trajectories There are of course instances where the quantitative discrepancy betweenthe Hartree and Hartree-Fock methods is more severe, e.g. when the ex-change interaction is crucial, cf. [20]. (which might be thought of as “fibers” of some kind), withrepulsive interactions along with the Exclusion Principle justmentioned. If there is non-ergodicity in the quantum systemof point particles, it would be associated with non-ergodicityin this classical system of line-like objects. In other words, theergodicity of our simplified hopping model corresponds to theindependence of expectation values on the initial configura-tion of particle trajectories in the Euclidean time formulation.We think this is essentially equivalent to a test of the eigen-state thermalization hypothesis in the corresponding real timeevolution of a quantum state.In our computation, we calculated (i) the hop susceptibility;and (ii) the probability that a particle remains in the initialsite after a Euclidean time lapse t . These are computed as afunction of κ for various U values, inverse temperature (theextension of the lattice in the inverse time direction), and fordifferent densities.Let us define the average hopping number at time t hop ( t ) = n p n p ∑ n = ( j ( n , t − ) + j ( n , t )) (7)and the corresponding hopping susceptibility χ hop = n p N t N t ∑ t = ( h hop ( t ) i − h hop ( t ) i ) (8)This observable has been defined so as to be local in time, andin a real time quantization could presumably be computed, atzero temperature, from the ground state wave functional.In order to investigate localization, we define n ( t , t + t ) = no. of particles which are at the samelattice site at times t and t + t (9)and from this quantity we compute the probability that a par-ticle will remain in the same position for t time steps P ( t ) = * N t N t ∑ t = n ( t , t + t ) n p + (10) III. RESULTSA. Non-ergodicity at half-filling
For half-filling (particle density 50% of maximum), wecomputed observables with two different starting configura-tions, which we term “random” and “minimum energy” (orjust “minimum”) respectively. Initially all particle trajectoriesare constant in time, but in the random configuration the x , y position of each trajectiory is chosen at random, apart fromthe constraint of the Exclusion Principle. So with this initial-ization a certain fraction of sites are doubly occupied at eachtime. In the minimum energy configuration there is one parti-cle per site, and no doubly occupied sites, with up and down spins alternating in an antiferromagnetic pattern. In these ini-tial configurations we have K ( n ) = U = χ hop must obviously be zero for all κ values, at such alarge value of U , and of course that is what one finds (Fig. 1). χ hop κ U = 100, N t = 100, minimun energy configuration start L 10 L 20 L 30 L 40 L 50 FIG. 1. A trivial result: hop susceptibility ( χ hop of eq. (8)) vs. κ atlow temperature ( N t = 100), strong repulsion U = L , for 50% density, initialized at the minimum energyconfiguration. χ hop κ U = 100, N t = 100, random configuration start L 10 L 20 L 30 L 40 L 50 FIG. 2. χ hop vs. κ at low temperature ( N t = 100), strong repulsion( U = L , 50% density with a randominitial configuration. On the other hand, in the random configuration start, the sit-uation is dramatically different. Data was taken on L × L lat-tices for L =
10 up to L =
50, and N t = χ hop at U =
100 obtained ina numerical simulation with 1000 thermalizing sweeps, fol-lowed by 9000 sweeps, with data taken every 100 sweeps.One might wonder if the system simply needs more MonteCarlo time to find its way to the minimum energy configura-tion. In Figure 3 the number of thermalizing and data takingsweeps have been increased by an order of magnitude, withvery little change in the hop susceptibility. Presumably, at asufficiently large number of sweeps on a finite lattice, even-tually the minimal energy configuration would be obtained,if there is any allowable path to that configuration. We thinkit is likely however, given these results, that the number ofsweeps required to reach a minimum energy configurationwould rapidly diverge with volume, while keeping the densityat half-filling fixed. χ hop κ U = 100, N t = 100, random configuration, longer simulation L 10 L 20 L 30 L 40 L 50 FIG. 3. Same as Fig. 2 but with both thermalizing and subsequentMonte Carlo sweeps increased by a factor of 10, to at total of 10 sweeps. This seems to make little difference to the final results. At a lower value of U =
10 there can be particle hoppingeven with a minimum energy initialization, as seen in Fig.4(a). Since the height of the susceptibility peak is volumedependent, there is an apparent suggestion here of a quan-tum phase transition. By contrast, a random initialization (Fig.4(b)) shows quite a different pattern, with no indication what-ever of a phase transition.In addition to the random and minimum initial configura-tions, we have also investigated what could be described as“annealing” initializations. For both random and minimumconfigurations, we initialize to a random or minimum con-figuration at each κ value, before evolving the system via theMetropolis algorithm. Instead, for the annealing initialization,the system is only initialized, with a random initial configu-ration, at κ =
0. After computing the hop susceptibility at κ =
0, we use the configuration obtained at the last MonteCarlo sweep as the initial configuration at the next value of κ , i.e. κ = δκ , with δκ = .
02. The last configuration of thesimulation at that κ value is then used as the initial configura-tion at the next κ value at κ = δκ , and so on, with the lastconfiguration at κ = n δκ used as the initial configuration forthe κ = ( n + ) δκ simulation, over the full range of κ usedin the calculation.The data with annealing initialization at U =
10 and 50%density produced the surprising results shown in Fig. 4(c).The reason this is surprising is that we see a substantial hop-ping susceptibility even at rather large κ , which should havesuppressed all hopping. But this is clearly an effect of the ini-tialization. At small values of κ it is not surprising that there χ hop κ U = 10, N t = 100, minimum energy configuration L 10 L 20 L 30 L 40 L 50 (a) χ hop κ U = 10, N t = 100, random configuration start L 10 L 20 L 30 L 40 L 50 (b) χ hop κ U = 10, N t = 100, annealing start L 10 L 20 L 30 L 40 L 50 (c) FIG. 4. χ hop vs κ at weaker repulsion U = , N t = L . Initialization in subfigures at: (a) minimumenergy; (b) random; (c) annealing start configurations. will be some non-negligible number of jumps along any giventrajectory. But since the last configuration at κ = n δκ is theinitial configuration of the next simulation at κ = ( n + ) δκ ,there is the possibility that the multiple jumps found in tra-jectories at low κ become “frozen in” as κ is raised to highervalues, and this appears to be what happens. At large κ theenergetics would prefer a low K ( n ) value, with the minimum K ( n ) = κ simply indicates that the system of trajectories cannot find itsway, at large κ , to anywhere near that minimum. We will re-turn to this point in section III D.These results at both U =
100 and U =
10, which are sovery clearly dependent on the starting configuration, are a firstindication of non-ergodicity (or exceptionally long relaxationtimes) in the classical system of particle trajectories, and a cor-responding non-ergodicity in the associated quantum systemof point particles.
B. Non-ergodicity at 60% filling
1. U = 100
We now increase the particle density to 60% of maximum,and compute the hop susceptibilities for random, minimum,and annealing initializations. At 60% we have investigatedtwo types of minimum energy initializations. After placingthe up and down spins alternately on each site, so that thereis one particle at each lattice site, we can place the remainingparticles at random (with the constraint of the Exclusion Prin-ciple), to form doubly occupied sites distributed randomly inthe lattice. We will refer to such an initialization as “minimumrandom”. An alternative is to add remaining particles in onecorner of the lattice, so that the doubly occupied sites are lo-cated in one connected region of the lattice. We will refer tothis initialization as simply “minimum.”It is interesting that at 60% filling, with random, minimum,and minimum random initializations, we seem to see (judg-ing from the hop susceptibility) evidence of a quantum phasetransition. Yet the plots of susceptibilities differ for the threeinitializations. In Fig. 5(a), the random initial particle con-figuration shows there are two peaks which both grow withspatial volume, indicating two quantum phase transitions, butwith both minimum initialization Fig. 5(b)) and minimum ran-dom initialization (Fig. 5(c)) there is only a single peak. Wedoubt that the second peak in Fig. 5(a) is physical, for reasonsdiscussed in section III D below.We have nonetheless computed the critical exponent fromfinite size scaling for the two peaks seen in Fig. 5(a). Themaximum peak height vs. volume is displayed on a log-logplot in Fig. 6, and a best fit of ( χ hop ) max to aL ( γ / ν ) , where a isa constant, gives ( γ / ν ) first = 1.37 ± a = 0 . ± . ( γ / ν ) second = 1.45 ± a = 0 . ± .
001 for the second peak.As before, we check if our calculation depends on the num-ber of Monte Carlo sweeps. The corresponding results ob-tained by increasing both thermalization and subsequent up-date sweeps by a factor of 10 are shown in Fig. 7. The peaksremain in hop susceptibilities and their positions have shiftedslightly. But while the peak heights for the random and mini-mum initializations have decreased, the peak height for mini-mum random is about the same after increasing the number ofsweeps by an order of magnitude. χ hop κ U = 100, N t = 100, random initialization L 10 L 20 L 30 L 40 L 50 (a) χ hop κ U = 100, N t = 100, minimum initialization L 10 L 20 L 30 L 40 L 50 (b) χ hop κ U = 100, N t = 100, minimum random initialization L 10 L 20 L 30 L 40 L 50 (c) FIG. 5. χ hop vs κ at 60% filling. Strong repulsion U = N t = L . Initialization in the subfigures isat (a) random; (b) mimimum; (c) minimum random configurations. The most dramatic difference is seen with the annealing ini-tializations, with hop susceptibilities displayed in Fig. 8. Withthis initialization there is essentially no dependence of peaksize on lattice area, and no evidence whatever of a quantumphase transition. But we observe that the data on a 10 × l og ( s u s pea k ) log(L)Critical Exponent N t = 100, U = 100, random initialization1st peak2nd peak1st peak fit2nd peak fit FIG. 6. critical exponent fittings for low temperature (T = 100) hop-ping susceptibility at U = 100, 60% filling random initial configura-tion laxation times, is a phenomenon which increases very rapidlywith lattice size. Presumably this translates, on a very largelattice, to a complete breakdown of eigenstate thermalizationin the corresponding real-time quantum system of point par-ticles on the lattice. In fact the situation is not so differentfrom what we have already seen at half-filling. Given a suf-ficiently lengthy simulation, the system must eventually findthe minimal energy configuration. But for a sizable lattice,the required simulation time to reach that minimal energy isprobably beyond the reach of any realizable computation.Our second observable is the localization probability P ( t ) ,which represents the probability that a particle at some time t , at a location ( x , y ) , will be found at that same site af-ter t units of Euclidean time, i.e. at time t + t . The results,for the mimimum, minimum random, and annealing initial-izations are shown in Fig. 9 for a 50 ×
50 spatial lattice. A P ( t ) which reaches a plateau at large t is indicative of stronglocalization, and this is what is seen at the larger κ values inFigs. 9(a) and 9(b). The sudden onset of localization with in-creasing κ is especially apparent in the minimum random ini-tialization (Fig. 9(b)). Note that the positions of the first peakin each plot in the hop susceptibilities in Fig. 5(b) and Fig.5(c) approximately coincide with the positions of the onset ofthe particle localization. In contrast, with the annealing ini-tialization (Fig. 9(c)) we see P ( t ) → t evenat large κ , but this is not hard to understand. In the annealedcase (here, and in all calculations of P ( t ) we use δκ = . κ are initialized to disorderedconfigurations that are inherited from the simulations at small κ . Thus there appears to be no localization, even if (as mustbe the case at large κ ) particle hopping occurs only rarely inthe course of the simulation.The plots shown in Fig. 10 are for the random initialization,where the onset of localization with κ is not quite as abrubt.Fig. 10(a) shows the numerical results for 1000 thermaliza-tions followed by 9000 sweeps, as usual, with data taken every100 sweeps. Fig. 10(b) is the same computation, but with ther-malization and subsequent Monte Carlo sweeps increased bya factor of 10; this increase seems to make little difference to χ hop κ U = 100, N t = 100, random initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (a) χ hop κ U = 100, N t = 100, minimum initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (b) χ hop κ U = 100, N t = 100, minimum random initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (c)
FIG. 7. Same as Figure 5, but with both thermalizing and subse-quent Monte Carlo sweeps increased by a factor of 10, to at total of10 sweeps. Initialization at (a) random (b) minimum; (c) minimumrandom starting configurations. the result. Again we find that the onset of localization with κ seems to coincide with the start of the first peak in hop suscep-tibility. Overall, it appears that the quantum phase transitionsseen in Figs. 5 and 7 (at least the first peak, in the case of ran-dom initialization) are associated with a localization transitionof some kind, reminscent of a glass transition.The localization data for the annealing initial configuration χ hop κ U = 100, N t = 100, annealing start L 10 L 20 L 30 L 40 L 50 FIG. 8. Low temperature ( N t = 100) hopping susceptibility as a func-tion of κ , U = is shown in Fig. 9(c) where we see that P ( t ) → t interval, indicating a lack of localization, and thereforeno abrupt transition to localization. This is consistent with ourhop susceptibilities calculations in Fig. 8 where we found thatno phase transition occurred.
2. U = , We have repeated the previous computations for U = , , U =
10 arequalitatively quite similar to the previous U =
100 case.The susceptibility and localization data for U = U = ×
10 lattice is not farfrom the result for the annealed lattice. On the other hand, wefind that the random, minimum random, and minimum ini-tializations give roughly consistent results (in particular thedouble peak structure is gone), and it is only annealing initial-ization that shows no sign of a quantum phase transition.We also repeated the susceptibility calculations by increas-ing the number of thermalization and data taking sweeps byan order of magnitude, i.e. to 10 total sweeps, in Fig. 13.The random and minimum initial configurations lead to sim-ilar results, seen in Fig. 13(a) and 13(b). Yet there is stillnon-ergodicity, because the susceptibility plot correspondingto the minimum random initialization still shows a strong sin-gle peak, as high as the one seen in Fig. 11(c), which was ob-tained with an order-of-magniture fewer Monte Carlo sweeps.In the case of no Coulomb repulsion whatever, i.e. U = P(t) plot U = 100, density 60% minimum start 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (a) minimum P(t) plot U = 100, density 60% minimum random start 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (b) minimum random P(t) plot U = 100, density 60% annealing start 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (c) annealing FIG. 9. Localization probability P ( t ) vs. t and κ on a 50 ×
50 spatiallattice at 60% filling, at U = , N t =
100 and spatial areas L . Ini-tialization at (a) minimum; (b) minimum random; (c) annealing startconfigurations. ever, if we again increase the Monte Carlo sweeps by a factorof 10, to 10 total sweeps, the situation is different, as shownin Fig. 15. The data for the larger volumes seem to converge,for this larger number of sweeps, to the data which was foundon the 10 ×
10 area, which itself fits with the simulated an-nealing data. It is true the larger lattice areas still have somepeculiar spikes; we might guess that these will disappear after
P(t) plot for U = 100, density 60% random start 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (a) P(t) plot U = 100 density 60% random start, longer simulation 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (b) FIG. 10. Same as Fig. 9 for random initializations. (a) 10 total Monte Carlo sweeps; (b) 10 total Monte Carlo sweeps. χ hop κ U = 1, N t = 100, random initialization L 10 L 20 L 30 L 40 L 50 (a) random χ hop κ U = 1, N t = 100, minimum initialization L 10 L 20 L 30 L 40 L 50 (b) minimum χ hop κ U = 1, N t = 100, minimum random initialization L 10 L 20 L 30 L 40 L 50 (c) minimum random χ hop κ U = 1, N t = 100, annealing start L 10 L 20 L 30 L 40 L 50 (d) annealing FIG. 11. χ hop vs. κ at weaker ( U =
1) repulsion and 60% filling, N t = still longer simulations. If so, the comparison of U = U cases suggests that relaxation times, assuming theyare finite at higher U , increase both with lattice size and withCoulomb repulsion. Again the analogy to a glass transition issuggestive. C. Ergodicity at a higher temperature: N t = , U = , filling The next question is whether ergodicity is restored (or re-laxation times reduced) at high temperature. In this connec-tion we decrease the time extension from N t =
100 to N t = U =
100 and 60% filling. The results
P(t) plot U = 1, N t = 100, random start 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (a) random P(t) plot U = 1, N t = 100, minimum start 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (b) minimum P(t) plot U = 1, N t = 100, minimum random start 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (c) minimum random P(t) plot U = 1, N t = 100, annealing start 5 10 15 20 25 30 35 40 45 50t 0 1 2 3 4 5 6 7 8 κ (d) annealing FIG. 12. Localization probability plots at low temperature ( N t = 100) as a function of κ and t on a 50 ×
50 lattice, U = χ hop κ U = 1, N t = 100, random initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (a) random χ hop κ U = 1, N t = 100, minimum initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (b) minimum χ hop κ U = 1, N t = 100, minimum random initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (c) minimum random
FIG. 13. Same as Figure 11(a-c), with Monte Carlo sweeps increased by a factor of 10, to 10 total. Initialization: (a) random; (b) minimum;(c) minimum random starting configurations. are shown in Fig. 16. All of susceptibility plotting results aredifferent from each other, indicating as before lack of ergod-icity, at least up to 10000 total sweeps. The one difference weobserve here is in random initial configuration result in Fig.16(a), there the first peak existing in the low temperature plotin Fig. 5(a) disappeared in high temperature plot. The singlepeaks arising in the minimum and minimum random initial-izations at low temperature remained in the high temperatureplots. The positions of those single peaks at the higher andlower temperatures are found about the same κ values. Withthe annealing initialization, the volume dependent peaks donot exist in Fig. 16(d).Now we again increase both the thermalization and data taking sweeps by a factor of 10, to 10 total update sweeps,and repeat our calculations in Fig. 17. We find that this longersimulation time eliminates the volume dependent single peaksin the hopping susceptibility plot in minimum and minimumrandom configurations in Fig. 17(b) and Fig. 17(c), and bringsboth of these plots close to the annealed result in Fig. 16(d). Yet even in this case, the data for random initialization hasnot converged to the annealed start; in that the second peak There is still a spike in the minimum random data, but its height is greatlyreduced compared to Fig. 16(c), but we may guess that this spike will dis-appear entirely in a still longer simulation. χ hop κ U = 0, N t = 100, random initialization L 10 L 20 L 30 L 40 L 50 (a) random χ hop κ U = 0, N t = 100, minimum initialization L 10 L 20 L 30 L 40 L 50 (b) minimum χ hop κ U = 0, N t = 100, minimum random initialization L 10 L 20 L 30 L 40 L 50 (c) minimum random χ hop κ U = 0, N t = 100, annealing initialization L 10 L 20 L 30 L 40 L 50 (d) annealing FIG. 14. χ hop vs. κ at zero repulsion ( U = N t = remains. Its position has shifted to a larger κ value, as com-pared to the fewer sweeps result in Fig. 16(a), but the heightof the second peak has stayed about the same.So the data regarding ergodicity at this much higher tem-perature are not unambiguous. The annealed, minimum, andminimum random initializations tend, eventually, to the sameresult. But the random initialization, with the peculiar secondpeak, still differs, and the question is whether this differenceis significant. D. The second peak at random initialization
To study the nature of the second peak, we have computedthe average hopping probability for nearest neighbor and next-nearest neighbor (diagonal) separately, with hopping contribu-tions denoted j + and j × respectively, corresponding to a hopat time t to a nearest or next nearest site at time t +
1. We workagain at U =
100 and 60% filling at the higher and lower tem-peratures of N t = N t = ( t ) = n p n p ∑ n = ( j × ( n , t − ) + j × ( n , t )) (11)NN hop ( t ) = n p n p ∑ n = ( j + ( n , t − ) + j + ( n , t )) (12) These quantities are then averaged over all time slices, andgive us the diagonal and nearest neighbor hopping proba-bilities, per particle per time, respectively. The results for10,000 total Monte Carlo sweeps are shown in Fig. 18. Near-est and next nearest neighbor hopping probabilities are aboutthe same up to the onset of the second peak, after which thenext nearest hopping seems preferred, both at low and hightemperatures.What is the explanation of this strange second peak whichappears with the random initialization, and of the hoppingprobability which remains surprising high at large values of κ ? And why is this peak, and this surprisingly high hoppingprobability at large κ not also seen for the minimum and min-imum random initializations? We think the reason is actuallystraightforward. With a random initialization, in contrast tothe minimum and minimum random configurations, the sys-tem begins the simulation with some number of unoccupiedsites. At U =
100 there is a very strong tendency for a particlein a double-occupied site to move to an empty site; the costin K ( n ) is more than made up for by the drop in the very highpotential energy. Of course this introduces some degree ofhopping along each trajectory, even at rather large κ . At somepoint there are no more unoccupied sites, or at least none thatcan be readily accessed by a single hop from double occupiedsites, and the trajectories are at that stage “frozen,” becauseany further hopping, even from a double occupied site to a sin-gle occupied site (which costs nothing in potential energy) isstrongly disfavored by the kinetic part of the action. By com-1parison, in the minimum and minimum random configurationsthere are no unoccupied sites in the initial configuration, andall hopping is strongly disfavored at large κ values.The conclusion is that the second peak is rather unphysical.It originates from particles moving towards the minimal en-ergy configuration in the early part of the simulation, and thehopping which is essentially frozen into particle trajectories isalmost entirely due to the motion at that period. E. Lower density
We would expect that the non-ergodicity that we have ob-served at 50% and 60% filling is a high density phenom-ena, which one would not see at much lower densities. Asa check, we have repeated our calculations at 30% filling at U = , N t =
100 for several different initializations. In thissituation, the initializations we have described as “minimum”and “minimum random” do not apply, because there are sim-ply not enough particles to have at least one particle at everysite. We introduce instead a new initial configuration wherewe start with sites which are either unoccupied or double oc-cupied, with the double occupied sites chosen randomly onthe lattice. We refer to this initialization as “random double.”As before, there is no variation in time in the initial randomdouble configuration. The results are shown in Fig. 19. Ofcourse there are the striking double peaks at large κ which,as just explained, are a phenomenon from the initial periodof the simulation in the random and random double starts.Looking aside from the second peak, at lower κ , the random,random double, and annealed initializations are quite similar,with the relatively small volume dependence seen in Fig. 19(a)no longer visible when the Monte Carlo sweeps are increasedby a factor of 10, in Fig. 19(d). There are no indications of a quantum phase transition. We see, as expected, ergodicityat lower density and moderate values of κ , while fluctuationsare essentially frozen (this is the origin of the second peak) atlarge κ . IV. CONCLUSIONS
The correspondence between the quantum mechanics ofpoint-like particles, and the statistical mechanics of line-like particle trajectories, suggests a possible source of non-ergodicity, if a dense system of line-like trajectories begins toexhibit “glassy” behavior, where the system becomes stuck ina rather localized region in the space of relevant configura-tions. In this article we have investigated the possible loss ofergodicity at high densities in a simple hopping model, in-spired by certain features of the Hubbard model, in which“spin up” and “spin down” particles can hop in Euclidean timeon a two dimensional lattice. We have seen that the behav-ior of this model depends very strongly on the initial startingconfiguration of the Monte Carlo simulation, and in fact wehave seen apparent quantum phase transitions for some start-ing configurations, but not for others. This dependence oninitialization occurs at high densities, and we have illustratedthe situation at 50% and 60% of maximum filling. We see thatergodicity is not an issue at lower densities. Whether this non-ergodic behavior at high densities might be relevant in someexperimental situations, e.g. systems of ultra cold atoms, re-mains to be seen.
ACKNOWLEDGMENTS
I thank Jeff Greensite for helpful discussions. [1] J. P. Garrahan, Physica A: Statistical Mechanics and its Ap-plications , 130 (2018), Lecture Notes of the 14th Interna-tional Summer School on Fundamental Problems in StatisticalPhysics.[2] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod.Phys. , 021001 (2019).[3] E. Altman and R. Vosk, Annu. Rev. Condens. Matter Phys. ,383 (2015).[4] R. Nandkishore and D. A. Huse, Annu. Rev. Condens. MatterPhys. , 15 (2015).[5] Z. Lan, M. van Horssen, S. Powell, and J. P. Garrahan, Phys.Rev. Lett. , 040603 (2018).[6] M. van Horssen, E. Levi, and J. P. Garrahan, Phys. Rev. B ,100305 (2015).[7] J. M. Deutsch, Phys. Rev. A , 2046 (1991).[8] M. Srednicki, Phys. Rev. E , 888 (1994).[9] H. Tasaki, Phys. Rev. Lett. , 1373 (1998).[10] M. Rigol, V. Dunjko, and M. Olshanii, Nature , 854 (2008).[11] G. Aarts, E. Seiler, and I.-O. Stamatescu, Phys. Rev. D ,054508 (2010), arXiv:0912.3360.[12] M. Cristoforetti, F. Di Renzo, A. Mukherjee, and L. Scorzato,Phys. Rev. D , 051501 (2013), arXiv:1303.7204. [13] K. Langfeld, B. Lucini, and A. Rago, Phys. Rev. Lett. ,111601 (2012), arXiv:1204.3243.[14] C. E. Berger et al., (2019), arXiv:1907.10183.[15] M. Korner, K. Langfeld, D. Smith, and L. von Smekal, (2020),arXiv:2006.04607.[16] M. Fukuma, N. Matsumoto, and N. Umeda, Phys. Rev. D ,114510 (2019), arXiv:1906.04243.[17] A. Mukherjee and M. Cristoforetti, Phys. Rev. B , 035134(2014), arXiv:1403.5680.[18] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. , 170201 (2005),arXiv:cond-mat/0408370.[19] S. Gasiorowicz, Quantum Physics, 3rd ed (Wiley, 2007).[20] D. Pfannkuche, V. Gudmundsson, and P. A. Maksym, Phys.Rev. B , 2244 (1993). χ hop κ U = 0, N t = 100, random initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (a) random initial configuration, 10 times longerthermalization time χ hop κ U = 0, N t = 100, minimun initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (b) minimum initial configuration χ hop κ U = 0, N t = 100, minimun random initialization, longer simulation L 10 L 20 L 30 L 40 L 50 (c) minimum random initial configuration
FIG. 15. 10 times longer thermalization time. Hopping susceptibility plots of low temperature ( N t = 100) as a function of κ , U = 0 of variousvolume, 60% filling. (a) random initial configuration. (b) minimum initial configuration. (c) minimum random initial configuration. χ hop κ U = 100, N t = 5, random initialization L 10 L 20 L 30 L 40 L 50 (a) random initial configuration at high temperature χ hop κ U = 100, N t = 5, minimum initialization L 10 L 20 L 30 L 40 L 50 (b) minimum initial configuration at high temperature χ hop κ U = 100, N t = 5, minimun random initialization L 10 L 20 L 30 L 40 L 50 (c) minimum random initial configuration at hightemperature χ hop κ U = 100, N t = 5, annealing initialization L 10 L 20 L 30 L 40 L 50 (d) annealing at high temperature FIG. 16. Effect of increasing temperature by a factor of 20, to N t =
5. As in Figure 5 we set the density at 60% with strong repulsion ( U = χ hop κ U = 100, N t = 5, random initialization long simulation L 10 L 20 L 30 L 40 L 50 (a) random χ hop κ U = 100, N t = 5, minimum initialization L 10 L 20 L 30 L 40 L 50 (b) minimum initial configuration χ hop κ U = 100, N t = 5, minimun random initialization L 10 L 20 L 30 L 40 L 50 (c) minimum random initial configuration FIG. 17. Same as Figure 16 at high temperature, with Monte Carlo sweeps increased by an order of magnitude, to 10 total sweeps. Initializa-tions: (a) random; (b) minimum; (c) minimum random initial configuration. H opp i ng p r obab ili t y κ U = 100, N t = 100, random initial configuration, L = 50 Diagonal Nearest (a) low temperature H opp i ng p r obab ili t y κ U = 100, N t = 5, random initial configuration, L = 50 Diagonal Nearest (b) high temperature FIG. 18. Diagonal and nearest-neighbor hopping probabilities, as defined in eqs. (11) and (12), at U = N t = 100 (b) high temperature, N t = 5. χ hop κ U = 100, N t = 100, random initialization L 10 L 20 L 30 L 40 L 50 (a) random χ hop κ U = 100, N t = 100, random double initialization L 10 L 20 L 30 L 40 L 50 (b) random double χ hop κ U = 100, N t = 100, annealing initialization L 10 L 20 L 30 L 40 L 50 (c) annealing χ hop κ U = 100, N t = 100, random initialization, long simulation L 10 L 20 L 30 L 40 L 50 (d) random, 10 MC sweeps
FIG. 19. Lower density, 30% filling. χ hop vs. κ at strong repulsion ( U = N t ==