Explicit construction of quasiconserved local operator of translationally invariant nonintegrable quantum spin chain in prethermalization

We numerically construct translationally invariant quasiconserved operators with maximum range $M$, which best commute with a nonintegrable quantum spin chain Hamiltonian, up to $M=12$. In the large coupling limit, we find that the residual norm of the commutator of the quasiconserved operator decays exponentially with its maximum range $M$ at small $M$, and turns into a slower decay at larger $M$. This quasiconserved operator can be understood as a dressed total “spin-$z$” operator, by comparing with the perturbative Schrieffer-Wolff construction developed to high order reaching essentially the same maximum range. We also examine the operator inverse participation ratio of the operator, which suggests its localization in the operator Hilbert space. The operator also shows an almost exponentially decaying profile at short distance, while the long-distance behavior is not clear due to limitations of our numerical calculation. Further dynamical simulation confirms that the prethermalization-equilibrated values are described by a generalized Gibbs ensemble that includes such quasiconserved operator.

Figure 1

Behavior of the squared residual norm ${\lambda }_0\left(M\right)$ (in units of $J^2$) vs maximum range $M$ on (a) log-log plot and (b) semilogarithmic plot, for model parameters $J=1.0$, $h=1.5$, and varying $g$. For small $g$, ${\lambda }_0\left(M\right)$ decays as power-law in $M$. For large $g$, it first decays exponentially as one increases $M$, and then turns into a slower trend at larger $M$. (b) shows additional data from the Schrieffer-Wolff construction of quasiconserved quantity (see Sec. 4 for details), which can be viewed as a variational bound. The residual norm ${\mu }_n$ from the SW construction of order $n$, which corresponds to $M=n+1$ maximum range, shows a classic asymptotic expansion behavior for the smaller $g$ values, where it starts to increase at large order. While this behavior is not manifest yet for the larger $g$ values, from the observed trends we suspect that ${\mu }_n$ will also start to increase eventually beyond some order.

Figure 2

Behavior of the squared residual norm $\lambda \left(M\right)$ for the first five slowest operators in the “TePe” and “ToPo” sectors. (a) For $g=1.0$, the slowest operator in the “TePe” sector shows similar dependence on $M$ as the other nearby slow operators; no particularly slow degrees of freedom exist in this case. On the other hand, in (b) for $g=3.0$ and (c) for $g=5.0$, the slowest operator has exponential dependence on $M$ up to some range, while the other operators decrease more slowly throughout, which suggests that the slowest operator has parametrically more slow dynamics compared to other degrees of freedom.

Figure 3

Comparison between the residual Frobenius norm and operator norm measures of the slowest operator ${\mathcal{Q}}_0^{\left(M\right)}$; the operator is obtained from the minimization of the residual Frobenius norm as described in Sec. 2. The inverse of $\frac{\parallel [H,{\mathcal{Q}}_0^{\left(M\right)}]{\parallel }_{\text{op}}}{\parallel {\mathcal{Q}}_0^{\left(M\right)}{\parallel }_{\text{op}}}$ gives the thermalization time scale of ${\mathcal{Q}}_0^{\left(M\right)}$. For large coupling, cases $g=3.0$ and $g=5.0$, we find that the numerical values of the residual Frobenius and operator norm measures are close to each other up to some $M$ and then start deviating (see text for some discussion).

Figure 4

(a) The overlap between the full numerical optimization ${\mathcal{Q}}_0^{\left(M\right)}$ with $M=11$ and the perturbative SW construction ${\tilde{I}}^{\left(n\right)}$ with order $n=1$ to 10. (b) One minus the overlap on the log-linear plot. At large $g$, the overlap between the two operators is almost $100\%$, which means that the slowest operator we found is essentially the dressed spin operator coming from the solvable limit $H_0$. On the other hand, for small $g$, the slowest operator does not look like the perturbative SW construction operator anymore. Interestingly, there is apparently a strong change in behavior around $g_c\approx 2$; however, we do not know if there is a true transition.

Figure 5

Operator inverse participation ratio of the slowest operator vs maximum range $M$ for different $g$. For large $g\gtrsim 2$, the OIPR appears to converge to a finite value, which suggests its locality in the operator space. On the other hand, in the ergodic regime, $g\lesssim 2$, the OIPR does not converge and instead grows strongly with $M$ (the behavior on the linear-log plot suggests exponential growth).

Figure 6

(a) The weight $W_r$ of range-$r$ operators contained in ${\mathcal{Q}}_0^{\left(M\right)}$ with maximum range $M=12$ for various $g$. For large $g\gtrsim 2$, the weight $W_r$ decays exponentially at short distance $r$. The decay length grows as $g$ decreases. For small $g\lesssim 2$, the decay of $W_r$ is naively better described by a Gaussian, with the curves almost independent of $g$. (b)–(f) The weight $W_r$ of range-$r$ operators in ${\mathcal{Q}}_0^{\left(M\right)}$ when varying $M$ from $M=6$ to $M=12$ for fixed $g$ indicated in each panel. For large $g$, the exponentially decaying part at short distances is essentially converged in $M$; however, the long-distance behavior is not clear. For small $g$, the weight distribution is pushed to larger $r$ and shows significantly slower decay as a function of $r$ when one increases $M$; this suggests that these operators are not normalizable in the large $M$ limit.

Figure 7

TEBD simulations with bond dimensions $\chi =256$ and $\chi =512$ of the evolution of various “magnetizations” $\langle M_{x,y,z}\rangle$ upon quench from the initial state $|Y+\rangle$. The Hamiltonian is given by Eq. (1) with parameters $J=1.0$, $h=1.5$, and different $g$ indicated in each panel. (a) Evolution of the magnetizations for $g=1$. The magnetizations appear to approach the thermal value ${\langle O\rangle }_{\text{th}}=0$ expected for any traceless observable $O$. (b) Evolution of the magnetizations for $g=3$. The magnetizations are approaching values described by the generalized Gibbs ensemble that includes also the quasiconserved operator (see text for details); the expected prethermalized values are marked with subscript “pth.” Insets in both panels show truncation error $1-\langle \psi \left(t\right)\left|\psi \right(t)\rangle$ of the matrix-product states. We set the cutoff for the $\chi =256$ simulation as ${\mathfrak{s}}_0={10}^{-6}$, while for the $\chi =512$ simulation the cutoff is ${\mathfrak{s}}_0={10}^{-8}$.

Figure 8

Numerical calculations of the iterative bounds on $\parallel V_m{\parallel }_1$: (a) ${\tilde{\mu }}_m$ generated by Eq. (D1) and (b) ${\mu }_m$ generated by Eq. (B7), for different SW order $n$. For convenience, the one-norms of $\parallel H_0{\parallel }_1$ and ${\parallel T\parallel }_1$ are taken to be one, which does not affect the functional dependence of the convergence radius ${\rho }_n$ on $n$. The curve $m=n+1$ in (b) denotes the bound on the infinite-SW $\parallel V_m{\parallel }_1$ since $V_m$ does not depend on $n$ once $n\ge m-1$. (Insets) The inverse convergence radius ${\rho }_n^{-1}$ as a function of $n$. By assuming ${\tilde{\mu }}_m=A_n{\left({\rho }_n\right)}^{-m}$, or $ln\left({\tilde{\mu }}_m\right)=ln{A}_n-mln\left({\rho }_n\right)$ for $m>n$, we can extract $ln(1/{\rho }_n)$ from the slope of $ln\left(\widetilde{{\mu }_m}\right)$ vs $m$ and plot ${\rho }_n^{-1}$ in the inset. For ${\rho }_n$ extracted from ${\tilde{\mu }}_m$, we suspect ${\rho }_n^{-1}\sim n$; while for ${\mu }_m$, we observe ${\rho }_n^{-1}\sim n^2$ as expected.

Figure 9

Actual Frobenius norms and one-norms (denoted by “F” and “one” respectively) of operators $V_m$ that appear in the SW transformation, to be compared with bounds in Fig. 8. The model is defined using $H_0$ in Eq. (A1) with $\mathrm{\Gamma }=1$ and $T$ in Eqs. (A2, A3, A4, A5, A6) with $t_1=-12/121$, $t_2=16/121$, $u_0=9/121$, and $w_0=t_2$. The parameters are chosen such that $\parallel H_0{\parallel }_1={\parallel T\parallel }_1=1$ and the ratios among $t_1$, $t_2$, $w_0$, $u_0$ corresponding to the case with $h=1.5$ and $g=2.0$ in the main text. Note that the actual $\parallel V_m{\parallel }_{\text{F}}$ and $\parallel V_m{\parallel }_1$ are still decreasing for the accessible $m$, in stark contrast with the bounds that show very fast increase (at least $m^m$) starting already at $m=1$.