# Boltzmann relaxation dynamics in the strongly interacting Fermi-Hubbard model

###### Abstract

Via the hierarchy of correlations, we study the Mott insulator phase of the Fermi-Hubbard model in the limit of strong interactions and derive a quantum Boltzmann equation describing its relaxation dynamics. In stark contrast to the weakly interacting case, we find that the scattering cross sections strongly depend on the momenta of the colliding quasi-particles and holes. Therefore, the relaxation towards equilibrium crucially depends on the spectrum of excitations. For example, for particle-hole excitations directly at the minimum of the (direct) Mott gap, the scattering cross sections vanish such that these excitations can have a very long life-time.

## I Introduction

The laws of thermodynamics are very powerful tools in physics with far reaching consequences. However, understanding the microscopic origin of thermal behavior can be a very challenging question – which is also the origin of the famous debate between Loschmidt and Boltzmann L1876 ; B1877 ; B1872 . For classical many-body systems, the relaxation to a thermal equilibrium state is typically understood in terms of an effective description in the form of a Boltzmann equation B75 . When and where such an effective description is adequate can still be a non-trivial question KWW06 ; GME11 ; BCH11 ; PSSV11 ; RDO08 ; Getal12 ; Ketal16 ; Netal16 , related to the BBGKY hierarchy K46 ; B46 ; BG46 and chaotic versus integrable behavior.

For quantum many-body systems, the question of whether and how these systems relax to a thermal equilibrium state can be even more involved and is being widely discussed in the literature, see, e.g., EKW09 ; EKW10 ; WEFWBH18 ; PSP18 ; D91 ; S94 ; RDO02 ; CR10 ; RS12 ; R13 . For example, the interplay between disorder and interactions can have a non-trivial impact on the relaxation dynamics, see, e.g., BAA06 ; CRFSS11 ; NH15 . In the following, we focus on closed quantum lattice systems without disorder and dissipation, whose unitary dynamics describes thermalization induced by the intrinsic interactions. Still, their relaxation and thermalization dynamics can show non-trivial features, e.g., it can undergo several stages with different time scales, see, e.g., BBS03 ; BBW04 ; KWE11 .

The thermalization of weakly interacting quantum many-body systems is typically understood in terms of a quantum version of the Boltzmann equation, derived by means of suitable approximation schemes such as the Born-Markov approximation BP02 ; RK02 .

There are several investigations for one-dimensional systems, see, e.g., MWNM98 ; R09a ; R09b ; RDYO07 ; BKL10 ; KISD11 ; SVPH14 . However, due to energy and momentum conservation and potential further conservation laws (chaotic versus integrable behavior), the relaxation dynamics in one dimension displays peculiar features and is qualitatively different from that in higher dimensions. Thus, these one-dimensional systems are of limited help for understanding higher dimensional cases.

## Ii The Model

In order to start filling this gap, we consider the Fermi-Hubbard Hamiltonian as a prototypical model for strongly interacting fermions which move on a regular lattice given by the hopping matrix and repel each other via the local interaction

(1) |

As usual, and are the fermionic creation and annihilation operators for the lattice sites and and the spin with the corresponding number operators . Furthermore, denotes the coordination number of the translationally invariant lattice, i.e., the number of nearest neighbors.

In one spatial dimension, the Fermi-Hubbard Hamiltonian (1) is integrable via the Bethe ansatz LW68 and thus would not display full thermalization in view of the infinite number of conserved quantities (in addition to the impossibility of thermalization via two-body collisions due to energy and momentum conservation, as mentioned in the Introduction). Thus, we focus on higher-dimensional lattices (with large ) in the following.

In the limit of small interactions , the ground state of (1) can be described by a Fermi gas and is thus metallic for . For large interactions , however, the structure of the ground state changes. Assuming half filling , the repulsion generates a gap and we obtain the Mott insulator state containing one fermion per site (plus virtual tunneling corrections), cf. H63 ; IFT98 .

## Iii Hierarchy of Correlations

For weak interactions , a perturbative expansion in allows us to simplify the equations of motion and to justify the Markov approximation (see the Appendix). For strong interactions , however, this procedure is no longer applicable and thus one has to find an alternative approach.

Here, we employ the hierarchy of correlations NS10 ; QNS12 ; QKNS14 ; KNQS14 ; NQS14 ; NQS16 ; QS19 and consider the reduced density matrices for one site and for two sites etc. After splitting off the correlations via and so on, we obtain the following hierarchy of evolution equations NS10

(2) | |||||

(3) | |||||

(4) | |||||

(5) |

and in complete analogy for the higher correlators.

In order to truncate this infinite set of recursive equations, we exploit the hierarchy of correlations in the formal limit of large coordination numbers . With the arguments outlined in NS10 , it can be shown that the two-site correlations are suppressed via in comparison to the on-site density matrix . Furthermore, the three-site correlators are suppressed even stronger via , and so on. This hierarchy of correlations facilitates the following iterative approximation scheme: To zeroth order in , we may approximate (2) via which yields the mean-field solution . As the next step, we may insert this solution into (3) and obtain to first order in the approximation which gives a set of linear and inhomogeneous equations for the two-point correlations . From this set, we obtain the quasi-particle excitations and their energies.

Since this set of equations is linear in , it does not describe interactions between the quasi-particles and hence we do not obtain a Boltzmann collision term to first order in . To this end, we have to go to higher orders in and study the impact of the three-point correlators in (3). As one might already expect from the well-known derivation for weak interactions (see the Appendix), it is not sufficient to truncate the set of equations (2)-(5) at this stage – we have to include the four-point correlators in order to derive the Boltzmann equation (see below).

Finally, the back-reaction of the quasi-particle fluctuations onto the mean field can be derived by inserting the solution for back into equation (2).

## Iv Mott Insulator State

As explained above, the starting point of the hierarchy is the on-site density matrix or its zeroth-order (mean-field) approximation . Assuming a spatially homogeneous state at half filling footnote , we get the simple solution of equation (2)

(6) |

where denotes the double occupancy and measures the deviation from the ideal Mott insulator state for .

Now we may insert this solution into Eq. (3) and study the two-point correlations. In order to describe the relevant correlators describing the dynamics of the quasi-particles (also called doublons) and holes (or holons), we introduce the short-hand notation which is just for but for (see the Appendix). Then we may define the uppercase operators via

(7) |

where is the spin index opposite to . For , they correspond to the annihilation of a fermion with spin at the lattice site when there is another fermion with opposite spin at that site. Thus, this case corresponds to a quasi-particle (doublon) excitation. In analogy, the case corresponds to the absence of another fermion with opposite spin at that site, i.e., a hole (holon) excitation.

In terms of these operators (7), the quasi-particle and hole correlators can be written as

(8) |

where denotes the difference between the positions and of the lattice sites and . Here, we have assumed spatial homogeneity. In principle, one could also consider inhomogeneous excitations, where these functions which enter the Boltzmann equation would acquire an additional position coordinate, i.e., instead of . Then, the Boltzmann equation would also contain terms describing the propagation of the excitations. However, here we are mainly interested in the collision terms in the Boltzmann equation and hence we assume spatial homogeneity for simplicity.

## V Dispersion Relation

In terms of the , the evolution equation for the two-point correlators (8) obtained from Eq. (3) reads

(9) |

where the source term contains the three-point correlators and is suppressed as . Apart from this source term, the set of equations (V) is linear and can be can be diagonalized by means of an orthogonal transformation matrix , see the Appendix. We denote the transformed (rotated) correlation functions by lowercase superscripts via . Thus, the set of equations (V) simplifies to

(10) |

with the quasi-particle () and hole () energies LPM69

(11) |

The functions are rapidly oscillating for but slowly varying for because of . Thus, the -expansion (hierarchy of correlations) employed here naturally provides a separation of time scales: We have rapidly varying quantities whose rate of change is given by the eigen-energies (11) or linear combinations thereof, while the rate of change of the slowly varying quantities is suppressed with (or even higher). As in the weakly interacting case, this separation of time-scales will be used to justify the Markov approximation.

In the (Mott insulating) ground state, these correlation functions assume the values , and , see, e.g., QS19 . Hence any deviation from these values indicates a departure from the ground state, i.e., an excitation. As a result, the correlation functions determine the excitations present in our system. Accordingly, we denote the slowly varying quantities as our quasi-particle distribution functions for () with and the hole distribution function for () with .

## Vi Higher Correlations

As shown above, the rate of change of is determined by the source term containing the three-point correlation functions

(12) | ||||

(13) | ||||

(14) |

which are of order . The evolution equations for these correlators (12)-(13) can be derived from equation (4) and read after the rotation with (see the Appendix)

(15) | ||||

(16) | ||||

(17) |

The source terms , , and in the above equations (15)-(16) contain various combinations of two-point correlators and the four-point correlators which are indispensable for the Boltzmann collision terms

(18) |

Finally, their evolution equation can be derived from Eq. (5). After a rotation with , we find (see the Appendix)

(19) |

where the source term contains three-point and two-point correlations as well as terms of higher order in , such as the five-point correlator, which we neglect.

## Vii Markov Approximation

In order to arrive at a time-local Boltzmann equation, the differential equations (15)-(16) and (VI) are integrated within the Markov approximation. All these equations are of the general form and thus have formally the solution

(20) |

The source terms containing the distribution functions are slowly varying, with their rate of change being suppressed by or even more, in comparison with the rapid oscillations . Hence we may approximate in the above integral (20) which gives

(21) |

with the infinitesimal shift selecting the retarded solution. As usual, this Markov approximation effectively neglects memory effects. It allows the elimination of all three-point and four-point correlators such that finally only the slowly varying distribution functions remain. After some algebra (see the Appendix) we arrive at

(22) |

This is the quantum Boltzmann equation and represents our main result. It has the same general form as in the weakly interacting case. Let us first discuss the common features. The describe the scattering cross sections for the various processes. For example, corresponds to the collision of two quasi-particles with initial momenta and , which are scattered to the final momenta and , thus satisfying momentum conservation (with the momentum transfer ). Energy conservation is incorporated via the Dirac delta function in the second line of Eq. (22). The last line of Eq. (22) corresponds to the inverse process, which ensures the conservation of probability.

As another analogy to the weakly interacting case, the structure of the last two lines of Eq. (22) reflects the fermionic character of the quasi-particles and holes. (For bosons, one would have instead of .) Related to this fermionic nature is the particle-hole duality where the distribution function describing quasi-particles is mapped to the distribution function of the holes. Thus, in addition to processes such as the collision between two quasi-particles or two holes or a quasi-particle with a hole , the above equation (22) does in principle also contain processes: E.g., corresponds to the inelastic scattering of one quasi-particle via the simultaneous creation of a new particle-hole pair (or the inverse process). However, here we are mainly interested in the strongly interacting limit , where such processes are forbidden by energy conservation: The initial particle energy is not large enough to create a final state with an energy of nearly .

As the final analogy to the weakly interacting case, we note that only quasi-particles (or holes) of opposite spins and scatter, at least to the leading order considered here. For weak interactions, this is a simple consequence of the structure of the on-site interaction term , but for strong interactions, the situation is a bit more complex (see below).

## Viii Strongly Interacting Limit

As the most crucial difference to the weakly interacting case, the scattering cross sections acquire a non-trivial momentum dependence. To illustrate this, let us consider the limit of strong interactions . In this limit, the Boltzmann equation (22) describing collisions of two quasi-particles simplifies to

(23) | |||||

For the collision of two holes, the equation has the same form after replacing all the with . The equations describing the collision of a quasi-particle and a hole have a very similar structure (see the Appendix).

For weakly interacting systems (see the Appendix), the scattering cross section is momentum independent and given by . Here, we find that the interaction does not occur in the Boltzmann equation (23) at all, where the scattering cross section reads and is thus depends on the momenta and of the incoming quasi-particles. This difference can be understood in terms of the following simplified and intuitive picture: In the Mott insulator state, all lattice sites are occupied by one fermion and thus a quasi-particle roughly corresponds to a doubly occupied lattice site. As a consequence, two quasi-particles cannot occur at the same lattice site and thus they cannot directly interact via the strong on-site repulsion . Instead, they can “feel” each other via virtual tunneling processes (which are Pauli blocked if the neighboring lattice site is also occupied by a quasi-particle). These virtual tunneling processes explain the scaling with and the momentum dependence.

This momentum dependence can have strong implications for the relaxation dynamics: If we consider momentum conserving excitation process such as a long-wavelength pump laser, the energy cost of creating a particle-hole pair is given by the direct gap

(24) |

which assumes its minimum value at those points where vanishes. Now, a weak enough pump laser with a frequency sufficiently below the gap would predominantly create excitations near those minimum-energy wave-numbers where . On the other hand, for these quasi-particle excitations, the scattering cross sections in the Boltzmann equation (23) vanish and thus they would relax very slowly. This behavior is also shown by the other channels (such as particle-hole collisions) in the strongly interacting limit.

## Ix Back-Reaction

Finally, via inserting the correlation functions back into equation (2), we may calculate the back-reaction of the quasi-particle and hole fluctuations onto the mean field . This determines the double occupancy in Eq. (6) via

(25) |

However, this small double occupancy does not affect our leading-order results, such as the scattering cross sections in the Boltzmann equation (23).

## X Conclusions and Outlook

As a prototypical example for strongly interacting quantum many-body system on a lattice, we consider the Fermi-Hubbard model (1) in the Mott insulator state. Via the hierarchy of correlations, we derive a quantum Boltzmann equation (22) describing the relaxation dynamics of the quasi-particle (doublon) and hole (holon) excitations. As the most crucial difference to the weakly interacting case, we find that the scattering cross sections display a strong momentum dependence, cf. Eq. (23), which has profound consequences for the relaxation dynamics. In analogy to the weakly interacting case, the Boltzmann equation (23) facilitates the derivation of an -theorem.

Our method can be generalized to other lattice systems, such as the Bose-Hubbard model or spin lattices KS18 ; WT11 . It can also be used to study higher-order correlators such as the spin modes in the Fermi-Hubbard model (such as with ), which are of bosonic nature. Considering the extended Fermi-Hubbard model including long-range Coulomb interactions, one would expect that they generate additional scattering cross sections in the Boltzmann equation (23) and thus also influence the relaxation dynamics.

###### Acknowledgements.

This work was funded by DFG, grant # 278162697 (SFB 1242) and 398912239.## Appendix A Boltzmann equations for weakly interacting fermions

For weakly interacting fermions, the Boltzmann evolution equation can be derived via time-dependent perturbation theory. The Hamiltonian for interacting fermions reads

(26) |

where and are spin indices and denotes the interaction potential. In order to apply perturbation theory, we shall transform (26) to Fourier space in order to diagonalize the kinetic part. Note that the hierarchical expansion starts from the atomic limit and the hopping Hamiltonian introduces the correlation between lattice sites, see below.

The Hamiltonian (26) has the Fourier representation

(27) |

from which one can obtain the equation of motion of the fermion distribution function , i.e.,

(28) |

As can be seen from (A), the dynamics is solely governed by the correlation functions

(29) |

The equation of motion of the correlators (A) can be integrated within Markov approximation. Substituting the resulting expression into (A) we arrive at

(30) |

## Appendix B Boltzmann equations for the strongly interacting Hubbard model

It is clear that for strongly interacting systems, the derivation of the Boltzmann dynamics cannot be based on an expansion in powers of the interaction strength between the electrons. As explained in the paper, we employ therefore a hierarchical expansion for large coordination numbers .

In the following we give a step-by-step derivation of the Boltzmann kinetic equation (22). We consider the simplest possible case and assume that the system is always in an unpolarized state at half filling which is metallic for and insulating for . We demand that the initial state has -symmetry, such that the density matrix commutes with for all times.

### b.1 Operator equations

We introduce a compact notation in order to make the calculation tractable. Therefore we define the operators

(31) | ||||

(32) |

where denotes the lattice site and is the spin index. Using the Heisenberg equations for the Hubbard Hamiltonian (1), we find

(33) |

and

(34) |

with and . The operator evolves according to

(35) |

the spin-flip operator satisfies the equation

(36) |

and the doublon creation (annihilation) operators have the equation of motion

(37) |

and

(38) |

In the following we shall use the above relations to evaluate the evolution equations of the hierarchical correlation functions.

### b.2 Double occupancy and two-site correlation functions

Due to the -symmetry, any expectation value which contains an odd number of creation operators and annihilation operators for a fixed spin index vanishes identically. This implies, for example,

(39) |

The zeroth order equation of the hierarchical expansion (2) determines the double-occupancy , i.e.

(40) |

From the first order equation (3) follows the dynamics of the two-point correlation functions