Thermalization and light cones in a model with weak integrability breaking

We employ equation-of-motion techniques to study the nonequilibrium dynamics in a lattice model of weakly interacting spinless fermions. Our model provides a simple setting for analyzing the effects of weak integrability-breaking perturbations on the time evolution after a quantum quench. We establish the accuracy of the method by comparing results at short and intermediate times to time-dependent density matrix renormalization group computations. For sufficiently weak integrability-breaking interactions we always observe prethermalization plateaus, where local observables relax to nonthermal values at intermediate time scales. At later times a crossover towards thermal behavior sets in. We determine the associated time scale, which depends on the initial state, the band structure of the noninteracting theory, and the strength of the integrability-breaking perturbation. Our method allows us to analyze in some detail the spreading of correlations and in particular the structure of the associated light cones in our model. We find that the interior and exterior of the light cone are separated by an intermediate region, the temporal width of which appears to scale with a universal power law $t^{1/3}$.

Figure 1

Green's function $\mathcal{G}(L/2,L/2+1;t)$ on an $L=256$ site system for a quench where the system is prepared in the ground state of $H(0,0.8,0)$ and time-evolved with $H(0,0.4,0.4)$. EOM results (red line) are in excellent agreement with t-DMRG computations [67] (triangles); the green dotted line shows the result of the “integrable quench” [i.e., time evolution generated by $H(0,0.4,0)$; cf. Sec. 3a1]. Inset: Behavior on a larger time interval.

Figure 2

Green's function $\mathcal{G}(L/2,L/2+1;t)$ for a system of size $L=320$ prepared in state with density matrix ${\rho }_0(2,0)$ (14) and time-evolved with $H(0,0.1,0.4)$. The insets show the position of the ED thermal values.

Figure 3

Green's function $\mathcal{G}(L/2,L/2-1;t)$ for a system of size $L=320$ prepared in state with density matrix ${\rho }_0(2,0)$ (14) and time-evolved with $H(J_2,0.1,0.4)$. Dashed lines indicate the thermal values computed by ED of $L=16$ sites. Different colors correspond to different values of $J_2=0,0.25,0.32,0.375,0.425,0.5,0.55,0.6$ (top to bottom).

Figure 4

Green's function $\mathcal{G}(L/2,L/2-1;t)$ for a system of size $L=320$ prepared in the density matrix ${\rho }_0(2,0)$ (14) and time-evolved with $H(0.375,0.1,U)$ for several values of $U$. Solid lines show results from the second-order EOM (34) for $U=0.1,0.2,0.4$ (top to bottom). For $U=0.1$ we compare with the first-order EOM (dashed).

Figure 5

Dispersion relations for the two bands of Bogoliubov fermions in the noninteracting model with $J_1=1,\delta =0.1$ and $J_2=0$ (violet, dotted), $J_2=0.25$ (blue, dashed), $J_2=0.375$ (red, solid), $J_2=0.5$ (green, dot-dashed). Increasing $J_2$ leads to additional crossings at a fixed energy for ${\epsilon }_{+}\left(k\right)$. For $J_2>0.6$ additional crossings at fixed energy appear as the ranges of ${\epsilon }_{+}\left(k\right)$ and ${\epsilon }_{-}\left(k\right)$ overlap.

Figure 6

Green's function $\mathcal{G}(L/2,L/2\pm 1;t)$ for a system of size $L=320$ prepared in state with density matrix ${\rho }_0(2,0)$ and time-evolved with $H(J_2,0.1,0.4)$. Expected steady state thermal values are shown as dotted lines, while the black dashed lines are exponential fits to Eq. (41). Data in the lower panel have been previously reported in Ref. [75].

Figure 7

Real and imaginary (inset) parts of the Green's function $\mathcal{G}(L/2,L/2+2;t)$ for a system of size $L=320$ prepared in state with density matrix ${\rho }_0(2,0)$ and time-evolved with $H(J_2,0.1,0.4)$. Expected steady state thermal values are shown by dotted lines while the black dashed lines are exponential fits to (41). These data have been previously reported in Ref. [75].

Figure 8

Same as Fig. 7 for $\mathcal{G}(L/2,L/2+4;t)$.

Figure 9

Inverse relaxation times ${\tau }_{L/2L/2+j}^{-1}\left(J_2\right)$ with $j=1,-1,2$ (top to bottom) characterizing the late-time behavior of the Green's function [cf. Eq. (41)]. The system is prepared in the initial density matrix ${\rho }_0(2,0)$ (14) and time-evolved under the Hamiltonian $H(J_2,0.1,0.4)$. The system size is $L=320$ (cf. Fig. 3). Errors are estimated by varying the initial time at which the exponential fit is applied.

Figure 10

$\mathcal{G}(L/2,L/2-1;t)$ for a system with Hamiltonian $H(0.5,0.1,0.4)$ and sizes $L=256,320$ initially prepared in the density matrix ${\rho }_0(\beta ,0)$ for four different values of $\beta$. The expected steady state thermal values are indicated by dotted lines and the exponential fit (41) by a dashed line.

Figure 11

Inverse relaxation times ${\tau }_{ii-1}^{-1}(J_2,\delta ,U,\beta )$ (above) and ${\tau }_{ii+2}^{-1}(J_2,\delta ,U,\beta )$ (below) obtained by fitting the data in Fig. 10 with Eq. (41). Error bars are estimated by varying the initial time at which the exponential fit is applied [107].

Figure 12

$\mathcal{G}(L/2,L/2-1;t)$ for a system with Hamiltonian $H(0.5,0.1,0.4)$ and size $L=320$ initially prepared in a density matrix ${\rho }_0(\beta ,{\delta }_i)$, where $\beta$ is chosen such that all initial states have the same energy density. Dashed lines show the exponential fit (41). Gray dotted lines indicate the thermal values computed by ED. The insets show the behavior of $d_j(i;0,0.5)$ [cf. Eq. (43)].

Figure 13

Dynamics of the full Green's function $\mathcal{G}(L/2,L/2+j;t)$ as a function of time. The system is initialized in the density matrix ${\rho }_0(2,0.5)$ and time-evolved with the integrable (free) Hamiltonian $H(0,0,0)$. The light cone's edge is spreading with velocity $v_{\max }=2{\max }_k{\epsilon }_{+}^{\prime }{\left(k\right)|}_{J_2=0;\delta =0}=2J_1$; see Eq. (11). The signal outside of the light cone decays exponentially in time.

Figure 14

Time evolution of $\mathcal{G}(0,111;t)$. The full line is obtained numerically integrating the exact expression (45), while the dashed lines are obtained using the asymptotic expansions (48, 49, 50).

Figure 15

Comparison of the Green's function $\mathcal{G}(0,j;t)$ for an infinite system as a function of time at fixed distance $j=401$ with the approximations (53) and (54). The system is initialized in the density matrix ${\rho }_0(2,0.5)$ and time-evolved with the free Hamiltonian $H(0.5,0,0)$.

Figure 16

Dynamics of the full Green's function $\mathcal{G}(L/2,L/2+j;t)$ as a function of time. The system is initialized in the density matrix ${\rho }_0(2,0.5)$ and time-evolved with the nonintegrable Hamiltonian $H(0.5,0,0.4)$. The signal outside of the light cone decays exponentially in time.

Figure 17

Double-logarithm plot of ${\max }_t\mathcal{G}(L/2,L/2+j;t)$ as a function of the odd separation $j$, for a system of length $L=1920$ $(U=0)$ and $L=320$ $(U\ne 0)$, initialized in the state ${\rho }_0(2,0.5)$ and evolved with $H(0.5,0,U)$ for $U=0,0.1,0.2,0.4$. Different symbols correspond to different values of the interaction and the red dashed line is $\propto j^{-1/3}$.

Figure 18

Airy Function fit of the Green's function $\mathcal{G}(L/2,L/2+j;t)$ as a function of time at fixed distance $j=137$. The system is initialized in the density matrix ${\rho }_0(2,0.5)$ and time-evolved with the nonintegrable Hamiltonian $H(0.5,0,0.4)$.

Figure 19

Double-logarithm plot of the fitting parameter $b$ [cf. Eq. (56)] as a function of $j$. The system, of length $L=1920$ for $U=0$ and $L=320$ for $U\ne 0$, is initialized in the density matrix ${\rho }_0(2,0.5)$ and time-evolved with the nonintegrable Hamiltonian $H(0.5,0,U)$ with $U=0,0.1,0.2,0.4$. Different symbols correspond to different values of $U$ and the red dashed line is $\propto j^{-1/3}$.

Figure 20

The time dependence of the Bogoliubov mode occupation numbers $n_{++}\left(k\right)$ for a system of size $L=320$ initialized in the density matrix $\rho (2,0.5)$ and time-evolved with $H(0.375,0,0.4)$. The different lines are different $k$ modes [we restrict our attention to $0\le k\le \pi /2$, as $n_{\mu \nu }(k,t)=\mu \nu n_{\mu \nu }(\pi -k,t)$, and plot every fourth $k$ mode].

Figure 21

Same as Fig. 20 for $n_{–}\left(k\right)$.

Figure 22

Same as Fig. 20 for $|n_{+-}\left(k\right)|$.

Figure 23

$\mathcal{G}(L/2,L/2+j;t)$ with (upper) $j=1$ and (lower) $j=2$ for a system of size $L=320$ initially prepared in a state with density matrix $\rho (2,0.5)$, and time-evolved with $H(J_2,0,0.4)$. Full lines show results obtained by integrating the EOM (34), and black dashed lines indicate the QBE results.

Figure 24

Green's functions ${\mathcal{F}}_{L/2L/2+1}(\tau ,U)$ (top) and ${\mathcal{F}}_{L/2L/2+2}(\tau ,U)$ (bottom) obtained by numerical solution of the full EOM (34) for a system prepared in the state with density matrix ${\rho }_0(2,0.5)$ and time-evolved with $H(0.5,0,U)$. Results for several values of $U$ are plotted against the rescaled variable $\tau =U^{2}t$.

Figure 25

Double-logarithmic plot of the $U$ dependence of the inverse relaxation times ${\tau }_{L/2L/2+1}^{-1}\left(U\right)$ for a system initialized in the density matrix ${\rho }_0(2,0.5)$ and time-evolved with $H(0.5,0,U)$. Relaxation times are obtained by fitting the EOM results to the form (41). Errors are estimated by varying the initial time at which the exponential fit is applied.

Figure 26

Mode occupation numbers $n_{++}(k,t)$ (top) and $n_{–}(k,t)$ (bottom) for a system of size $L=320$ that is initialized in the thermal state ${\rho }_0(2,0.5)$ and time-evolved with $H(0.375,0,0.4)$. Solid (dotted) lines are obtained by integrating the EOM (QBE). The solid black line is the thermal value found by means of second-order perturbation theory in $U$.

Figure 27

${\mathcal{F}}_{L/2L/2+1}(\tau ,U)=\mathcal{G}(L/2,L/2+1;t,U)$ (top) and ${\mathcal{F}}_{L/2L/2-3}(\tau ,U)=\mathcal{G}(L/2,L/2-3;t,U)$ (bottom). The system is initially prepared in the state ${\rho }_0(2,0)$ and evolved with the Hamiltonian $H(0.5,0.4,U)$, for different values of $U$. The time evolution is obtained by numerical solution of the full EOM (34) and plotted as a function of the rescaled variable $\tau =U^{2}t$.

Figure 28

The $U$ dependence of the exponents ${\tau }_{L/2L/2-1}^{-1}\left(U\right)$ (top) and ${\tau }_{L/2L/2+2}^{-1}{\left(U\right)}_{\mathrm{r}}$ (bottom) obtained from the exponential fit (41), plotted in logarithmic scale [107]. In the bottom panel the error bars are contained within the symbols. Data are presented for the time evolution from the initial state ${\rho }_0(2,0)$ with Hamiltonian $H(0.5,0.1,U)$ for $U=0.15,0.2,0.25,0.3,0.4$ (see also Fig. 27). The errors are estimated by varying the initial time at which the exponential fit is applied.

Figure 29

$\mathcal{G}(L/2,L/2-1;t)$ (top) and $\mathcal{G}(L/2,L/2+2;t)$ (bottom) for two systems with Hamiltonians $H(0.25,0.1,0.2),H(0.5,0.4,0.4)$ and size $L=320$ initially prepared in a thermal state (14) with density matrix ${\rho }_0(2,0)$. The full lines are obtained by integrating the EOM (34) and the black dashed lines are found by means of the QBE.

Figure 30

${\mathcal{G}}_{+-}(L/2,L/2+1;t)$ (top) and ${\mathcal{G}}_{++}(L/2,L/2+1;t)$ (bottom) for a system prepared in the state $\sigma (\infty ,0.2,0)$ and time-evolved with $\mathcal{H}(0.5,0.5,0.1)$. Dotted lines denote the noninteracting GGE.

Figure 31

$\mathcal{S}(1;t)$ [cf. Eq. (83)] evolved by $\mathcal{H}(2,0.005,0.005)$, starting from the ground state of the Majumdar-Ghosh Hamiltonian [131]. The full line is obtained by integration of the first-order EOM (reported in Appendix pp4), while the dashed lines are the prediction of the mean-field approach developed in Ref. [71]. The inset shows a shorter time window where the observable lies on the first quasistationary plateau (black dotted line).

Figure 32

Green's function $\mathcal{G}(L/2,L/2+1;t)$ for a systems of sizes $L=320,384$ that is prepared in the density matrix ${\rho }_0(2,0)$ and time-evolved with $H_3(J_3,0.1,0.4)$ [cf. (A1)]. The expected steady state thermal values computed by ED for $L=16$ are shown by dotted lines. We note that the cases with $J_3\ne 0$ exhibit pronounced finite-size effects in the ED data.

Figure 33

Dispersion relation ${\epsilon }_{+}(k,\delta ,J_3)=-{\epsilon }_{-}(k,\delta ,J_3)$ [see Eq. (A2)] for the two bands of Bogoliubov fermions in the noninteracting model with $J_1=1,\delta =0.1$ and $J_3=0$ (dotted), $J_3=-0.25$ (dashed), $J_3=-0.5$ (solid). Increasing $J_3$ leads to additional crossings at a fixed energy.

Figure 34

The Green's function $\mathcal{G}(L/2,L/2-1;t)$ for a system of size $L=320$ that is prepared in the density matrix ${\rho }_0(2,0)$ and time-evolved with $H(0.5,{\delta }_f,0.4)$ [cf. (1)] with ${\delta }_f=0.1,0.3,0.5,0.7$ (top to bottom). The expected steady state thermal values are shown by dotted lines.

Figure 35

Time evolution of $\mathcal{G}(L/2,L/2\pm 1;t)$ for a system initialized in the density matrix ${\rho }_0(2,0,5)$ (14) and time-evolved with the translationally invariant Hamiltonian $H(0.5,0,0.4)$.

Figure 36

The Green's function $\mathcal{G}(L/2,L/2+1;t)$ with time evolution generated by (blue) $H(0.5,0.1,0.4)$ and (red) $H(0.5,0.1,-0.4)$. The system starts from the thermal state (14) with ${\beta }_i=2,J_2=0,{\delta }_i=0$, and $U=0$.