Test of quantum thermalization in the two-dimensional transverse-field Ising model

We study the quantum relaxation of the two-dimensional transverse-field Ising model after global quenches with a real-time variational Monte Carlo method and address the question whether this non-integrable, two-dimensional system thermalizes or not. We consider both interaction quenches in the paramagnetic phase and field quenches in the ferromagnetic phase and compare the time-averaged probability distributions of non-conserved quantities like magnetization and correlation functions to the thermal distributions according to the canonical Gibbs ensemble obtained with quantum Monte Carlo simulations at temperatures defined by the excess energy in the system. We find that the occurrence of thermalization crucially depends on the quench parameters: While after the interaction quenches in the paramagnetic phase thermalization can be observed, our results for the field quenches in the ferromagnetic phase show clear deviations from the thermal system. These deviations increase with the quench strength and become especially clear comparing the shape of the thermal and the time-averaged distributions, the latter ones indicating that the system does not completely lose the memory of its initial state even for strong quenches. We discuss our results with respect to a recently formulated theorem on generalized thermalization in quantum systems.

In recent years the non-equilibrium dynamics of isolated many-body quantum systems, defined by the unitary time evolution starting from a generic non-eigenstate of the Hamiltonian, has gained tremendous interest  . A central point of fundamental importance is the nature of the stationary state of the observables of the system. Here one is especially interested in answering the question whether the system thermalizes, i.e. whether time-averaged observables and their probability distributions can be described by the canonical Gibbs ensemble (CGE) 6,16,[30][31][32][33][34][35][36][37][38][39][40][41] : The Lagrange multiplier β = 1/T is the inverse temperature and is determined under the constraint of maximizing the entropy.
Thermalization is closely linked to the conserved quantities of the system. While in non-integrable systems only the energy is conserved, integrable systems possess additional conserved quantities which avoid thermalization. The canonical Gibbs ensemble thus cannot be applied to integrable systems, but it has been shown that their stationary state can be well described by the generalized Gibbs ensemble (GGE) [2][3][4][7][8][9][10][11][12]14,22,[42][43][44][45][46] , which includes all the conserved quantities of the system, the so-called charges: with Tr (2) H I H I GGE GGE GGE n n n n n n The conserved charges are taken into account by the operators Î n and the Lagrange multipliers β  and λ n are uniquely determined by the initial conditions. However new results show that care has to be taken in the definition of the GGE [47][48][49][50][51] . While in the past only local charges of the system were considered in the GGE, recently the existence of previously unknown quasilocal charges for different field theoretical models 52 as well as for the the spin-1/2 Heisenberg chain 53,54 has been shown. These quasilocal charges have to be included in the GGE, too, to give an adequate description of the stationary state of the systems. An universal mathematical framework for the construction of the GGE including quasilocal charges has been recently formulated along with the conditions under which generalized thermalization should occur in systems of any dimensionality 55 .
In non-integrable systems the energy is the only conserved quantity, but non-integrability is not synonymous to thermalization. Counterexamples are the driven Rabi model 56 or a non-integrable model of hard-core bosons on connected triangular lattices 57 . For the latter one it has been shown that the applicability of the CGE to the stationary state depends on the symmetries in the system: In case of an extensive number of local symmetries the GGE has to be applied rather than the CGE.
Predictions on whether non-integrable systems thermalize or not usually rely on the eigenstate thermalization hypothesis (ETH) [2][3][4][58][59][60][61][62] , which computes matrix elements of observables in energy eigenstates of the system. The ETH is a sufficient but not a necessary condition for thermalization and has been applied to a wide variety of systems. In contrast to this there are only few studies on thermalization in non-integrable systems for system sizes larger than those accessible with exact diagonalization which compute the time evolution of the systems. Among them are studies of the one-dimensional Bose Hubbard model with time-dependent density matrix renormalization group theory (t-DMRG) 6 and of the antiferromagnetic anisotropic Heisenberg chain with a Chebyshev polynomial expansion 63,64 . In higher dimensions there are results for a D-dimensional effective O(N)-Hamiltonian close to dynamical critical points based on a renormalization-group method [65][66][67] .
Here we present for the first time a systematic investigation of the relaxation dynamics of a two-dimensional, non-integrable model going beyond system sizes accessible with exact diagonalization and applicable to large areas of the parameter space. We study the transverse-field Ising model in two dimensions (2D-TFIM) after global interaction quenches in the paramagnetic phase and global field quenches in the ferromagnetic phase. In contrast to the Ising chain in one dimension, the 2D model is non-integrable and cannot be solved analytically. We use a real-time variational Monte Carlo (rt-VMC) method to compute the time evolution of observables like magnetization and correlation functions with high accuracy for long time scales and large system sizes. To answer the question whether the 2D-TFIM thermalizes or not we compare the time-averaged distributions of the observables to their thermal distributions for the system in equilibrium at temperatures determined by the excess energy after the quench. For a system that thermalizes one would expect the asymptotic time-averaged distributions and the thermal distributions to be identical. We discuss our results with respect to the theorem on generalized thermalization 55 . For the interaction quenches in the paramagnetic phase, for which the conditions of the theorem are fulfilled, we indeed observe thermalization. For the field quenches on the other hand we find a continuously increasing degree of non-thermalization with increasing quench strength.

Model and Methods
The model. We study the 2D-TFIM with nearest neighbour interactions on a square lattice of size L × L with periodic boundary conditions (PBC). The system is described by the Hamiltonian with coupling strength J and external transverse field h. The total number of spins in the system is N = L 2 , thus the dimension of the Hilbert space  is 2 N . As basis of  we choose the x-basis, in which the operator σ x R measures the orientation of the spin at site R, while σ z R inverts it. The model is highly symmetric. Its Hamiltonian is invariant under the global  2 spin flip transformation σ σ → −x x R R and σ σ →ẑ z R R generated by the unitary operator σ Σ = Πẑ z R R . Due to the square lattice with PBC the Hamiltonian also shows translation, rotation and reflection symmetries. Their generators can be constructed from the unitary transposition operators These symmetries can also be found in the eigenstates of the Hamiltonian and are conserved under unitary time evolution.
Ferromagnetic long-range order exists in the thermodynamic limit (L → ∞ ) at low temperatures and fields due to spontaneous symmetry breaking of the global spin flip symmetry and is indicated by a non-vanishing ground state magnetization in In the ferromagnetic phase the system has no energy gap in contrast to the paramagnetic phase. For h = 0 the equilibrium phase transition occurs at (T/J) crit ≈ 1.135 68 and for T = 0 at (h/J) crit ≈ 3.044 [69][70][71] . As there is no spontaneous symmetry breaking in the finite system, here due to the global spin flip symmetry irrespective of the values of T/J and h/J. For this reason we use for the finite systems in our numerical studies the modulus of the magnetization as order parameter. In order to get rid of finite size effects it is renormalized according to the expectation value of the modulus of the magnetization in the completely uncorrelated state. In the thermodynamic limit it is Figure 1(a) shows the phase diagram of µ  x CGE for a 16 × 16 system. Already for this system size one observes only small deviations from the results for the system in the thermodynamic limit, whose critical values for T/J and h/J can be obtained from the Binder cumulant applying finite size scaling 72 .
The Hamiltonian of the 1D-TFIM can be diagonalized by a transformation to a system of free fermions 73 . For the 2D-TFIM this is not possible as its Hamiltonian is non-local after the 2D-Jordan-Wigner transformation [74][75][76] , so that there is no canonic transformation to a system of free fermions. For this reason its relaxation process after a quench cannot be described with the semiclassical theory of non-interacting quasiparticles introduced for the 1D-TFIM either 20 .
of the two fully magnetized basis states of the x-basis with the energy . Both the initial state of the interaction quenches and the field quenches are invariant under global spin inversion, i.e. there is no symmetry breaking of the global spin flip symmetry and thus . For the determination of the effective temperature after the quenches the thermal energy (see Fig. 1(b)) of the final Hamiltonian has to be equal to the ground state energy of the system before the quench. This is illustrated in Fig. 2, which shows  E as a function of T/J for different values of h/J for the system with L = 16. The results for T eff for the interaction and the field quenches are shown in Fig. 3(a) and (b) respectively for different system sizes. The energy argument predicts that the end points of the interaction quenches always lie in the paramagnetic phase, i.e. starting from the completely uncorrelated state the system cannot be driven into the ferromagnetic phase switching on a coupling between the spins. There is a minimum of T eff /J for h/J ≈ (h/J) crit . A dependency of the effective temperature on the system size can only be observed for interaction quenches ending close to the phase transition. This is due to the increasing correlation length in the vicinity of the phase transition. Field quenches are predicted to drive the system out of the ferromagnetic phase when the external field is quenched to values larger than ≈ h J h J / ( / ) 1 2 crit . At the phase transition the energy isoline = −  E h J T J ( / ; / ) 1 , which defines the end points of the field quenches, shows a turning point, which becomes more pronounced with increasing system size. Within the ferromagnetic phase the system size has almost no effect on the effective temperature attributed to the quench. Only in the vicinity of the phase transition larger deviations between the effective temperature can be observed for different system sizes. These deviations decrease again in the paramagnetic phase.
The shape of T eff /J as function of h/J can be understood considering the effects of J, h and T onto the order in the system. In the ground state at T = 0 for h/J < (h/J) crit a parallel orientation of the spins in x-direction is energetically preferred, while for h/J > (h/J) crit the orientation of the spins in z-direction parallel to the external transverse field is favourable. The temperature causes fluctuations of the orientation of the spins and thus disturbs the described order. The strength of this effect depends on h/J. For strong couplings or strong fields higher temperatures are necessary to disturb the order of the spins. This can be seen considerng the energy as a function of T/J for fixed h/J in Fig. 2. The energy increases monotonically as a function of T/J with . The high temperature limit follows from the symmetry of the energy eigenvalues with respect to 0. For T = 0 on the other hand the energy  E is close to − 1 for small values of h/J, while for large ratios of h/J it approaches − h J 2 from below. For low temperatures there is a temperature interval in which the energy is almost not affected by the increase of the temperature. The length of this interval depends on h/J. The closer h/J is to (h/J) crit , the less ordered the spins are and thus can be more easily disturbed by the temperature. As for the interaction quenches the expectation value of the energy of the system in thermal equilibrium at T eff has to be equal to − h J 2 , this causes the minimum of T eff /J close to (h/J) crit . Away from (h/J) crit the effective temperature of the system after the interaction quenches increases. Although for large ratios h/J the energy  E at T = 0 is only slightly smaller than the energy − h J 2 in the conditional equation for the effective temperature, due to the strong transverse field high temperatures are necessary to disturb the order of the spins. For interaction quenches with small ratios h/J the same argument holds. In this case the order is caused by the coupling between the spins. For the field quenches the effective temperature is determined by the condition that = −  E 1. Obviously it is T eff /J = 0 for h/J = 0. In the limit h/J → ∞ on the other as the ground state energy of the system at T = 0 is lowered almost linearly with h in this case.  Time evolution. We prepare the system in the ground state of the initial Hamiltonian Ĥ i before the quench, i.e.
i,0 with ketΨ i,0 according to equations (9) and (10) respectively. In general the initial state of the system is not an eigenstate of the final Hamiltonian, so that the system evolves unitarily in time according to the Schrödinger equation In terms of the eigenbasis of Ĥ f the time evolution of the expectation value of an arbitrary operator  reads The initial state of the system is a pure state and its state remains a pure state under unitary time evolution. A thermal system on the other hand is described by a mixed state. For this reason time averages of the observables after the quench have to be compared to the thermal values as will be shown in the following. Considering the time-evolved expectation value in equation (13), one observes that the diagonal part is time-independent, while the non-diagonal contributions consist of harmonic oscillations. Averaging over time the non-diagonal part vanishes for long time intervals if there are no degenerate energy eigenvalues, so that the stationary state is determined only by the diagonal part: The stationary state of the system can thus be described by a mixed state in the so-called diagonal ensemble: As the dimension of the Hilbert space  grows exponentially with the system size, the above computations in the eigenbasis of Ĥ f with exact diagonalization are only possible for small systems. For this reason we use a rt-VMC method to give an accurate description of the time evolution of the 2D-TFIM for larger system sizes.

Real-time variational Monte Carlo.
Rt-VMC was introduced by Carleo et al. for the Bose-Hubbard model 77,78 and has also been successfully applied to lattice bosons and spin systems with long-range interactions 79 as well as to strongly correlated electron systems 80 . Its idea is the existence of a set of variational parameters which are sufficient to describe the physical properties of the system while their number is much smaller than the dimension of the Hilbert space. The equations of motion of the variational parameters are determined minimizing the Euclidian distance  t ( ) between the exact time evolution Ψ  t ( ) exact of the variational state and its variational time evolution Ψ  t ( ) var , which reads in the x-basis x exact v ar 2  A common choice for the variational state is the Jastrow ansatz [77][78][79] , which is well suited to describe the time evolution of the 2D-TFIM after interaction quenches in the paramagnetic phase. For the field quenches in the ferromagnetic phase we introduce a new ansatz, which makes use of the symmetries of the model and the high degree of order in this phase. Both ansatz functions reduce the number of parameters in the wave function of the system from a number growing exponentially with the system size to a number growing algebraically with the number of sites.
Paramagnetic phase. In the paramagnetic phase we use the Jastrow ansatz for the variational function. This ansatz is constructed from the completely uncorrelated state. Correlations are taken into account by the Jastrow factor. For the 2D-TFIM the Jastrow ansatz reads 79 The operators measure the correlations between all spin pairs of the system with distance r normalized by their number N r . The sum over r runs over all independent directions in the lattice, whose number also determines the number of variational parameters α r (t), which is N/8 + 3L/4 for the square lattice with edge length L. For a given r, the average over all dependent directions is taken (see the supplementary material). The state ↑↑ … ↑↑ z is the completely uncorrelated state of the system, which is the exact ground state in the limit h/J → ∞ , i.e. which represents the completely paramagnetic state. Inserting the Jastrow ansatz from equation (18) into equation (17), one gets the following equations of motion for the variational parameters 77,78 : Ferromagnetic phase. Due to its construction from the completely paramagnetic state, the Jastrow ansatz is well suited for the paramagnetic phase, but fails to describe the time evolution of the system after field quenches in the ferromagnetic phase. For example the initial state cannot be represented by it. We thus derive in the Appendix the ansatz , m n , for the field quenches in the ferromagnetic phase. |Ψ m,n 〉 is the normalized symmetric superposition of all basis states with m spin down and n broken bonds, so-called kinks. The ansatz separates the Hilbert space of the system into subspaces  m n , of states with the same magnetization per site n xx The dimension of the subspace (m, n) is N m,n . Transitions between subspaces are induced by spinflips. A spinflip increases or decreases m by one and keeps n untouched or changes it by ± 2 or ± 4. Thus each subspace is linked to up to 10 other subspaces. We denote the total number of transitions between the subspaces (m, n) and (m′ , n′ ) with T m,n;m′,n′ . The T m,n;m′,n′ are symmetric with respect to (m, n) ↔ (m′ , n′ ). The determination of the N m,n and the T m,n;m′,n′ is a pure combinatorics problem. As for the 2D-TFIM there are no closed-form expressions for their values and due to the high dimensionality of the Hilbert space we determine them with rare event sampling (RES) 82 described in the supplementary material. As the N m,n and T m,n;m′,n′ are independent of the quench parameters, they have to be determined only once for each system size. The possible values of n depend on m in a non-trivial way. In case of the 2D-TFIM not just the maximal number of kinks is a function of m like in one dimension, but also their minimal number. The PBC have to be taken into account, too. In leading order the number of possible values of n grows linearly with m. As the number of possible values of m grows linearly with the system size N, the number of variational parameters thus grows in leading order with N 2 . Due to the symmetry of the variational parameters with respect to m ↔ N − m and the conservation of this symmetry under time evolution, we can reduce the number of independent variational parameters using α m,n (t) = α N−m,n (t). For the system sizes we studied the numbers of the variational parameters are listed in Table 1.
The equations of motion of the variational parameters are derived in the Appendix: The sum over m′ and n′ runs over all subspaces that can be reached from any basis state of the subspace (m, n) by flipping one single spin. The equations of motion of the variational parameters are thus a system of coupled linear differential equations of first order with constant (time-independent) coefficients, which are known from RES. As each subspace is linked to only up to 10 other subspaces, the system is sparse. We solve it with a fourth order Runge-Kutta scheme. The initial values of the variational parameters are For t > 0 the α m,n (t) are in general complex.
Observables. As observables we consider the rescaled modulus of the magnetization according to equation (4) and the correlation function between two spins at distance r according to equation (19), which reads for nearest neighbours For the interaction quenches the expectation values of the observables and their distributions are computed in the course of the single spin flip quantum Monte Carlo algorithm for the coefficients of the equations of motion at each time step, while for the field quenches there is a direct functional relationship to the variational parameters. For the modulus of the magnetization this relationship reads with the eigenvalues µ m x of μ x according to equation (24) and for the correlation function between nearest neighbours it is with the eigenvalues ε n xx of Ĉ xx nn according to equation (25). The distributions of µ m x and ε n xx at time t are given by The correlation function between spins that are not nearest neighbours can be computed according to The C m n ( , ) xx r also have to be determined with RES. To decide whether the system thermalizes or not the expectation values of the observables as well as their distributions in the stationary state have to be compared to their counterparts for the thermal system. As the exact computation in the diagonal ensemble would require the knowledge of the full spectrum of the Hamiltonian after the quench, we use time averages to approximate the expectation values and distributions in the stationary state according to   with ρ = Ψ Ψ t t t ( ) ( ) ( ) and  j the (possibly degenerate) eigenvalues of  . We computed the time averages for different interval lengths Δ t and different t 0 in the range of times that we can simulate (t 0 + Δ t < 25 for the 16 × 16 system, longer times for smaller systems) and found that the results are stable with respect to our tests.
The time averages of the observables are compared to the thermal expectation values of the system in equilibrium at the effective temperature T eff attributed to the quench, which are given by The thermal distribution of the eigenvalues  j of  reads The expectation values and distributions for the system in thermal equilibrium are computed with a cluster Monte Carlo algorithm in continuous imaginary time 71 .

Results and Discussion
Before we apply the rt-VMC algorithm to large system sizes, we consider a system of size 4 × 4, whose time evolution can also be computed with numerical integration of the Schrödinger equation, and compare the rt-VMC results to the exact results to benchmark the algorithm. We do this exemplarily for the rescaled modulus of the magnetization. Results are presented in Fig. 4 for (a) interaction quenches and (b) field quenches. (a) i. and (b) i. contain a comparison between the exact time evolution (green) and the rt-VMC time evolution (red). Time averages are represented by the dashed lines of the respective colours, while the black dashed line is the thermal expectation value for the system in equilibrium at the temperature attributed to the quench. In (a) ii. and (b) ii. we compare the time-averaged distributions of µ m x of the exact diagonalization (green) and the rt-VMC (red) to the thermal distribution (black). We observe that after the interaction quenches the shape of the curves of the time evolution is close to harmonic oscillations with time-dependent variations of the amplitude. As long as ( / ) crit , the rt-VMC algorithm with the Jastrow ansatz allows a good description of the time evolution. There are differences between the frequencies and between the amplitudes, which increase for stronger interaction quenches when the system is driven closer to its phase transition. Despite of these deviations the time averages as well as the time-averaged distributions of the exact and the rt-VMC time evolution still show a very good agreement except for the quench (0; 3.5) → (1; 3.5), which drives the system close to its phase transition. Deviations from the thermal values on the other hand can already be observed beginning from the quench (0; 5) → (1; 5). For the field quenches the rt-VMC results for the time evolution after the quench show an even better agreement with the exact time evolution. The frequency of the oscillations is reproduced with high accuracy even for strong quenches. With increasing quench strength small deviations of the amplitude of the oscillations can be observed, but the time averages and the time-averaged distributions of the rt-VMC and the corresponding results of the exact time evolution coincide for all the quenches we studied. As in case of the interaction quenches there are increasing deviations between the time averages after the quench and the results for the thermal system in equilibrium at T eff with increasing quench strength. In order to quantify the deviations between the thermal expectation values and the time averages after the quench as a function of the quench parameters we compare in Fig. 5(a) i. and (b) i. the thermal expectation values for the system in equilibrium at the temperature T eff attributed to the quench (black) and the time-averaged values of the exact time evolution (green) and of the rt-VMC time evolution (red). The relative error between the time-averaged expectation values of the exact and the rt-VMC time evolution is shown in (a) ii. and (b) ii.. For the interaction quenches we observe an excellent agreement while  h J / 4 with a deviation of less than 2%. For the field quenches we find that beginning from h/J ≈ 0.75 deviations increase, but even for h/J = 1.5 they do not exceed 1.5%. For both quench protocols our rt-VMC method thus allows us to compute the time averages of the observables as well as their distributions with high accuracy for a wide range of ratios h/J. We may derive two main results from our studies of the 4 × 4 system: First we observe that for this small system size significant deviations between time-averaged results and the thermal results exist when the system is quenched close to its phase transition. Second there is a very good agreement between the time-averaged values of the exact time evolution and the rt-VMC time evolution for a wide range of ratios h/J. This agreement does not just concern the time averages of the observables, but also the underlying distributions. Only for strong quenches deviations between the time averages of the exact time evolution and the rt-VMC time evolution can be observed, but these deviations are much smaller than the deviations from the values of the system in thermal equilibrium.
We now apply the rt-VMC method to larger system sizes. To make predictions for the system in the thermodynamic limit, we computed results for the rescaled modulus of the magnetization and the correlation function between nearest neighbours after interaction quenches and after field quenches for systems of size L = 8, 12 and 16. Figure 6 shows our results for the time evolution of the observables for the different quench protocols we considered as well as a comparison between their time averages and the thermal expectation values for the system in equilibrium at T eff as a function of h/J after the quench. For the interaction quenches we confine us to values h/J > (h/J) crit , as for smaller ratios h/J the accuracy of the rt-VMC with the Jastrow ansatz decreases, while for the field quenches we study quenches with < h J h J / ( / ) 1 2 crit as for stronger quenches the system is predicted to be driven from the ferromagnetic to the paramagnetic phase. We observe that after the quenches the observables quickly approach a stationary value and in the following oscillate around it. This stationary value is in good approximation constant for the simulated time intervals. For small quenches both the amplitude and the frequency of the oscillations are almost constant. For larger quenches modulations of the amplitudes occur.
Scientific RepoRts | 6:38185 | DOI: 10.1038/srep38185 Considering the time averages and the thermal expectation values as function of h/J after the quench we observe that the relative shape of the curves is preserved when the system size is increased. There is a good agreement between the time averages and the thermal expectation values for small quenches with continuously increasing deviations for stronger quenches. While for the interaction quenches the observed deviations are small, in case of the field quenches the deviations become significant when the system is quenched closer to the phase transition. Comparing the deviations for strong field quenches one observes that they do not decrease with the system size, so that we assume them not to be caused by finite size effects.
Up to this point we have only considered the time-dependent expectation values of the observables and compared their time averages to the thermal expectation values. To decide whether the stationary state of the system after the quenches can be described by the CGE not only the expectation values of the observables have to coincide, but also their distributions. The study of the distributions is thus especially important for the interaction quenches and small field quenches, for which we have found a good agreement between the thermal and the time-averaged expectation values of the rescaled modulus of the magnetization as well as of the correlation function between nearest neighbours. [The strong variations in the distributions of the correlation function between nearest neighbours for ε n xx close to + 1 are not due to deficiencies of the algorithm, but are intrinsic to the system. Values of ε n xx close to + 1 correspond to small kink numbers n. x and ε n xx grows linearly with the system size, the resolution of the distributions becomes higher with increasing system size. For this reason we show in Fig. 7 the distributions of µ m x and ε n xx for the 16 × 16 system, i.e. the largest system size we can simulate. The quench protocols for the interaction and the field quenches are the same as in Fig. 6. The distance of the end point of the quenches from the phase transition is reduced from top to bottom, i.e. for the interaction quenches the ratio h/J after the quench is decreased, while for the field quenches it is increased. For small interaction quenches we observe a very good agreement between the thermal and the time-averaged distributions. Deviations increase when the system is quenched closer to the phase transition, but the shape of the thermal curves is well reproduced by the time averages after the quenches. Comparing the distributions for the 16 × 16 system to those of the 4 × 4 system in Fig. 4 we observe that the agreement is better than for the smaller system. For the field quenches we find strong deviations between the thermal distributions and the time-averaged distributions. In the distributions deviations can already be observed for quench protocols for which the time averages still agree with the thermal expectation values. These deviations could not be seen in the distributions of the 4 × 4 system in Fig. 4 due to the lower resolution of the distributions caused by the small system size. In general we observe that the time-averaged distributions after the field quenches are wider than the thermal distributions. The positions of their maxima are almost the same as for the thermal distributions, but the maxima are less pronounced. In addition the time-averaged distributions after the quenches show an increased probability to find the system in the fully ordered state (µ = ±1 m x or ε = +1 n xx respectively) compared to the thermal system. As these states are the initial state of the system this means that the system does not lose the memory of its initial state which contradicts to thermalization. We can thus state that for the system sizes and time scales that we can simulate the system does not thermalize after field quenches.
Up to this point we have only given a qualitative discussion of the distributions of the observables. In order to compare the deviations as a function of the system size and make predictions for the system in the thermodynamic limit we introduce the measure for the deviations between the thermal and the time-averaged distributions j is the number of different eigenvalues of  . The normalization is with respect to the maximum of the thermal distribution for the considered quench protocol and system size. We computed µ ∆  ( ) x and ∆Ĉ ( ) xx nn for the system sizes L = 4, 8, 12 and 16 for the described quench protocols and did a finite size scaling to conclude to the system in the thermodynamic limit. Figure 8 shows i. µ ∆  ( ) x and ii. ∆Ĉ ( ) xx nn after (a) the interaction and (b) the field quenches for the different quench protocols. The results are plotted as a function of the inverse system size 1/N. We find that for the interaction quenches the deviations between the thermal and the time-averaged distribution decrease with increasing system size. We thus conclude that the observed deviations of the expectation values and the distributions between the thermal and the time-averaged system will further decrease with increasing system size and that the system will thermalize in the thermodynamic limit. Before we make a statement on the results of the finite size scaling for the field quenches we discuss the effect of the spontaneous symmetry breaking in the ferromagnetic phase for the system in the thermodynamic limit. For the finite system sizes considered in our numerical studies the initial state of the field quenches is the symmetric superposition of the two fully magnetized states according to equation (10). In the thermodynamic limit on the other hand µ Ψ Ψ = ± 1 For the rescaled modulus of the magnetization µ  x and the correlation function Ĉ xx r the first and the second expectation value in the sum give the same result. Their sum thus just corresponds to the time evolution starting from ↑↑ … ↑↑ x . Applying a series expansion of the time evolution operator one can easily show that the two remaining matrix elements vanish in the thermodynamic limit, so that the time-evolved expectation values of µ  x and Ĉ xx r are independent of the initial state for the system in the thermodynamic limit. Simulations for finite For the interaction quenches the system does not leave the paramagnetic phase, while field quenches will drive the system from the ferromagnetic into the paramagnetic phase for ⪆ h J h J / ( / ) 1 2 crit . We observe continuously increasing deviations between the thermal expectation values and the time averages after the quenches the closer the system is quenched to its phase transition.
Scientific RepoRts | 6:38185 | DOI: 10.1038/srep38185 system sizes indeed show that the difference between the time evolutions starting from the two different initial states decrease with increasing system size. Our results for the field quenches can thus be used to conclude to the system in the thermodynamic limit. The finite size scaling in Fig. 8(b) shows that both µ ∆  ( ) x and ∆Ĉ ( ) xx nn increase with the system size or at least do not decrease. We interpret this as indication that the system will not completely thermalize in the thermodynamic limit either.
As the last point of our studies we consider correlation functions between spins which are not nearest neighbours. Due to the long-range order in the ferromagnetic phase the correlation function of two spins will not vanish at large distances but reach a constant non-vanishing value given by the square of the expectation value of the rescaled modulus of the magnetization. In the paramagnetic phase on the other hand there is no long-range order and the Hamiltonian is gapped, i.e. there is an energy gap between the ground state and the first excited state. It has been shown that in the ground state (T = 0) of gapped quantum systems correlations decay exponentially with the distance 83,84 . The same holds for thermal states in the paramagnetic phase 85 . In contrast to this there are no analytic expressions for the decay of the correlations in phases with long-range order, i.e. the ferromagnetic phase in case of the 2D-TFIM. We compute the long-range correlations for the 4 × 4 and the 12 × 12 system. The constraint to the 12 × 12 system is due to the computational effort in the computation of the C m n ( , ) xx r . The rt-VMC results for the 4 × 4 system are compared to exact results. Figure 9 shows the results for the long-range correlations for i. the 4 × 4 system and ii. the 12 × 12 system after (a) interaction quenches and (b) after field quenches. The distance d between two sites is measured in the Manhattan metric, i.e. if r = R − R′ defines the relative position of the sites we have with D the dimensionality (here D = 2). Thus for a given distance d there may be several r. The strength of the correlation between two sites depends on the number of shortest paths between them, which depends on r. We compare the time-averaged rt-VMC results after the quenches (red) to the correlation functions for the system in equilibrium at the temperature T eff attributed to the quench (black). For the 4 × 4 system we additionally show results of the exact time evolution (green). For the interaction quenches we do a least-square fit with a · e bd . We observe that the decay of the correlations in the thermal system as well as in the quenched system is well described by the fit curves in agreement with the analytic results. For the field quenches we have added the square of the rescaled modulus of the magnetization. We observe that within the numerical error its value and the correlation function for large d coincide as predicted. Comparing the decay of the correlations for the thermal system to the decay in the quenched system we find for the interaction quenches a good agreement both for the 4 × 4 and the 12 × 12 system. The curves show only small deviations and have the same shape going to 0 for large distances d. The rt-VMC results show a very good agreement to the exact results apart from the interaction quench (0; 3.5) → (1; 3.5). For the field quenches the differences between the time-averaged results after the quenches and the thermal results are more significant. For small quenches the shapes of the curves are very similar and the values coincide within the numerical error. For larger quenches we still observe a good agreement between the rt-VMC results and the exact results, but here larger deviations from the thermal curves arise. Although the absolute differences between the curves for a given value of d are not too large, the decay of the correlations as function of d is different for strong field quenches. After the field quenches the correlations decay faster with d than for the system in thermal equilibrium and the stationary long-range value of the correlations is reached at larger distances d and is lower than for the thermal system.

Related work
Recently an exact theorem on generalized thermalization in D-dimensional quantum systems in the thermodynamic limit has been formulated 55 . The theorem states that generalized thermalization can be observed if the state of the system is algebraically sizably clustering. It also holds for exponentially sizably clustering states. Then the stationary state of the system can be described by a GGE which has to take into account all local and quasilocal charges of the system. For non-integrable systems, for which the total energy is the only conserved quantity, the generalized thermalization reduces to thermalization with a stationary state according to the CGE. We now discuss the exact theorem with respect to the 2D-TFIM. Considering the 2D-TFIM one has to distinguish between the ferromagnetic and the paramagnetic phase. In the paramagnetic phase the 2D-TFIM is gapped and there is no symmetry breaking both for finite system sizes as well as in the thermodynamic limit. As the ground state of gapped quantum systems is sizably exponentially clustering 83,84 the exact theorem can be applied to the 2D-TFIM in the thermodynamic limit after the interaction quenches in the paramagnetic phase and predicts thermalization. In our numerical studies we have indeed observed a very good agreement both between the time-averaged observables and their thermal counterparts as well as between the distributions for small quenches, i.e. large ratios h/J, also for the finite system sizes that we can simulate. For larger quenches closer to the phase transition we have found deviations between the time averages and the thermal values, but the finite size scaling shows that they decrease with the system size. Our results for the interaction quenches in the paramagnetic phase are thus in agreement with the exact theorem. In the ferromagnetic phase on the other hand the Hamiltonian of the system is not gapped in the thermodynamic limit. The spin flip symmetry is spontaneously broken and long-range order exists, so that all spins of the system are correlated to each other and the correlations do not cluster. An analytic expression for the shape of the decay of the correlations has not been found yet, thus it cannot be decided whether the exact theorem can be applied in the ferromagnetic phase or not.

Summary and Conclusion
We have studied the quantum relaxation of the 2D-TFIM after global interaction quenches in the paramagnetic phase and global field quenches in the ferromagnetic phase using a newly developed rt-VMC method which allowed us to explore system sizes and time scales that have not been accessible before. In order to answer the question whether this two-dimensional, non-integrable system thermalizes or not we compared time-averaged results after the quenches to results for the system in thermal equilibrium at a temperature defined by the excess energy after the quench. We found that the presence or absence of thermalization depends crucially on the quench protocol and the initial state. For the interaction quenches there is a good agreement between the results for the quenched system and the thermal system in accordance with a recently formulated exact theorem for systems in  1.25). One observes that for the interaction quenches the differences between the distributions decrease with increasing system size, while for the field quenches they increase or remain almost constant. the thermodynamic limit. Deviations are only observed for strong quenches ending in the vicinity of the phase transition. Finite size scaling suggests that these deviations should vanish in the thermodynamic limit. In contrast to this we have found significant deviations between the thermal results and the time-averaged results after the field quenches in the ferromagnetic phase. These deviations become especially clear comparing the distributions which show deviations already for small quenches for which the thermal expectation values and the time averages still agree. The shape of the distributions indicates that the system does not completely lose the memory of its initial state during the relaxation process, which is a clear contradiction to thermalization. Finite size scaling shows that the deviations do not decrease with the system size either. Although we currently cannot give an explanation for the observed deviations between the thermal system and the time averages after the field quenches in the ferromagnetic phase, we assume that they might be related to the long-range order or the structure of the spectrum of the Hamiltonian without an energy gap between the ground state and the first excited state.