Probing thermalization in quenched integrable and nonintegrable Fermi-Hubbard models

Using numerically exact methods we examine the Fermi-Hubbard model on arbitrary cluster topology. We focus on the question of which systems eventually equilibrate or even thermalize after an interaction quench when initially prepared in a state highly entangled between system and bath. We find that constants of motion in integrable clusters prevent equilibration to the thermal state. We discuss the size of fluctuations during equilibration and thermalization and the influence of integrability. The influence of real-space topology and in particular of infinite-range graphs on equilibration and thermalization is studied.

Figure 1

Comparison of the results for the thermal expectation value (5) of the double occupancy $d_i$ stemming from exact diagonalization (solid lines) and TPQ states (circles) for a one-dimensional chain with periodic boundary conditions at $U_i=U=3J=3J_{ij}$. All TPQ results are calculated using a Krylov space dimension of $s=10$ except for the open circles in the $N=10$ case ($s=50$).

Figure 2

Thermal expectation value of the double occupancy in a half-filled $N=10$ Hubbard chain at $U=3J$ calculated by KPM using $k$ moments ${\mu }_n$. The number of moments increases from top to bottom. A higher number of moments improves the accuracy, but it does not change the overall behavior for $\beta J>1$.

Figure 3

Comparison of results from KPM ($k=900$) and TPQ states (Krylov space dimension $s=50$) with results from exact diagonalization (ED). KPM results show notable deviations for increasing $\beta$ while TPQ does not. Results for a half-filled Hubbard chain with periodic boundary conditions of $N=10$ sites at $U=3J$.

Figure 4

Time evolution of the double occupancy on the integrable periodic Hubbard chain after interaction quench $U=3J$ and on the nonintegrable (generic) cluster of $N=12$ sites at half filling. For the nonintegrable case, the time evolution of cluster (l) of Appendix pp1 and its site $i=10$ is depicted and a 1% randomization around the average $U\approx 3J$ is chosen. Solid lines denote the average values (27), and dashed lines denote the average plus and minus the standard deviation ${\sigma }_i$, both calculated at $\tau =0.6$.

Figure 5

Averages of the double occupancy of cluster (l) at the sites $i$ determined according to (27). A tendency of the dynamics to converge towards an essentially constant value around $\tau \approx 0.5$ to $\approx 0.6$ is discernible for all sites.

Figure 6

Global standard deviations $\sigma$ as derived from (30) of the double occupancies $d_i\left(t\right)$ fluctuating around their individual average values ${\overline{d}}_i$. Results for both integrable (PBC) and nonintegrable (generic) clusters are shown. Three different least-square fits to the numerical data are displayed using either $\sigma =a_0+b_0/N$ (upper panel), $log\left(\sigma \right)=log\left(a_1\right)-b_{1}N$ (lower panel, solid lines), or $\sigma =a_{2}exp(-b_{2}N)$ (lower panel, dashed lines). The latter two fits seem to be the same, but this is not the case because the condition of least squares depends on the functional form and leads to differing weights and thus to differing optimum sets $(a_i,b_i)$.

Figure 7

Comparison of the actual global standard deviations $\sigma$ (filled symbols) to the upper bounds (open symbols) given by Eq. (31) on logarithmic scale. The same fits and parameters as in the lower panel of Fig. 6 are shown; they appear here as straight lines.

Figure 8

Infinite-range clusters for $N=4$ and 8 denoting clusters with the maximum number of hoppings possible, also called complete graphs. For each of the $N$ sites the coordination number is $z=N-1$ leading to a total of $K=1/2N(N-1)$ hopping links.

Figure 9

Global standard deviations $\sigma$ of integrable chains with coordination number $z=2$ and of infinite-range clusters $G_c$ with $z=N-1$ and a 1% randomization. Results are to be compared with Fig. 6. The amount of fluctuations depends on the number of bonds and decreases upon increasing coordination number so that the infinite-range clusters display only small fluctuations in the limit $N\to \infty$ relative to the fluctuations in the PBC clusters. The available upper bounds (31) are inserted using open symbols. Dashed and solid lines are fits (cf. Fig. 6).

Figure 10

Global deviation of the time averages ${\overline{d}}_i$ from the thermal predictions ${\langle d_i\rangle }_{\mathrm{th}}$ at the effective temperature for $U=3J$. This deviation ${\mathrm{\Delta }}_{\mathrm{therm}}$ is shown in dependence on the inverse cluster size $N$. For the generic clusters the shown values are averaged over various clusters of the same size and the error bar indicates the spread within this set of clusters. The lines represent linear regressions to the data.

Figure 11

Results for $U=1J$ showing the global standard deviation $\sigma$ as derived from (30) of the double occupancies $d_i\left(t\right)$ fluctuating around their average values ${\overline{d}}_i$ after interaction quenches. In accordance with Fig. 6 fluctuations are becoming exponentially smaller with increasing cluster size $N$. Upper bounds (31) are computed by complete exact diagonalization and shown using open symbols. The bounds are indeed well above the actual data, but they are clearly not tight. Solid and dashed lines denote fits (cf.  Fig. 6).

Figure 12

 Results for an interaction quench of the strength $U=1J$ showing the global deviation ${\mathrm{\Delta }}_{\mathrm{therm}}$ of the time averages ${\overline{d}}_i$ from the thermal predictions ${\langle d_i\rangle }_{\mathrm{th}}$ at the effective temperature of the quench. Integrable (PBC) and nonintegrable (generic) clusters are shown. A pronounced spread between the different generic clusters can be noticed. Since the system is only weakly quenched no clear tendency of thermalization or the absence thereof is visible in the data for both the PBC and the generic systems.

Figure 13

Global standard deviations $\sigma$ for $U=6J$ as derived from (30). The tendency of the fluctuations to decrease exponentially with increasing system size $N$ is obvious. Open symbols show upper bounds (31).

Figure 14

Global deviation of the time averages ${\overline{d}}_i$ from the thermal predictions ${\langle d_i\rangle }_{\mathrm{th}}$ at the effective temperature for $U=6J$. The results agree qualitatively with the ones shown in Fig. 10 where a detailed analysis can be found.