Ergodicity in Strongly Correlated Systems

We present a concise, but systematic, review of the ergodicity issue in strongly correlated systems. After giving a brief historical overview, we analyze the issue within the Green's function formalism by means of the equations of motion approach. By means of this analysis, we are able to individuate the primary source of non-ergodic dynamics for a generic operator and also to give a recipe to compute unknown quantities characterizing such a behavior within the Composite Operator Method. Finally, we present examples of non-trivial strongly correlated systems where it is possible to find a non-ergodic behavior.


Historical Overview
The issue of ergodicity in condensed matter physics is well known since fifties [ 1]. Given two operators A and B, describing physical quantities (e.g., charge, spin, pair densities or currents), one can study the physical response of a system described by a certain Hamiltonian H through the generalized susceptibility where h is an external field entering the Hamiltonian of the system under study in a coupling term of the type −hB and · · · stands for the statistical average in some ensemble over the perturbed system. Kubo [ 1] immediately noticed that the static isolated susceptibility χ I (0), defined for an isolated system perturbed by an external field turned on adiabatically, and the isothermal susceptibility χ T , defined for a system in thermal equilibrium in presence of a time-independent external field, can be generally different. In particular, Falk [ 2] has shown that the static isolated susceptibility is just a lower bound for the isothermal one χ I (0) ≤ χ T .
We recall that the static isolated susceptibility χ I (0), within the linear response theory, is defined through the related retarded Green's function [ 1] where · · · 0 stands for the statistical average in the microcanonical ensemble (i.e., fixing energy) on the unperturbed system and F for the Fourier transform.
On the other hand, the isothermal susceptibility χ T can be computed as In fact, starting from the expression of the thermal average in the canonical ensemble (i.e., fixing temperature) A and by the ratio where Z 0 denotes the partition function of the unperturbed system. Taking the derivative after (1.1) and exploiting the cyclic property of the trace, we obtain the isothermal susceptibility as in (1.3). Now, if we rewrite both expressions by means of the general formulas for the retarded Green's functions and the correlation functions given in the companion article [ 3] (see Section 3), present in this same issue, we get χ I (0) = 1 V k,l σ (l,−1) (k) ω l (k) (1.9) and where [ 3,4] V is proportional to the volume of the system, the sum over l ranges over the number of fields in the chosen basis, σ (l) are the spectral density functions, ω (l) are the poles of the propagator, Γ AB is an unknown function appearing in case of poles with zero value. We can immediately see that the two susceptibilities differ for the following expression Now, one can check that rewriting the expression lim t→∞ AB(t) by means of the general formula for the correlation functions given in the companion article [ 3] (see Section 3), present in this same issue, we just get and, accordingly, the difference at the r.h.s. of (1.11) is just what enters Khintchin's theorem [ 5]: a dynamics is ergodic (i.e., phase space equilibrium averages are equal to ensemble microcanical averages, which are much easier to compute) if an only if lim t→∞ AB(t) = A B (1.14) In other words a dynamic is ergodic if correlations attenuate in time. In particular, for B ≡ A, the dynamics of A is ergodic if, during its time evolution, it has non-zero matrix elements only between states within a zero-volume region of the phase space of the system [ 6]. It is clear now the link between ergodicity and response theory: the two definitions of susceptibility differ when the dynamics of the system is not ergodic. Two little, but important, notices: finite systems are not ergodic by definition, just because of the inequivalence of the ensembles; non-ergodicity at zero temperature is just the result of a degeneracy in the ground state. Several years later it was shown [ 7,8] that the difference between the two definitions of susceptibility is related to the zero-frequency anomaly exhibited by bosonic correlation functions: the presence of undetermined constants in the bosonic correlation functions. This is exactly what the relations derived above and the results of the companion article [ 3] predict establishing a definite link between the ergodicity of the dynamics and the Green's function formalism. It was first put in evidence in [ 9] and then studied by many other authors [ 10,11,7,8,6,12,13,14,15,16]. There is a general believe that this problem is of academic interest and in the last years no much attention has been dedicated to it. The main reason is that the response functions, the experimentally observed quantities, are given by retarded bosonic Green's function which formally do not depend on the such undetermined constants, which are, therefore, considered of no physical interest. The general attitude [ 1,10] is to believe that in macroscopic real systems at equilibrium at a temperature T , the fluctuations are very small and the interaction between the system and the reservoir would introduce an irreversible relaxation and decouple the correlation functions. Then, as suggested in [ 10], these constants should be always determined by requiring the ergodicity. This procedure is some how an artifice and may lead to serious problems because it might break the internal self-consistency of the entire formulation. As remarked in [ 10], the zero-frequency anomaly is a manifestation of the difficulty in extracting irreversible behavior from the statistical mechanics. This is true, but as long as we use the scheme of statistical mechanics we must be careful in doing self-consistent calculations. Breaking the self-consistency might bring to serious errors.
According to the well-known relations existing between casual (C), retarded (R) Green's functions and correlation functions the zero-frequency excitations do not contribute explicitly to the imaginary part of the retarded Green's functions and, consequently, Γ does not explicitly appear in the expressions of susceptibilities. At any rate, susceptibilities retain an implicit dependence on Γ through the matrix elements. Then, the right procedure to compute both correlation functions and susceptibilities is clearly the one that starts from the causal Green's function, which is the only Green's function that explicitly depends on Γ. It is worth noticing that the value of Γ dramatically affects the values of directly measurable quantities (e.g., compressibility, specific heat, magnetic susceptibility, . . . ) through the values of correlation functions and susceptibilities. According to this, whenever it is possible, Γ should be exactly calculated case by case.
If we do not have access to the complete set of eigenstates and eigenvalues of the system, which is the rule in the most interesting cases, we have to compute correlation functions and susceptibilities within some, often approximated, analytical framework. Now, since no analytical tool can easily determine Γ (e.g., the equations of motion cannot be used to fix Γ as it is constant in time), one usually assumes the ergodicity of the dynamics of ψ and simply substitutes Γ by its ergodic value (i.e., by the r.h.s. of (1.14)): Unfortunately, this procedure cannot be justified a priori (i.e., without computing Γ through its definition (2.58)) by absolutely no means. The existence of just one integral of motion and, more generally, of any operator that has a diagonal part with respect to the Hamiltonian [ 6] (i.e. by any operator that has a diagonal entries whenever written in the basis of eigenstates of the Hamiltonian) divides the phase space into separate subspaces not connected by the dynamics. This latter, in turn, becomes non ergodic: time averages give different results with respect to ensemble averages. This latter consideration also clarifies why the ergodic nature of the dynamics of an operator mainly depends on the Hamiltonian it is subject to. It is really remarkable that Γ is directly related to relevant measurable quantities such as compressibility and specific heat trough the dissipation-fluctuation theorem. For instance, we recall the formula that relates the compressibility to the total particle number fluctuations whereN is the total particle number operator, N is its average and V is proportional to the volume of the system. We see that a compressibility different from zero requires the non-ergodicity of the system with respect to total particle number operator. According to this, in the case of infinite systems too the correct determination of Γ cannot be considered as an irrelevant issue (e.g., (1.19) holds in the thermodynamic limit too).
In the next section, we provide some examples of violation of the ergodic condition (1.14). It is necessary pointing out, in order to avoid any possible confusion to the reader, that we are using full operators and not fluctuation ones (i.e., we use operators not diminished of their average value, in contrast with what it is usually done for the bosonic excitations like spin, charge and pair). According to this, the Γ can be different from zero (i.e., be equal to the squared average of the operator), and still indicate an ergodic dynamics for the operator.

Two-site Hubbard model
The two-site Hubbard model is described by the following Hamiltonian where the summation range only over two sites at distance a from each other and the rest of notation is standard [ 4]. The hopping matrix t ij is defined by where α(k) = cos(ka) and k = 0, π/a. We now proceed to study the system by means of the equation of motion approach and the Green's function formalism [ 17]. A complete set of fermionic eigen-operators of H is the following one We define ψ α (i) = j α ij ψ(j) and use the spinorial notation for the field operators.
is the charge (µ = 0) and spin (µ = 1, 2, 3) operator; greek (e.g., µ, ν) and latin (e.g., a, b, k) indices take integer values from 0 to 3 and from 1 to 3, respectively; sum over repeated indices, if not explicitly otherwise stated, is understood; σ µ = (1, σ) and σ µ = (−1, σ); σ are the Pauli matrices. In momentum space the field ψ(i) satisfies the equation of motion where the energy matrix ε(k) has the expression Straightforward calculations, according to the scheme traced in [ 17], show that two correlators Then, the Green's functions depend on three parameters: µ, ∆ and p. The correlator ∆ can be expressed in terms of the fermionic correlation function C(i, j) = ψ(i) ψ † (j) ; the chemical potential µ can be related to the particle density by means of the relation n = 2 The parameter p cannot be calculated in the fermionic sector; it is expressed in terms of correlation functions of the bosonic fields n µ (i) and c ↑ (i) c ↓ (i). According to this, the determination of the fermionic Green's functions requires the parallel study of bosonic Green's functions.
After quite cumbersome calculations, it is possible to see [ 17] that a complete set of bosonic eigenoperators of H in the spin-charge channel is given by with the definitions: The field B (µ) (i) satisfies the equation of motion where the energy matrix κ(k) has the expression The energy spectra are given by Straightforward calculations show that the correlation function has the expression where f (n,µ) (0) = 0 f or n = 3, 4, 5, 6 (2.29a) Owing to the fact that zero-energy modes appear for n = 1, 2 and k = 0 [cfr.
Eq. (2.21)], Γ appear in the correlation functions One might think, as is often done in the literature, to fix this constant by its ergodic value. However, this is not correct as we are in a finite system in the grandcanonical ensemble and the ergodicity condition does not hold. For the moment, we can state that this constant remains undetermined. The spectral density functions depends on a set of parameters which come from the calculation of the normalization matrix I (µ) (k) = F B (µ) (i, t), B (µ) † (j, t) . In particular, for the (1,1)-component the following parameters appear: The parameters C α and C α 12 are related to the fermionic correlation function C(i, j) = ψ(i) ψ † (j) . The parameter χ α s can be expressed in terms of the bosonic correlation function C (µ) (i, j) = B (µ) (i) B (µ) † (j) . In order to use the standard procedure of self-consistency, we need to calculate the parameter d. For this purpose we should open both the pair channel and a double occupancy-charge channel (i.e., we will need the static correlation function n ↑ (i) n ↓ (i) n α (i) ). The corresponding calculations are reported in Ref. [ 17] where is shown that these two channels do not carry any new unknown Γ. The self-consistence scheme closes; by considering the four channels (i.e., fermionic, spin-charge, pair and double occupancy-charge) we can set up a system of coupled self-consistent equations for all the parameters. However, Γ (µ) (0) has not been determined yet: we have not definitely fixed the representation of the Green's functions.
In conclusion, the standard procedure of self-consistency is very involved and is not able to give a final answer because of the problem of fixing the Γ. We will now approach the problem by taking a different point of view. The proper representation of the Green's functions must satisfy the condition that all the microscopic laws, expressed as relations among operators must hold also at macroscopic level as relations among matrix elements. For instance, let us consider the fermionic channel. We have seen that there exists the parameter p, not explicitly related to the fermionic propagator, that can be determined by opening other channels. However, we know that at the end of the calculations, if the representation is the right one, the parameter p must take a value such that the symmetries are conserved. By imposing the algebra constraints (??) and by recalling the expression for ∆ we get three equations This set of coupled self-consistent equations will allow us to completely determine the fermionic Green's functions. Calculations show [ 17] that this way of fixing the representation is the right one: all the symmetry relations are satisfied and all the results exactly agree with those obtained by means of Exact Diagonalization. We do not have to open the bosonic channels; the fermionic one is self-contained.
Next, let us consider the spin-charge Green's functions. In the spin-charge sector we have the parameters C α , C α 12 , χ α s , d and two Γ   Fig. 1 for various temperatures. It is worth noting that they assume their ergodic values (i.e. n 2 and 0, respectively) only in some regions of the parameter space: (at zero temperature) at n = 1 (both b 0 and b k ) and at n = 0.5 (b 0 only). In these regions, the grand-canonical ensemble is equivalent to the microcanonical one and the underlying ergodicity of the charge and spin dynamics emerges.
It is worth noting that b 0 is directly related to the compressibility by means of the following relation [ 17] According to this, if we erroneously set the value of b 0 to the ergodic one (i.e., n 2 ) we would get a constant zero compressibility.

Tight-binding model
A narrow-band Bloch system in presence of an external magnetic field is described by the following Hamiltonian where n 3 (i) is the third component of the spin density operator and h is the intensity of the external magnetic field. The indices i and j run on an infinite d-dimensional lattice. Straightforward calculations show that the causal Green's function G (µ) C (i, j) = T [n µ (i) n µ (j)] and the correlation function C (µ) (i, j) = n µ (i) n µ (j) of the chargespin operator n µ (i) = c † (i) σ µ c(i) have the following expressions comes from the proper fermionic loop and is the Fourier transform of Here G C (i, j) = T c(i) c † (j) is the causal fermionic function and has the expression Γ µ is fixed by the algebra constraints (??) which requires The loop Q (µ) (k, ω) can be calculated by means of (2.45). Calculations show By recalling the algebra constraints (2.39), Eq. (2.50) gives for the Γ in accordance with the ergodic nature of the spin and charge dynamics in this system. It is worth noting that the compressibility of this system can be computed by means of the general formula (1.19) that holds in the thermodynamic limit too and gives We can see that an ergodic charge dynamics can lead to a non-ergodic value of the Γ relatively to the total number operator, which is an integral of motion. Also in the infinite systems the decoupling inspired by the requirement of ergodicity cannot always be applied.

Heisenberg chain
We will now study [ 18,19] the ergodicity of the dynamics of the operator S z i , the z-component of the spin at site i, in the 1D anisotropic extended Heisenberg model described by the following Hamiltonian: where S x i , S y i and S z i are the x, y and z components of the spin-1/2 at site i, respectively. The model (2.57) is taken on a linear chain with periodic boundary conditions. We take the interaction term parameterized with J z ferromagnetic (J z > 0) and the next-nearest-neighbor interaction term, which is parameterized with J ′ , isotropic. In order to frustrate ferromagnetism, we have considered only the case with J ′ > 0, that is, with an antiferromagnetic coupling between next-nearest neighbors. According to this, only chains with even number of sites have been studied in order to avoid topological frustration that would be absent in the thermodynamic limit. Since it is possible to exactly map all results obtained for J ⊥ > 0 to those for J ⊥ < 0 by means of a simple canonical transformation, we have limited our study only to positive values of J ⊥ .
We have numerically diagonalized the Hamiltonian (2.57) for chains of size L ranging between 6 and 18 by means of Exact Diagonalization (ED) (divide and conquer algorithm) and for chains of size L ranging between 20 and 26 by means of Lanczos Diagonalization (LD). We have systematically taken into account translational symmetry and classified the eigenstates by the average value of S z = i S z i , which is a conserved quantity. Whenever we have used ED, all eigenvalues and eigenvectors of (2.57) have been calculated up to machine precision and, therefore, we have been able to determine the exact dynamics of the system for all temperatures. On the contrary, when we have used LD, we have been limited to the zero-temperature case since only the ground state can be considered exact in LD.
In this case, we have the opportunity to exactly compute Γ in terms of the exact eigenvalues E n and eigenstates |n of the system. As a matter of fact, it read as [ 4] En=Em e −βEn n|S z i |m m|S z i |n (2.58) As already discussed above, the dynamics of an operator (e.g., S z i ) is ergodic whenever (1.14) is satisfied, or equivalently, (2.58) is equal to its ergodic value: The dynamics of a finite system is hardly ergodic, since (2.58) and (2.59) unlikely coincide. In the thermodynamic limit, the sums in (2.58) and (2.59) become series and no conclusion can be drawn a priori. Since we have diagonalized the Hamiltonian (2.57) numerically (i.e., only for finite systems) and since L → ∞ is the most interesting case, we have analyzed our results through finite-size scaling in order to speculate on the properties of the bulk system. If the ground state of (2.57) is N-fold degenerate then, at T = 0, (2.58) and (2.59) read as follows: respectively. Thanks to the translational invariance enjoined by the system S z i is independent of i and proportional to the z-component of the total spin operator average S z tot . It is easy to show that, even if there is a finite magnetic moment per site in any of such N degenerate ground states, S z i at T = 0 is always zero in absence of magnetic field. Indeed, if a ground state with non-zero S z i = M exists, also another ground state with S z i = −M exists. Thus, at zero temperature, Γ erg is always zero and the only quantity of interest is Γ. A finite value of this latter implies non-ergodicity. Obviously, if N = 1 then both values coincide. Therefore, a non-ergodic phase corresponds to degenerate ground states with finite magnetization.
In the studied range of coupling constants (see Fig. 2) we have found two nonergodic phases (NE-I and NE-II), two ergodic ones (E-I and E-II) and a weird phase (W). Our computational facilities limit the range of chain sizes that we can analyze such that we could not establish, by means of finite-size scaling, whether the weird phase (W) is ergodic or not. In the non-ergodic phases (NE-I and NE-II), we were able not only to perform the finite-size scaling, but also to write down an analytic expression for Γ as a function of the chain size L. The weird phase (W) has exhibited a strong dependence of the ground state upon the particular values of the couplings. On the contrary, the other phases exhibit ground states that are independent of the particular values of the coupling constants.
In the standard Heisenberg model (J ′ = 0 and J ⊥ = J z ) at T = 0 the dynamics is non-ergodic for ferromagnetic coupling (J ⊥ = J z < 0) as the system has a L + 1 degenerate ground state It is clear from (2.60) that Γ remains non-ergodic also in the thermodynamic limit. This point (J ′ = 0 and J ⊥ = J z ) becomes a line in our phase diagram and is denoted as NE-I (see Fig. 2). In fact, the next-nearest-neighbor interaction J ′ may frustrate (J ′ > 0) or favor (J ′ < 0) the ferromagnetism. In the latter case, the ground state remains unchanged for any value of J ′ < 0. Therefore, we expect the line denoting the phase NE-I to extend also to negative J ′ . If, on the contrary, J ′ is positive and large enough to frustrate the system in such a way that the ground state loses its ferromagnetic character, the ergodicity is restored. This occurs at a finite critical J ′ ∼ 0.25J z . For values of J ′ larger than the critical one, we find a non-degenerate ground state with S z tot = 0. If J ⊥ = J z the rotational invariance is broken so that states with the same S 2 tot , but different S z tot , are not degenerate anymore. In the non-ergodic region (NE-II) of the phase diagram (see Fig. 2), the ground state is just doubly degenerate (not L + 1 degenerate as in (NE-I)): one ground state corresponds to a configuration with all spins up and the other to a configuration with all spins down. Hence, the value of Γ in this phase is 1/4 and does not depend neither on the Hamiltonian couplings nor on the number of sites in the chain. It is clear that also this phase extends to negative values of J ′ . This kind of ground state stands the frustration introduced by next-nearest-neighbor interaction up to J ′ ∼ 0.3J z (see Fig. 2).
The ergodic region (E-I) of the phase diagram (see Fig. 2) has Γ = 0 for all sizes of the system and values of the couplings: the unique ground state belongs to the sector with S z tot = 0. On the contrary, the other ergodic phase (E-II) has non-zero values of Γ for values of L not multiples of four. The ground state in this phase has average total spin equal to one and, therefore, Γ = 1/L 2 . We obviously conclude that (E-II) phase is ergodic in the thermodynamic limit.
The values of Γ in these four phases (NE-I, NE-II, E-I and E-II) exhibit perfect finite-size scaling as shown in Fig. 3a). This has allowed us to make definite statements also in the thermodynamic limit.
The weird phase (W) (see Fig. 2) is characterized by a quite strong size dependence, as shown in Fig. 3b) where a tentative finite-size scaling of Γ in the different points of the phase is presented. This region manifests a diverging finite-size scaling within the range of sizes we were able to handle. In this case, the behavior of Γ as a function of L strongly depends on the particular choice of the Hamiltonian couplings and is highly non monotonous when increasing L, according to the strong dependence on L of S z tot in the ground state. In this critical region the eigenvalues of (2.57) present many level crossings, which means that the maximum value of L we were able to reach (L max = 26) is not large enough to perform a sensible finite-size scaling analysis. However, we expect that this phase becomes ergodic in the thermodynamic limit, although still different from the ergodic phases E-I and E-II.
We can summarize our findings in the thermodynamic limit at zero temperature as follows:

Conclusions
In conclusion, we have analyzed the issue of ergodicity, after a brief historical overview, within the Green's function formalism by means of the equations of motion approach. We have individuated the primary source of non-ergodic dynamics for a generic operator in the appearance of zero-frequency anomaly in its correlation functions and given a recipe to compute the unknown quantities characterizing such a behavior within the Composite Operator Method. Finally, we have presented examples of non-trivial strongly correlated systems where it is possible to examine a non-ergodic behavior: two-site Hubbard model, tight-binding model, Heisenberg chain.