Boltzmann relaxation dynamics in the strongly interacting Fermi-Hubbard model

TL;DR

This paper derives a quantum Boltzmann equation for the Mott insulator phase of the strongly interacting Fermi-Hubbard model, revealing momentum-dependent scattering and long-lived excitations, contrasting with weakly interacting systems.

Contribution

It introduces a novel derivation of a Boltzmann equation for the strongly interacting Fermi-Hubbard model's Mott phase, highlighting momentum-dependent scattering effects.

Findings

Scattering cross sections depend strongly on quasi-particle momenta.

Excitations at the Mott gap minimum have vanishing scattering cross sections.

Relaxation dynamics are significantly influenced by the excitation spectrum.

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.

Equations

$$\hat{H}=-\frac{1}{Z}\sum _{\mu ,\nu ,s}{J}_{\mu \nu }{\hat{c}}_{\mu ,s}^{†}{\hat{c}}_{\nu ,s}+U\sum _{\mu }{\hat{n}}_{\mu }^{\uparrow }{\hat{n}}_{\mu }^{\downarrow }.$$

$${\partial }_t{\hat{\rho }}_{\mu }$$

$${\partial }_t{\hat{\rho }}_{\mu \nu }^{\mathrm{corr}}$$

$${\partial }_t{\hat{\rho }}_{\mu \nu \sigma }^{\mathrm{corr}}$$

$${\partial }_t{\hat{\rho }}_{\mu \nu \sigma \lambda }^{\mathrm{corr}}$$

$${\hat{\rho }}_{\mu }=\left(\frac{1}{2}-\mathfrak{D}\right)\left(|\uparrow ⟩⟨\uparrow |+|\downarrow ⟩⟨\downarrow |\right)+\mathfrak{D}\left(|\uparrow \downarrow ⟩⟨\uparrow \downarrow |+|0⟩⟨0|\right),$$

$${\hat{C}}_{\mu ,s}^X={\hat{c}}_{\mu ,s}{\hat{N}}_{\mu ,\bar{s}}^X,$$

$$f_{\mu \nu ,s}^{XY}=⟨({\hat{C}}_{\mu ,s}^X{)}^{†}{\hat{C}}_{\nu ,s}^Y⟩=\underset{\mathbf{k}}{\int }{f}_{\mathbf{k},s}^{XY}\exp \{i\mathbf{k}\cdot \Delta {\mathbf{r}}_{\mu \nu }\},$$

$$i{\partial }_t{f}_{\mathbf{k},s}^{XY}$$

$$+S_{\mathbf{k},s}^{XY},$$

$$i{\partial }_t{f}_{\mathbf{k},s}^{ab}=(E_{\mathbf{k}}^b-E_{\mathbf{k}}^a)f_{\mathbf{k},s}^{ab}+2S_{\mathbf{k},s}^{ab},$$

$$E_{\mathbf{k}}^{\pm }=\frac{1}{2}\left(U-J_{\mathbf{k}}\pm \sqrt{J_{\mathbf{k}}^2+U^2}\right).$$

$$⟨{{\hat{N}}_{\rho ,\bar{s}}}^X({\hat{C}}_{\mu ,s}^Y{)}^{†}{\hat{C}}_{\nu ,s}^Z{⟩}^{\mathrm{corr}}$$

$$⟨{\hat{c}}_{\rho ,s}^{†}{\hat{c}}_{\rho ,\bar{s}}({\hat{C}}_{\mu ,\bar{s}}^X{)}^{†}{\hat{C}}_{\nu ,s}^Y{⟩}^{\mathrm{corr}}$$

$$⟨{\hat{c}}_{\rho ,s}^{†}{\hat{c}}_{\rho ,\bar{s}}^{†}{\hat{C}}_{\mu ,\bar{s}}^X{\hat{C}}_{\nu ,s}^Y{⟩}^{\mathrm{corr}}$$

$$i{\partial }_t{G}_{\mathbf{pq},\bar{s}ss}^{Xab}$$

$$i{\partial }_t{I}_{\mathbf{pq},\bar{s}s}^{ab}$$

$$i{\partial }_t{H}_{\mathbf{pq},\bar{s}s}^{ab}$$

$$⟨({\hat{C}}_{\alpha ,\bar{s}}^X{)}^{†}{\hat{C}}_{\beta ,\bar{s}}^Y({\hat{C}}_{\mu ,s}^V{)}^{†}{\hat{C}}_{\nu ,s}^W{⟩}^{\mathrm{corr}}=\underset{\mathbf{p},\mathbf{q},\mathbf{k}}{\int }{J}_{\mathbf{pqk},\bar{s}\bar{s}ss}^{XYVW}\times$$

$$e^{i\mathbf{p}\cdot \Delta {\mathbf{r}}_{\beta \alpha }+i\mathbf{q}\cdot \Delta {\mathbf{r}}_{\mu \alpha }+i\mathbf{k}\cdot \Delta {\mathbf{r}}_{\nu \alpha }}.$$

$$i{\partial }_t{J}_{\mathbf{pqk},\bar{s}\bar{s}ss}^{abcd}$$

$$+S_{\mathbf{pqk},\bar{s}\bar{s}ss}^{abcd}$$

$$C(t)=-i\underset{-\infty }{\overset{t}{\int }}d{t}^{\prime }S(t^{\prime })e^{-i\Omega (t-t^{\prime })}.$$

$$C(t)\approx -\frac{S(t)}{\Omega -iϵ},$$

$${\partial }_t{f}_{\mathbf{k},s}^d$$

$$\delta \left(E_{\mathbf{p}+\mathbf{q}}^a-E_{\mathbf{p}}^b+E_{\mathbf{k}-\mathbf{q}}^c-E_{\mathbf{k}}^d\right)\times$$

$$\displaystyle\big{[}f_{\mathbf{k},s}^{d}f_{\mathbf{p},\bar{s}}^{b}\left(1-f_{\mathbf{k-q},s}^{c}\right)\left(1-f_{\mathbf{p+q},\bar{s}}^{a}\right)$$

$$\displaystyle-f^{c}_{\mathbf{k-q},s}f^{a}_{\mathbf{p+q},\bar{s}}(1-f^{d}_{\mathbf{k},s})(1-f^{b}_{\mathbf{p},\bar{s}})\big{]}.$$

$${\partial }_t{f}_{\mathbf{k},s}^{+}$$

$$\Delta E_{\mathbf{k}}=E_{\mathbf{k}}^{+}-E_{\mathbf{k}}^{-}=\sqrt{J_{\mathbf{k}}^2+U^2},$$

$$i{\partial }_t\mathfrak{D}=\sum _s{\int }_{\mathbf{k}}{J}_{\mathbf{k}}(f_{\mathbf{k},s}^{01}-f_{\mathbf{k},s}^{10}).$$

$$\hat{H}=-\frac{1}{Z}\sum _{\mu ,\nu ,s}{J}_{\mu \nu }{c}_{\mu ,s}^{†}{\hat{c}}_{\nu ,s}+\frac{1}{2Z}\sum _{\mu ,\nu ,s,s^{\prime }}{V}_{\mu \nu }^{s{s}^{\prime }}{\hat{n}}_{\mu ,s}{\hat{n}}_{\nu ,s^{\prime }}$$

$$\hat{H}$$

$$+\frac{1}{2N}\sum _{\mathbf{k},\mathbf{q},\mathbf{p}}\sum _{s,s^{\prime }}{V}_{\mathbf{k}}^{s{s}^{\prime }}{\hat{c}}_{\mathbf{q}+\mathbf{k},s}^{†}{\hat{c}}_{\mathbf{q},s}{\hat{c}}_{\mathbf{p}-\mathbf{k}{s}^{\prime }}^{†}{\hat{c}}_{\mathbf{p},s^{\prime }}$$

$$i{\partial }_t{n}_{\mathbf{k},s}$$

$$\displaystyle-\langle\hat{c}^{\dagger}_{\mathbf{k-q},s}\hat{c}^{\dagger}_{\mathbf{p+q},s^{\prime}}\hat{c}_{\mathbf{p},s^{\prime}}\hat{c}_{\mathbf{k},s}\rangle^{\mathrm{corr}}\bigg{)}\,.$$

$$⟨{\hat{c}}_{{\mathbf{p}}_1,s}^{†}{\hat{c}}_{{\mathbf{p}}_2,s^{\prime }}^{†}{\hat{c}}_{{\mathbf{p}}_3,s^{\prime }}{\hat{c}}_{{\mathbf{p}}_4,s}{⟩}^{\mathrm{corr}}=⟨{\hat{c}}_{{\mathbf{p}}_1,s}^{†}{\hat{c}}_{{\mathbf{p}}_2,s^{\prime }}^{†}{\hat{c}}_{{\mathbf{p}}_3,s^{\prime }}{\hat{c}}_{{\mathbf{p}}_4,s}⟩$$

$$+({\delta }_{s,s^{\prime }}{\delta }_{{\mathbf{p}}_1,{\mathbf{p}}_3}{\delta }_{{\mathbf{p}}_2,{\mathbf{p}}_4}-{\delta }_{{\mathbf{p}}_1,{\mathbf{p}}_4}{\delta }_{{\mathbf{p}}_2,{\mathbf{p}}_3})n_{{\mathbf{p}}_1,s}{n}_{{\mathbf{p}}_2,s^{\prime }}.$$

Full text

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 $J_{\mu\nu}$ and repel each other via the local interaction $U$

[TABLE]

As usual, $\hat{c}^{\dagger}_{\mu,s}$ and $\hat{c}_{\nu,s}$ are the fermionic creation and annihilation operators for the lattice sites $\mu$ and $\nu$ and the spin $s\in\{\uparrow,\downarrow\}$ with the corresponding number operators $\hat{n}_{\mu}^{s}=\hat{c}^{\dagger}_{\mu,s}\hat{c}_{\mu,s}$. Furthermore, $Z$ 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 $Z$) in the following.

In the limit of small interactions $U$, the ground state of (1) can be described by a Fermi gas and is thus metallic for $0<\langle\hat{n}_{\mu}^{s}\rangle<1$. For large interactions $U$, however, the structure of the ground state changes. Assuming half filling $\langle\hat{n}_{\mu}^{s}\rangle=1/2$, the repulsion $U$ 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 $U$, a perturbative expansion in $U$ allows us to simplify the equations of motion and to justify the Markov approximation (see the Appendix). For strong interactions $U$, 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 $\hat{\rho}_{\mu}$ for one site and $\hat{\rho}_{\mu\nu}$ for two sites etc. After splitting off the correlations via $\hat{\rho}_{\mu\nu}^{\rm corr}=\hat{\rho}_{\mu\nu}-\hat{\rho}_{\mu}\hat{\rho}_{\nu}$ and so on, we obtain the following hierarchy of evolution equations NS10

[TABLE]

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 $Z\to\infty$. With the arguments outlined in NS10 , it can be shown that the two-site correlations are suppressed via $\hat{\rho}_{\mu\nu}^{\rm corr}=\,{\cal O}(1/Z)$ in comparison to the on-site density matrix $\hat{\rho}_{\mu}=\,{\cal O}(Z^{0})$. Furthermore, the three-site correlators are suppressed even stronger via $\hat{\rho}_{\mu\nu\sigma}^{\rm corr}=\,{\cal O}(1/Z^{2})$, and so on. This hierarchy of correlations facilitates the following iterative approximation scheme: To zeroth order in $1/Z$, we may approximate (2) via $\partial_{t}\hat{\rho}_{\mu}\approx f_{1}(\hat{\rho}_{\nu},0)$ which yields the mean-field solution $\hat{\rho}_{\mu}^{0}$. As the next step, we may insert this solution $\hat{\rho}_{\mu}^{0}$ into (3) and obtain to first order in $1/Z$ the approximation $\partial_{t}\hat{\rho}_{\mu\nu}^{\rm corr}\approx f_{2}(\hat{\rho}_{\nu}^{0},\hat{\rho}_{\mu\nu}^{\rm corr},0)$ which gives a set of linear and inhomogeneous equations for the two-point correlations $\hat{\rho}_{\mu\nu}^{\rm corr}$. From this set, we obtain the quasi-particle excitations and their energies.

Since this set $\partial_{t}\hat{\rho}_{\mu\nu}^{\rm corr}\approx f_{2}(\hat{\rho}_{\nu}^{0},\hat{\rho}_{\mu\nu}^{\rm corr},0)$ of equations is linear in $\hat{\rho}_{\mu\nu}^{\rm corr}$, it does not describe interactions between the quasi-particles and hence we do not obtain a Boltzmann collision term to first order in $1/Z$. To this end, we have to go to higher orders in $1/Z$ and study the impact of the three-point correlators $\hat{\rho}_{\mu\nu\sigma}^{\rm corr}$ 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 $\hat{\rho}_{\mu}$ can be derived by inserting the solution for $\hat{\rho}_{\mu\nu}^{\rm corr}$ back into equation (2).

IV Mott Insulator State

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

[TABLE]

where $\mathfrak{D}$ denotes the double occupancy and measures the deviation from the ideal Mott insulator state for $U\gg J$.

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 $\hat{N}_{\mu,s}^{X}$ which is just $\hat{n}_{\mu,s}$ for $X=1$ but $1-\hat{n}_{\mu,s}$ for $X=0$ (see the Appendix). Then we may define the uppercase operators via

[TABLE]

where $\bar{s}$ is the spin index opposite to $s$. For $X=1$, they correspond to the annihilation of a fermion with spin $s$ at the lattice site $\mu$ when there is another fermion with opposite spin $\bar{s}$ at that site. Thus, this case $X=1$ corresponds to a quasi-particle (doublon) excitation. In analogy, the case $X=0$ corresponds to the absence of another fermion with opposite spin $\bar{s}$ 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

[TABLE]

where $\Delta\mathbf{r}_{\mu\nu}=\mathbf{r}_{\mu}-\mathbf{r}_{\nu}$ denotes the difference between the positions $\mathbf{r}_{\mu}$ and $\mathbf{r}_{\nu}$ of the lattice sites $\mu$ and $\nu$. 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., $f^{XY}(\mathbf{k},\mathbf{r},s)$ instead of $f^{XY}(\mathbf{k},s)$. Then, the Boltzmann equation would also contain terms $\partial f^{XY}(\mathbf{k},\mathbf{r},s)/\partial\mathbf{r}$ 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 $f^{XY}_{\mathbf{k},s}$, the evolution equation for the two-point correlators (8) obtained from Eq. (3) reads

[TABLE]

where the source term $S_{\mathbf{k},s}^{XY}$ contains the three-point correlators and is suppressed as $1/Z^{2}$. Apart from this source term, the set of equations (V) is linear and can be can be diagonalized by means of an orthogonal $2\times 2$ transformation matrix $O_{X}^{a}(\mathbf{k})$, see the Appendix. We denote the transformed (rotated) correlation functions by lowercase superscripts via $f_{\mathbf{k},s}^{ab}=2\sum_{XY}O_{X}^{a}(\mathbf{k})O_{Y}^{b}(\mathbf{k})f^{XY}_{\mathbf{k}}$. Thus, the set of equations (V) simplifies to

[TABLE]

with the quasi-particle ($a=+$) and hole ($a=-$) energies LPM69

[TABLE]

The functions $f^{ab}_{\mathbf{k},s}$ are rapidly oscillating for $a\neq b$ but slowly varying for $a=b$ because of $S_{\mathbf{k},s}^{XY}=\,{\cal O}(1/Z^{2})$. Thus, the $1/Z$-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 $1/Z$ (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 $f^{XY}_{\mathbf{k},s}$ assume the values $f^{01}_{\mathbf{k},s}=f^{10}_{\mathbf{k},s}=J_{\mathbf{k}}/(4\sqrt{U^{2}+J_{\mathbf{k}}^{2}})$, $f^{00}_{\mathbf{k},s}=1/4+U/(4\sqrt{U^{2}+J_{\mathbf{k}}^{2}})-\mathfrak{D}$ and $f^{11}_{\mathbf{k},s}=1/4-U/(4\sqrt{U^{2}+J_{\mathbf{k}}^{2}})+\mathfrak{D}$, 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 $f^{ab}_{\mathbf{k},s}$ determine the excitations present in our system. Accordingly, we denote the slowly varying quantities $f^{a=b}_{\mathbf{k},s}$ as our quasi-particle distribution functions for ($a=b=+$) with $f^{+}_{\mathbf{k},s}$ and the hole distribution function for ($a=b=-$) with $f^{-}_{\mathbf{k},s}$.

VI Higher Correlations

As shown above, the rate of change of $f^{a=b}_{\mathbf{k},s}$ is determined by the source term $S_{\mathbf{k},s}^{ab}$ containing the three-point correlation functions

[TABLE]

which are of order $1/Z^{2}$. The evolution equations for these correlators (12)-(13) can be derived from equation (4) and read after the rotation with $O_{X}^{a}(\mathbf{k})$ (see the Appendix)

[TABLE]

The source terms $S_{\mathbf{p}\mathbf{q},\bar{s}ss}^{G,Xab}$, $S^{I,ab}_{\mathbf{p}\mathbf{q},\bar{s}s}$, and $S^{H,ab}_{\mathbf{p}\mathbf{q},\bar{s}s}$ 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

[TABLE]

Finally, their evolution equation can be derived from Eq. (5). After a rotation with $O_{X}^{a}(\mathbf{k})$, we find (see the Appendix)

[TABLE]

where the source term $S_{\mathbf{p}\mathbf{q}\mathbf{k},\bar{s}\bar{s}ss}^{abcd}$ contains three-point and two-point correlations as well as terms of higher order in $1/Z$, 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 $i\partial_{t}C=\Omega\,C+S$ and thus have formally the solution

[TABLE]

The source terms $S$ containing the distribution functions are slowly varying, with their rate of change being suppressed by $1/Z$ or even more, in comparison with the rapid oscillations $\Omega=\,{\cal O}(Z^{0})$. Hence we may approximate $S(t^{\prime})\approx S(t)$ in the above integral (20) which gives

[TABLE]

with the infinitesimal shift $\epsilon>0$ 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

[TABLE]

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 $M^{abcd}_{\mathbf{p+q,p,k-q,k},\bar{s}\bar{s}ss}$ describe the scattering cross sections for the various processes. For example, $M^{++++}_{\mathbf{p+q,p,k-q,k},\bar{s}\bar{s}ss}$ corresponds to the collision of two quasi-particles with initial momenta $\mathbf{k}$ and $\mathbf{p}$, which are scattered to the final momenta $\mathbf{k-q}$ and $\mathbf{p+q}$, thus satisfying momentum conservation (with the momentum transfer $\mathbf{q}$). 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 $1+f^{d}_{\mathbf{k},s}$ instead of $1-f^{d}_{\mathbf{k},s}$.) Related to this fermionic nature is the particle-hole duality where the distribution function $f^{+}_{\mathbf{k},s}$ describing quasi-particles is mapped to the distribution function $1-f^{-}_{\mathbf{k},s}$ of the holes. Thus, in addition to $2\to 2$ processes such as the collision between two quasi-particles $M^{++++}_{\mathbf{p+q,p,k-q,k},\bar{s}\bar{s}ss}$ or two holes $M^{----}_{\mathbf{p+q,p,k-q,k},\bar{s}\bar{s}ss}$ or a quasi-particle with a hole $M^{--++}_{\mathbf{p+q,p,k-q,k},\bar{s}\bar{s}ss}$, the above equation (22) does in principle also contain $1\to 3$ processes: E.g., $M^{+-++}_{\mathbf{p+q,p,k-q,k},\bar{s}\bar{s}ss}$ 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 $U\gg J$, where such processes are forbidden by energy conservation: The initial particle energy $E^{+}_{\mathbf{k}}\approx U-J_{\mathbf{k}}/2$ is not large enough to create a final state with an energy of nearly $2U$.

As the final analogy to the weakly interacting case, we note that only quasi-particles (or holes) of opposite spins $s$ and $\bar{s}$ 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 $U\hat{n}_{\mu}^{\uparrow}\hat{n}_{\mu}^{\downarrow}$, 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 $M^{abcd}_{\mathbf{p+q,p,k-q,k},\bar{s}\bar{s}ss}$ acquire a non-trivial momentum dependence. To illustrate this, let us consider the limit of strong interactions $U\gg J$. In this limit, the Boltzmann equation (22) describing collisions of two quasi-particles simplifies to

[TABLE]

For the collision of two holes, the equation has the same form after replacing all the $f^{+}$ with $f^{-}$. 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 $U^{2}$. Here, we find that the interaction $U$ does not occur in the Boltzmann equation (23) at all, where the scattering cross section reads $(J_{\mathbf{k}}+J_{\mathbf{p}})^{2}$ and is thus depends on the momenta $\mathbf{k}$ and $\mathbf{p}$ 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 $U$. 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 $J^{2}$ 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

[TABLE]

which assumes its minimum value $\Delta E_{\mathbf{k}}^{\rm min}=U$ at those points where $J_{\mathbf{k}}$ vanishes. Now, a weak enough pump laser with a frequency sufficiently below the gap would predominantly create excitations near those minimum-energy wave-numbers $\mathbf{k}$ where $J_{\mathbf{k}}=0$. On the other hand, for these quasi-particle excitations, the scattering cross sections $(J_{\mathbf{k}}+J_{\mathbf{p}})^{2}$ 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 $\hat{\rho}_{\mu}$. This determines the double occupancy in Eq. (6) via

[TABLE]

However, this small double occupancy $\mathfrak{D}=\,{\cal O}(1/Z)$ 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 $H$-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 $\langle\hat{\sigma}_{\mu}^{x}\hat{\sigma}_{\nu}^{x}\rangle^{\rm corr}$ with $\hat{\sigma}_{\mu}^{x}=\hat{c}_{\mu,\uparrow}\hat{c}_{\mu,\downarrow}^{\dagger}/2+\rm h.c.$), 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

[TABLE]

where $s$ and $s^{\prime}$ are spin indices and $V_{\mu\nu}^{ss^{\prime}}$ 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

[TABLE]

from which one can obtain the equation of motion of the fermion distribution function $n_{\mathbf{k},s}=\langle\hat{c}_{\mathbf{k},s}^{\dagger}\hat{c}_{\mathbf{k},s}\rangle$, i.e.,

[TABLE]

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

[TABLE]

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

[TABLE]

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 $Z$.

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 $U\ll J$ and insulating for $U\gg J$. We demand that the initial state has $\sigma_{z}$-symmetry, such that the density matrix commutes with $\sum_{\mu}(\hat{n}_{\mu,\uparrow}-\hat{n}_{\mu,\downarrow})$ for all times.

B.1 Operator equations

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

[TABLE]

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

[TABLE]

and

[TABLE]

with $U^{0}=0$ and $U^{1}=U$. The operator ${\hat{N}_{\mu,s}}^{X}$ evolves according to

[TABLE]

the spin-flip operator satisfies the equation

[TABLE]

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

[TABLE]

and

[TABLE]

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 $\sigma_{z}$-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,

[TABLE]

The zeroth order equation of the hierarchical expansion (2) determines the double-occupancy $\langle{\hat{N}_{\mu,s}}^{1}{\hat{N}_{\mu,\bar{s}}}^{1}\rangle=\langle{\hat{N}_{\mu,s}}^{0}{\hat{N}_{\mu,\bar{s}}}^{0}\rangle=~{}\mathfrak{D}$, i.e.

[TABLE]

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

[TABLE]

The evolution equation (B.2) involves terms of order $\mathcal{O}(1/Z)$ which determine the free dynamics of the quasi-particles. In this order, each mode evolves independently. The three-point correlations of order $\mathcal{O}(1/Z^{2})$ couple different modes with each other and are crucial in the derivation of the Boltzmann dynamics, see below. In order to represent equation (B.2) momentum space, we define the Fourier components of the two-point correlation function and the various three-point correlation function to be (cf. equations (8) and (12-13))

[TABLE]

[TABLE]

[TABLE]

[TABLE]

With these definitions we find from equation (B.2) (cf. equation (V))

[TABLE]

with a source term determining the free quasi-particle dynamics,

[TABLE]

and a source term of order $\mathcal{O}(1/Z^{2})$ which contains the three-point correlators,

[TABLE]

We omitted in equation (48) the terms which do not contribute to the Boltzmann dynamics in leading order.

It is useful to employ a two-dimensional orthogonal transformation which transforms from the $X-Y$ basis to the particle-hole basis. A general tensor transforms as

[TABLE]

The orthogonal matrix $O^{a}_{X}(\mathbf{k})$ satisfies the eigenvalue equation

[TABLE]

and has the explicit form

[TABLE]

with

[TABLE]

and

[TABLE]

The excitation energies of quasi-particles and holes are (cf. equation (11))

[TABLE]

With the transformation (49) we can rewrite the equations (B.2) as (cf. equation (10))

[TABLE]

After the rotation into the particle-hole basis we can separate the slow degrees of freedom ($a=b$) from the fast degrees of freedom ($a\neq b$) which are changing on a time-scale $\sim 1/U$. Within Markov approximation (cf. equations (20) and (21)) we find

[TABLE]

The slow dynamics is then determined by the evolution of the diagonal elements,

[TABLE]

since the $1/Z$-contributions of the source term in (56) are vanishing for $a=b$.

The correlation functions and the quasi-particle- and hole-distribution functions (which contains also the on-site contribution of order $\mathcal{O}(1)$) are related by the algebraic relation

[TABLE]

The time-evolution for a negligible change of the double occupancy, $\partial_{t}\mathfrak{D}\approx 0$, is then given by

[TABLE]

The hierarchical method relies on a separation of expectation values into correlated and uncorrelated parts. Since we want to express our final result in terms of quasi-particle and hole distribution functions, we need the inversion of the relation (59). It can be checked that up to first order $\mathcal{O}(1/Z)$ we have

[TABLE]

B.3 Boltzmann part of the three-point correlation functions

The second order of the hierarchical expansion (cf. equation (4)) determines the evolution of the three-point correlation functions (43), (44) and (45). Since we are primarily interested in correlations among four lattice sites, we shall omit here the explicit form of the source terms which contain only two- or three-point correlation functions. Some of the equations below end therefore with “…”.

The three-point correlations (43) are the source terms for particle-number correlations. For them we find

[TABLE]

with

[TABLE]

Taking the Fourier transform, switching to the particle-hole basis and integrating within Markov approximation gives (cf. equation (15))

[TABLE]

where we introduced the Fourier components of the four-point correlations (cf. equation (18)),

[TABLE]

The correlation functions (44) are the source of spin-flip correlations and obey the differential equation

[TABLE]

with

[TABLE]

Again, after Fourier transformation and switching to the particle-hole basis, we find within Markov approximation (cf. equation (16))

[TABLE]

Finally, the correlation functions (45) generate the doublon-holon correlations and evolve according to

[TABLE]

with

[TABLE]

which leads to (cf. equation (17))

[TABLE]

All these three-point correlators determine the evolution of the particle- and hole-distribution functions (59). From (60) together with (64), (68) and (71) we find

[TABLE]

B.4 Three-point correlation functions up to $1/Z^{2}$

In the previous section we omitted the $1/Z^{2}$-contribution of the three-point correlation functions since we focused onto the Boltzmann-part which is of order $1/Z^{3}$. As will be shown below, the computation of the Fourier components $J^{abcd}_{\mathbf{q}_{1},\mathbf{q}_{2},\mathbf{q}_{3},\bar{s}\bar{s}ss}$ up to order $1/Z^{3}$ requires the knowledge of $G^{YXab,1/Z^{2}}_{\mathbf{p}_{1},\mathbf{p}_{1},\bar{s}sss}$, $I_{\mathbf{p}_{1},\mathbf{p}_{2},\bar{s}s}^{ab,1/Z^{2}}$ and $H_{\mathbf{p}_{1},\mathbf{p}_{2},\bar{s}s}^{ab,1/Z^{2}}$.

B.4.1 Three-point correlators $G^{YXab,1/Z^{2}}_{\mathbf{p}_{1},\mathbf{p}_{1},\bar{s}sss}$

We begin with the differential equation for the three-point correlations

[TABLE]

The source term reads

[TABLE]

We neglect the particle-number correlations which are of $\mathcal{O}(1/Z^{2})$ and transform the Fourier coefficients in the particle-hole basis. We find the symmetric and anti-symmetric combinations

[TABLE]

and

[TABLE]

with

[TABLE]

and

[TABLE]

B.4.2 Three-point correlators $I_{\mathbf{p}_{1},\mathbf{p}_{2},\bar{s}s}^{ab,1/Z^{2}}$

The inhomogeneity of order $1/Z^{2}$ in (B.3) reads

[TABLE]

The last four lines of equation (B.4.2) are of order $1/Z^{3}$ (the two-site correlations of the spin-flip operators are of order $1/Z^{2}$) and will be neglected in the following. Within this approximation we arrive at

[TABLE]

with

[TABLE]

B.4.3 Three-point correlators $H_{\mathbf{p}_{1},\mathbf{p}_{2},\bar{s}s}^{ab,1/Z^{2}}$

The term of order $1/Z^{2}$ which was omitted in equation (B.3) reads

[TABLE]

Again, the last four lines of equation (B.4.3) are of order $1/Z^{3}$ and will be neglected in the following. After Fourier transform we obtain within the Markov approximation

[TABLE]

with

[TABLE]

B.5 Four-point correlation functions up to $1/Z^{3}$

The differential equation of the four-point correlators originates from the third order of the hierarchical expansion (5) and is given by

[TABLE]

with

[TABLE]

At half filling we find after the Fourier transform in Markov approximation

[TABLE]

with the source term

[TABLE]

B.6 Boltzmann equations

From equations (B.3) and (B.5) we find after some tedious algebra the time-evolution of the distribution functions $f^{d}_{\mathbf{k},s}$, (cf. equation (22)),

[TABLE]

with

[TABLE]

The time-evolution of the double-occupancy is determined by (40). Within the Markov approximation, its Boltzmann-time-evolution reads (cf. equation (25))

[TABLE]

which is of order $\mathcal{O}(1/Z^{4})$ and becomes negligible for $J\ll U$.

B.7 Weak interactions

In equations (54) and (51), the rotation matrix was chosen such that the particle-hole excitation energy is always positive, $E^{+}_{\mathbf{k}}-E^{-}_{\mathbf{k}}>0$. This choice is useful in the limit of strong interactions, see below. However, in the limit of weak interactions, $U/J\ll 1$, the Hubbard bands are overlapping and the system is in a metallic state where the notion of quasi-particles and holes looses its meaning. In the weak-coupling limit, the calculation is simplified considerably if the rotation matrix orders the eigenvalues such that

[TABLE]

Note that we never used the explicit form of the rotation matrix in the above derivation of the Boltzmann equation, therefore we have some freedom as long as the eigenvalue equation (50) is satisfied. Equation (92) corresponds to the choice

[TABLE]

For $U/J\ll 1$, the dominating channel is $a=b=c=d=-$. The remaining matrix elements determine the dynamics of slower collisions with energies $\sim U^{2}/J$ or $\sim U$. Using the energy conserving delta-distribution for the dominating channel, we find from (B.6)

[TABLE]

The evolution equation (B.6) simplifies to

[TABLE]

In this limit the distribution function reads

[TABLE]

and we find

[TABLE]

which is the standard expression of the Boltzmann kinetic equations in the weak coupling limit. It coincides with the perturbative result (A) for $V_{\mathbf{q}}^{ss}=0$ and $V_{\mathbf{q}}^{s\bar{s}}=U$.

B.8 Strong interactions

In the limit of strong interactions $J/U\ll 1$, we choose the rotation matrix $O^{a}_{X}(\mathbf{k})$ such that $E^{+}_{\mathbf{k}}-E^{-}_{\mathbf{k}}>0$. From (51) we find then $O^{a}_{X}(\mathbf{k})\approx\delta^{a}_{X}$. The four-point correlator (B.6) simplifies to

[TABLE]

From (B.6) follows then the evolution equation of the hole modes

[TABLE]

and for the particle modes (cf. equation (23))

[TABLE]

Note that in the strong-coupling limit, the quasi-particle and hole distribution functions are related to the correlation functions via

[TABLE]