Strong and weak thermalization of infinite non-integrable quantum systems

Methods

We consider an infinite translationally invariant spin chain with nearest-neighbour interactions of the Ising type, plus a local magnetic field with transverse ($g$) and parallel ($h$) components to the two-body interaction.

With a parallel ($g=0$) or transverse ($h=0$) magnetic field, the model is exactly solvable, but at a different angle, we have a non-integrable model, i.e. the energy density is the only conserved local quantity The non-integrability of the Hamiltonian can be assessed also from the point of view of its spectral statistics izrailev90rep (see Supplementary Material).. We simulate numerically the time evolution of various initial configurations under fixed Hamiltonian parameters, $g=-1.05$ and $h=0.5$, far from any integrability limits. As initial configurations, we choose translationally invariant product states. They are determined by the state of an individual spin,

In particular, the representative states discussed above correspond to parameters $\theta=\frac{\pi}{2}$, $\phi=0$ ($|X+\rangle$), $\theta=\phi=\frac{\pi}{2}$ ($|Y+\rangle$) and $\theta=0$ ($|Z+\rangle$).

Using the recently developed folding method banuls09fold, we compute all the time dependent expectation values of one-, two- and three-body operators, for each initial state, what allows us to reconstruct the whole reduced density matrix for up to three sites. The thermal state $\rho_{th}(\beta)$ with the same energy as the initial state is also calculated numerically with the same method (see Supplementary Material). The distance between the evolved and thermal reduced density operators is then measured by the operator norm of their difference, $d(\rho_{1},\rho_{2})\equiv\|\rho_{1}-\rho_{2}\|_{op}$, which in this case coincides with the maximum eigenvalue in absolute value of the difference $\rho_{1}-\rho_{2}$.

References

I Supplementary material

II The Hamiltonian

The model we consider is an infinite translationally invariant spin chain with an Ising type nearest neighbour interaction, plus a magnetic field

When $g=0$, the Hamiltonian is trivially solvable, while for $h=0$ it reduces to the integrable Ising model with a transverse field, for which the exact solution can be found by fermionization. In any other case, the model is non-integrable.

We fix the initial values of the Hamiltonian parameters for our study to $g=-1.05$ and $h=0.5$, which is not close to either of the integrable situations. In this situation, the energy density is the only conserved quantity. We have additionally analyzed the spectral properties of the Hamiltonian, for a finite system of length $L$, to check the non-integrability in the sense of a spectrum with the characteristics of a random matrix ensemble. As shown in Fig. 5, already for $L=14$ the level spacing distribution evidences the non-integrability of the chosen Hamiltonian also from the point of view of its spectral properties. For comparison, the level spacing distributions in both integrability limits are also shown.

It could be argued that the non thermalization we observe occurs for states which lie close to the edges of the spectrum, as $|Z+\rangle$ and $|X+\rangle$, and at these energies there could be an integrable effective model distasio95, even for $g=-1.05$, $h=0.5$. The spectral properties discussed in Fig. 5 would be dominated by the central part of the spectrum and not reflect the properties at very low or high energies, while the spectrum in an energy interval in these regions should show a very different behavior. To discard the integrability of the system in the interesting cases, we have checked the level statistics of a small energy window around the energy per site of a weak thermalizing state (Fig. 6) and a non-thermalizing one (Fig. 7) for the case of a finite system ($L=14$) which can be exactly solved. In both cases we have found that the level spacing distribution is typical of a non-integrable system, also in these regions of energy.

III The numerical method

The time evolution of an infinite 1D quantum system is simulated numerically within the Matrix Product States aklt88; kluemper91; kluemper92; fannes92fcs; verstraete04dmrg; perez07mps (MPS) formalism, using the new technique introduced in banuls09fold. With this folding method, it is possible to study the out-of-equilibrium dynamics of the system after longer times than with other similar methods.

In this method, any time dependent expectation value $\langle\Psi(t)|O|\Psi(t)\rangle$ can then be expressed as a two dimensional tensor network (Fig. 8(a)), constructed from the a Suzuki-Trotter expansion of the evolution operator trotter59; suzuki90. Each discrete time step corresponds then to a sequence of matrix product operators (MPO) murg08mpo. In contrast to the standard MPS algorithm, in which the evolved state is approximated by an MPS after each time step by means of successive truncations vidal03eff; vidal07infinite, the new algorithm performs the contraction of the tensor network in the transverse direction, i.e. along space. The left and tight semi-infinite lattices can then be effectively substituted by the left and right dominant eigenvectors of the transfer matrix of the evolved state, $\langle L|$ and $|R\rangle$. Before contracting we apply a folding to the network along the time direction (Fig. 8(b)). The folding operation can be understood as performing the contraction $\langle\Psi(t)|O|\Psi(t)\rangle=\langle\Phi|(O|\Psi(t)\rangle\otimes|\bar{\Psi}(t)\rangle)$, where $|\Phi\rangle$ is a product of unnormalized maximally entangled pairs between each site of the chain and its conjugate (see banuls09fold for a detailed discussion of the algorithm). This is equivalent to grouping together tensors that correspond to the same time step in $\Psi$ and its Hermitian conjugate, and achieves a more efficient representation of the entanglement in the transverse direction, which in turn gives access to the simulation of longer evolution times.

The folding technique is appropriate for the simulation of time evolution, but using imaginary time evolution, it is also possible to efficiently compute any local expectation value in a thermal state. We obtain in this way the dependency of the thermal state energy with the inverse temperature $\beta$, $E_{th}(\beta)$ (Fig. 9). From these data we may compute, for each one of the initial states we have studied, the value of $\beta$ corresponding to the thermal state with the same energy. This determines the state to which the initial configuration would be expected to relax, being energy density the only conserved quantity constraining the relaxation. It is then possible to compare the $N$-particle reduced density matrix for the evolved state and for the thermal state at $\beta$, to compute the desired distance between density operators.

To determine the reduced density matrix for $N$ sites, both for the thermal and the evolved state, we need to compute the expectation values of all $N$-term products of Pauli matrices and the identity. Both in the cases of the thermal and the evolved state, the computation of the dominant left and right eigenvectors is common to all such expectation values. We show here the results for $N=3$ sites.

IV Accuracy of the results

Due to the numerical nature of the study, all the results are subject to errors. In order to assess the validity of our conclusions, we discuss here the character and magnitude of these errors using different criteria.

The numerical method has two sources of error. The first one is the Trotter decomposition. The approximation of the evolution operator by a product of exponentials introduces an error, which can be controlled by either reducing the parameter $\delta$, which determines the size of the time step, or by using a Suzuki-Trotter expansion which is exact to a higher order in $\delta$. Both decreasing the time step or increasing the order of the expansion involve longer vectors in the transverse direction, which will potentially worsen the truncation error. It is then convenient to find a trade-off between both factors. In our analysis we use a fourth order Suzuki-Trotter decomposition with time step $\delta=0.1$. We have checked the convergence of our results with this order and time step.

A second, generally less benign source of error is truncation, i.e. the fact of approximating a given vector by a MPS of fixed bond dimension, which constitutes the main source of numerical errors in any MPS algorithm. In our scheme the evolved state is not explicitly truncated, but truncation takes place in the transverse direction, when we approximate the dominant left and right eigenvectors of the transfer matrix by MPS. In a standard MPS algorithm for time evolution, truncation errors are dramatic, in the sense that, when they appear, the results deviate abruptly from the exact ones, and it becomes imposible to extract information from the computed quantities as soon as truncation error sets in schollwoeck06tDMRG. However, in the folding technique, we may extract information from longer time simulations, because errors come in smoothly and, even when some truncation error occurs, the method is still able to provide a qualitative description of the physics, since we expect that our predictions deviate smoothly from the exact values (see Fig. 10 and discussion in banuls09fold).

To bound the numerical errors in the non-integrable model we may compare our results to those from other approaches, check the convergence of the results as we increase the bond dimension, or make use of some external physical criterion to assess the consistency of the computed numbers. We have used all three kinds of tests. First, we have cross-checked the results of our simulations with some large bond dimension simulations using the iTEBD algorithm vidal07infinite, in which contraction is done in the time direction. As shown in Fig. 10, the folding results with $D=120$ are accurate to the longest times we can simulate with iTEBD, $t\simeq 9-10$. Most remarkably, this bond dimension is enough to get a qualitative description (precision $1\%$) of the dynamics to even longer times. Second, to witness the appearance of truncation errors, we have run the simulations with increasing bond dimension.The comparison of our results with highest bond dimension $D=240$ with those for $D=120$ gives us a bound on the error, which we represent on the plots as an error bar.

Finally, as a physical check of the consistency of the results, we test the conservation of energy along the evolution. We study the unitary evolution of a closed system, and the energy per particle must be constant. However, the numerical implementation does not enforce this condition. On the contrary, from the Suzuki-Trotter expansion to the truncation of the dominant eigenvectors, the numerical errors will in general violate this condition. A very large deviation between the time-evolved energy and the initial one would warn us about the validity of the results. As shown in Fig. 11 and 12 for the initial states $|Z+\rangle$ and $|X+\rangle$, with bond dimension $D=240$, the relative error in the energy for the range of times we are analysing (respectively $t\simeq 18$ and $t\simeq 12$) is kept to only a few percent, consistent with the estimated truncation error. For the initial state $|Y+\rangle$, with zero initial energy, we plot instead the expectation value of energy as a function of time (see Fig. 13).

One may think that this deviation of the energy could also introduce an error in the distance we are computing, as the thermal states corresponding to the computed energy density and to the initial one will be different. To bound this error, we have computed the distance between such pair of thermal states, corresponding to the initial state and to the largest value of the energy found in the evolution, to the range of times we are showing. We find that this distance is significantly smaller than the one we observe during the dynamical evolution. In particular, for the $|Z+\rangle$ initial state, the deviation in energy at times $t\approx 10$ reaches a $1\%$, which corresponds to a distance $d(\rho(\beta_{0}),\rho(\beta^{\prime}))\approx 9\times 10^{-3}$, while the observed distance between the thermal state and the evolved one oscillates around $0.15$. For the longest times we show $t\approx 18$, the largest distance grows to a maximum value of $0.06$. For the $|X+\rangle$ initial state, the maximum deviation in energy, at $t\approx 12$, corresponds to a distance $d(\rho(\beta_{0}),\rho(\beta^{\prime}))\approx 11\times 10^{-3}$, while the one we find at this same long time is $0.04$

V Detailed results

Here we compile our results using various initial states and Hamiltonian parameters, to show the survival of the different thermalization regimes over a range of parameters.

Our initial configurations, translationally invariant product states, are specified by the state of a single spin,

which can be represented on the Bloch sphere by the point with coordinates $(\theta,\phi)$. This state has energy per spin $E=-\left(\cos^{2}\theta+g\sin\theta\cos\phi+h\cos\theta\right)$. States with spins polarised in the three orthogonal directions, $|X+\rangle$, $|Y+\rangle$ and $|Z+\rangle$, behave very differently with respect to thermalization regarding the distance between the reduced evolved density matrix and the thermal one. We may additionally check that their dynamics are essentially different by looking at the individual expectation values (see also Fig. 21). For the initial state $|Y+\rangle$, showing strong thermalization, all of the individual expectation values converge fast to the thermal ones. Instead, for the initial state $|Z+\rangle$, we observe irregular oscillations, showing no sign of damping, to the longest times we are able to simulate. For the initial $|X+\rangle$ state, we check that only few expectation values deviate from the thermal average, and are those preventing thermalization of the whole reduced density matrix. In particular, for $N=1$, only $\langle\sigma_{x}\rangle$ is responsible for the lack of thermalization. The behavior for larger reduced density matrices $N=2,$ $3$ is qualitatively similar, although the time it takes for $|Y+\rangle$ to thermalize becomes longer.

By rotating the inital state on the Bloch sphere, we observe a transition from the strong thermalizing $|Y+\rangle$ to the weak one $|Z+\rangle$, and also to the apparently non-thermalizing $|X+\rangle$ (Fig. 15), the latter occurring around $\phi\in[\frac{\pi}{6},\frac{\pi}{4}]$ ($\beta\in[-0.5382,-0.3915]$ ). We have also analysed the transition from weak thermalization $|Z+\rangle$ to non-thermalization $|X+\rangle$ (Fig. 16). In this case we find some intermediate states for which strong thermalization occurs. By looking at the energy, we may infer the corresponding $\beta$ of every initial product state. We discover that all the strong thermalizing states we have found have energies, and thus $\beta$, close to zero.

Instead, for larger $\beta>0$, we observe oscillations and weak thermalization as in the $|Z+\rangle$ initial state. An interesting case is the state of maximum $\beta$ (the product state with minimal energy), which we can identify, by studying the energy landscape over the Bloch sphere, at $\theta\approx 0.43$, $\phi=\pi$. The dynamics of this initial state shows also weak thermalization, with strong oscillations of all the expectation values, but fast thermalization in average (Fig. 17). We may perform a similar analysis on some extra states, with initial energy densities close to $|X+\rangle$ and $|Y+\rangle$, respectively, to check whether they relax to similar thermal states. We find that they seem to show the same regime of thermalization (weak or strong) as the original states (Fig. 17), suggesting this has to do with the initial energy. The relaxation curves however differ, indicating that, even starting with the same $\beta$, the state of the system does not relax to the same state, at least during the long range of times we simulate.

V.2 Varying the Hamiltonian parameters

Finally, we have also studied how the Hamiltonian parameters affect the appearance of the non-thermalizing behaviour, to ensure that this behaviour is not singular to our particular choice.

As a reference, we may compare the behavior of the same initial states under the chosen Hamiltonian and the integrable one, corresponding to $h=0$. In Fig. 18 we compare the dynamics under both models, $h=0$ and $h=0.5$, for and the three most representative cases. We observe that in the integrable case, the $N=3$ reduced density matrix appears to relax fast to a state which is not the thermal one, since the distance converges to a value different from zero (This could be compatible with the generalized thermal ensemble rigol07gge).

We test also other values of the parameters $g$ and $h$ in the non-integrable regime, to study how the strong and weak thermalizing regimes appear. Keeping $h=0.5$ constant, we observe (Fig. 19) that both regimes are present for a large range of values of $g$, but as we decrease $g$, weak thermalization becomes dominant, while for higher values of $g$, the behaviour approaches strong thermalization in most cases.

If we do the same now for the transition to non-thermalizing as observed between $|Y+\rangle$ and $|X+\rangle$, we also observe that both types of behaviour survive over a wide range of values of $g$ (see Fig. 20).

To get a more detailed idea of the differences among the various types of thermalization behavior, we study the individual expectation values for different observables. For clarity, we show the plots only for the $N=1$ operators in Fig. 21, for constant parallel field $h=0.5$ and varying $g=-0.5,$ $-1.05$, $-1.5$, in two of the extreme cases, $|Z+\rangle$ and $|X+\rangle$. We observe that the oscillating behaviour of the initial state $|Z+\rangle$ appears clearly correlated with the value of $g$, and for big values the oscillations of the expectation values are clearly damped. The behaviour of the initial state $|X+\rangle$ is quite different, and in all cases we observe thermalization of some observables while others deviate from the thermal expectation value.