{'fraction_non_alphanumeric': 0.03567716294711098, 'fraction_numerical': 0.010394374808926934, 'mean_word_length': 4.654451166810717, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Strong and weak thermalization of infinite non-integrable quantum systems

When a non-integrable system evolves out of equilibrium for a long time, local observables are expected to attain stationary expectation values, independent of the details of the initial state. However, intriguing experimental results with ultracold gases have shown no thermalization in non-integrable settings, triggering an intense theoretical effort to decide the question. Here we show that the phenomenology of thermalization in a quantum system is much richer than its classical counterpart. Using a new numerical technique, we identify two distinct thermalization regimes, strong and weak, occurring for different initial states. Strong thermalization, intrinsically quantum, happens when instantaneous local expectation values converge to the thermal ones. Weak thermalization, well-known in classical systems, happens when local expectation values converge to the thermal ones only after time averaging. Remarkably, we find a third group of states showing no thermalization, neither strong nor weak, to the time scales one can reliably simulate.

In the 19th century, statistical mechanics was developed as a microscopic model explaining the fundamental results of thermodynamics. Starting in the early 20th century, quantum statistical mechanics has been developed to describe thermodynamics of quantum systems. However, current experiments [1, 2] with ultracold atoms have attained a level of isolation and control of parameters that forces us to address not just the equilibrium properties of quantum systems but also the question of how and why these systems relax to equilibrium starting from a nonthermal state. For classical systems, Boltzmann's molecular chaos assumption provides a quantitative tool to describe this relaxation, and ergodic properties of the system give the explanation. If the time dynamics of a system explores all states with a given energy with uniform probability, then the long-time average of any observed quantity will approach the expectation value for that quantity in the microcanonical ensemble.
Similarly, the emergence of a thermal bath [9] in a quantum system has been justified under the assumption that the long-time average of typical initial states produces a density matrix that approaches the thermal average [7,9,10]. It has been proposed, however, that thermalization may happen without any time average in the quantum case [11]; indeed, it is possible that, due to quantum entanglement, starting from fixed, non-thermal initial conditions, the reduced density matrix at a given time t on a given region A, such that A is small compared to the system size, will converge to the thermal expectation value at long times t. Such a phenomenon, which we call "strong thermalization", cannot occur for a classical system as it relies on the quantum mechanical fact that even if the global density matrix is a pure state, the reduced density matrix may be a mixed state.
For integrable systems the existence of local conserved quantities prevents relaxation to the thermal state, which is constrained only by energy. Instead, it has been sug-gested that these systems will relax to a state described by a generalized thermal ensemble, compatible with the set of conserved quantities [3,6,7]. In an infinite nonintegrable system, there exist an infinite number of conserved quantities, such as the powers of the Hamiltonian. However, not being local, they are not expected to prevent the thermalization of local observables [7,9,11].
Remarkably enough, a recent experiment [1] observed no signs of thermalization after long evolution in a nearly integrable case. The results could be in part understood with the hypothesis of relaxation to a generalized thermal ensemble [3] in an integrable case, but from the theoretical point of view it is not clear why no thermalization was present away from true integrability.
Various theoretical studies have tried to elucidate the question of whether or under which conditions thermalization will occur in non-integrable models [4,5,7,8,[12][13][14]. The study has to resort to numerical techniques, given the lack of analytical solutions. Moreover, the intrinsic computational complexity of simulating large quantum systems limits the affordable studies to finite systems or short times, so that results showing nonthermalization cannot be extrapolated to infinite times or system sizes.
Most of the studies have tried to decide whether a given model thermalizes or not as a function of the Hamiltonian parameters [4,5,8,12]. Recent works [15,16]. have also analyzed the link between appearance of quantum chaos and thermalization in these systems.
Here we consider the question of relaxation using a fixed, non-integrable Hamiltonian, but a range of initial states. We discover a rich phenomenology of different relaxation regimes. The analyzed system is an infinite translationally invariant spin chain with a nearestneighbor interaction. Thanks to a recently developed algorithm [17], we are able to explore its dynamics at relatively long times. Furthermore, the method gives us access to the whole reduced density matrix for a block of few particles. Instead of analyzing the behavior of individual expectation values, as previous numerical studies, we may then quantify the degree of thermalization as the distance between the reduced density matrix of the evolved state and the thermal state, ρ th (β) = e −βH tr(e −βH ) , corresponding to the same energy. We estimate this distance both for the instantaneous reduced state, ρ(t), and for the time averaged oneρ(t) ≡ 1 We find two clearly distinct thermalization regimes. Some states show strong thermalization. Their instantaneous expectation values converge (fast) to the thermal ones, without the need to consider the time average. For other states we encounter instead weak thermalization. We have strong numerical evidence that such states do not relax even at very long times. Nevertheless, the time averaged observables do relax to the thermal values. Finally, we also obtain weaker numerical evidence that for some third group of states no relaxation occurs, at least to the longest time scales we can reliably simulate. Strikingly, this situation occurs even in regions of energy for which the spectrum of the system is shown to be chaotic (see Supplementary material). Since our analysis is numerical, it may be argued that the conclusions are only valid for a limited range of time. However, the study of the properties of the method [17], as well as the detailed analysis of errors in the current study (described in the Supplementary material) ensure that the reliability of our results extends to longer times than that of previous studies, enough to give strong evidence about the existence of these different quantum thermalization regimes.
The first regime, or strong thermalization, is illustrated by the initial state |Y + (Fig. 1), in which all spins are aligned along the positiveŷ direction, corresponding to zero initial energy, and thus β = 0. The density matrix at any time converges fast to the thermal distribution at the same energy. As a consequence, the time averaged density matrix also converges to the thermal values.
We identify an example of the second regime, or weak thermalization in the initial state |Z+ (Fig. 2), in which all spins are initially aligned along the positiveẑ direction. For this state, whose energy matches that of a thermal state with β ≈ 0.7275, no relaxation at all is observed for long times. Instead, the distance between the evolved and the thermal reduced density matrices seems to oscillate strongly for arbitrarily long times. All expectation values of the form σ a ⊗ σ b ⊗ σ c show similar non-damped irregular oscillating behavior. Since the entropy increases linearly with time in the |Z+ state, the behavior is not due to a lack of propagating excitations (see Supplementary material). The time averaged density matrix shows instead an only slightly oscillating behavior, compatible with a very slow convergence to the thermal state, at a rate 1/ √ t. This rate is characteris- The main plot shows the distance between the three-body reduced density matrix at instant t and the corresponding thermal state (β = 0). The error bars show the difference between the result with the largest bond dimension (D = 120) and with a lower one (D = 60), which gives us an estimate of the truncation error. The right inset shows the distance between the timeaveraged reduced density matrix and the thermal state. We superimpose a fit to a curve decaying like b/ √ t for long times (solid black line). The left inset shows the difference between the thermal expectation values and the time dependent single body observables σx (blue solid line), σy (dashed green) and σz (dash-dotted black line). Convergence to the thermal state is observed in all three plots.
tic of diffusive relaxation in one dimension, and has been seen elsewhere [12].
Finally, if the system is started in the |X+ state ( Fig. 3), with all spins aligned along the positivex direction, positive energy per particle, and corresponding to a Gibbs state with β ≈ −0.7180, the reduced density matrix shows initially a fast relaxation. However at long times the distance between the evolved state and the thermal one is different from zero, signaling that one or more expectation values of local observables have not converged to the thermal averages. We observe that the density matrix has not reached a stationary value, either. Given the numerical character of the study, this observation could indicate an absence of thermalization, but also a much slower one. It is an intriguing case, worthy of further, ideally experimental, study. Even in average, we detect no thermalization of the |X+ state, up to the time scales we are able to explore.
The appearance of the different thermalization regimes does not require any fine tuning of the initial conditions or Hamiltonian parameters. On the contrary, by analyzing the evolution of different sets of product initial states, rotating the initial polarization of the spins fromŷ tô z, we see both strong and weak thermalization regimes over a range of parameters, separated by a transition. As shown in Fig. 4, we observe that weak thermaliza- tion appears approximately half-way (corresponding to β ≈ 0.3502), with strong thermalization for those states closer to |Y + . Similarly, we observe a transition between the strong and the non-thermalizing cases (see Supplementary material).
The different regimes survive also under changes to the Hamiltonian parameters, showing that they are also not an isolated phenomenon for the particular values we chose. We have checked (see Supplementary material) that weak thermalization becomes prevalent if we decrease the magnitude of the transverse magnetic field g, while as we increase it, strong thermalization shows up in a larger fraction of the initial states. The same effect occurs when we decrease the magnitude of the parallel magnetic field, h.
Our data show that the thermalization process of a quantum system is a much richer phenomenon than its classical analogue. The appearance of the various thermalization regimes cannot be linked exclusively to the integrability of the Hamiltonian, since we observe these different regimes for a fixed model, not close to any integrable limit. In particular, some initial states thermalize strongly, i.e. at the level of the instantaneous expectation values, while others do it weakly, or only in average. Some other states seem to retain memory of the initial configuration for much longer and no type of relaxation could be proved. The non-relaxing configurations would No thermalization is observed in the instantaneous density matrix, nor in the time averaged one (right inset). On the latter, we superimpose a fit of the computed values to a time dependent function, which asymptotically tends to a constant (0.03). The plotted results correspond to a bond dimension D = 240, while the distance to the results with D = 120 is shown as error bars. In the left inset we plot, for the observables that determine the one site reduced density matrix, σx (blue solid line), σy (dashed green) and σz (dash-dotted black), the difference with respect to the thermal values. We observe that σx is the responsible for the lack of thermalization, while all the other expectation values converge to the thermal averages. Studying the evolution of only a few local expectation values may not suffice then to detect the nonthermalization. be good candidates for an experimental study of thermalization.

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 [34]. 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 θ = π 2 , φ = 0 (|X+ ), θ = φ = π 2 (|Y + ) and θ = 0 (|Z+ ).
Using the recently developed folding method [17], 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 ρ th (β) 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(ρ 1 , ρ 2 ) ≡ ρ 1 − ρ 2 op , which in this case coincides with the maximum eigenvalue in absolute value of the difference ρ 1 −ρ 2 .
We are grateful to Maciej Lewenstein and Jens Eisert for suggestions. We acknowledge support by the DFG through Excellence Cluster

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+ and |X+ , and at these energies there could be an integrable effective model [20], even for g = −1.05, h = 0.5. The spectral properties discussed in The left plot shows that the level spacing distribution, in the finite system of length L = 14, for a window of energy centered on E/N = 0.91 and of width ∆E/N = 0.21, is closer to Wigner-Dyson than Poisson. Nevertheless, the product state with the same energy per site in the infinite chain shows no thermalization, as the right plot illustrates. The magenta diamonds represent the distance between the thermal reduced density matrix and the time averaged ρ for this particular state as a function of time. The strongly thermalizing state |Y + (black triangles) and the non thermalizing |X+ (red crosses) are also shown for comparison. 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.

THE NUMERICAL METHOD
The time evolution of an infinite 1D quantum system is simulated numerically within the Matrix Product States [21-26] (MPS) formalism, using the new technique introduced in [17]. 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 Ψ(t)|O|Ψ(t) can then be expressed as a two dimensional tensor network ( Fig. 8(a)), constructed from the a Suzuki-Trotter expansion of the evolution operator [27,28]. Each discrete time step corresponds then to a sequence of matrix product operators (MPO) [29]. 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 [30,31], 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, L| and |R . 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 Ψ(t)|O|Ψ(t) = Φ|(O|Ψ(t) ⊗ |Ψ(t) ), where |Φ is a product of unnormalized maximally entangled pairs between each site of the chain and its conjugate (see [17] for a detailed discussion of the algorithm). This is equivalent to grouping together tensors that correspond to the same time step in Ψ 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 β, E th (β) (Fig. 9). From these data we may compute, for each one of the initial states we have studied, the value of β 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 β, 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.

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 δ, which determines the size of the time step, or by using a Suzuki-Trotter expansion which is exact to a higher order in δ. 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 δ = 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 [32]. 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 [17]).
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 [31], 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 ≃ 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+ and |X+ , with bond dimension D = 240, the relative error in the energy for the range of times we are analysing (respectively t ≃ 18 and t ≃ 12) is kept to only a few percent, consistent with the estimated truncation error. For the initial state |Y + , 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+ initial state, the deviation in energy at times t ≈ 10 reaches a 1%, which corresponds to a distance d(ρ(β 0 ), ρ(β ′ )) ≈ 9 × 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 ≈ 18, the largest distance grows to a maximum value of 0.06. For the |X+ initial state, the maximum deviation in energy, at t ≈ 12, corresponds to a distance d(ρ(β 0 ), ρ(β ′ )) ≈ 11 × 10 −3 , while the one we find at this same long time is 0.04

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.

Different initial states
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 (θ, φ). This state has energy per spin E = − cos 2 θ + g sin θ cos φ + h cos θ . States with spins polarised in the three orthogonal directions, |X+ , |Y + and |Z+ , 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 + , showing strong thermalization, all of the individual expectation values converge fast to the thermal ones. Instead, for the initial state |Z+ , we observe irregular oscilla- tions, showing no sign of damping, to the longest times we are able to simulate. For the initial |X+ 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 σ x 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 + to thermalize becomes longer.
From the iTEBD simulations we may analyze how the entropy of the half-chain increases in time for the nonthermalizing states (see Fig. 14). We notice that for the weak thermalizing initial state |Z+ the entropy of the half-chain grows linearly with time, so that the oscillatory behavior is not due to the absence of propagating excitations. For the initial state |X+ the entropy also grows linearly with time, but it does so at a faster pace. As shown in the figure, the iTEBD simulation can reproduce this growth until the time when the truncation error becomes dominant. A closer look at the data of this simulation shows that the distribution of Schmidt coefficients is very different for both cases, even at times For β = 0, we compare |Y + (red upwards-pointing triangles) and (θ = 0.4398,φ = π) (blue downwards-pointing triangles). For β = −0.7180, we compare |X+ (red x) to (θ = 1.90,φ = 0.031) (blue crosses). We show also the behaviour of the state with maximal β = 1.5859 (green stars), whose expectation values oscillate like for |Z+ , but in average converges very fast to thermal. All data correspond to D = 120 simulation.
when they attain a comparable entropy. If we analyze the Schmidt decompositions of both states at the time at which their entropy is S ≈ 0.8 (t = 1.1875 for |X+ and t = 8 for |Z+ ), we observe that the coefficient distribution for |Z+ has a much longer tail. This difference in the distribution can also be quantified by the 2-Rényi entropy S 2 = −log(trρ 2 ), which attains a lower value for the weak thermalizing state, S 2 (Z+) = 0.348, while S 2 (X+) = 0.611.
By rotating the inital state on the Bloch sphere, we observe a transition from the strong thermalizing |Y + to the weak one |Z+ , and also to the apparently nonthermalizing |X+ (Fig. 15), the latter occurring around φ ∈ [ π 6 , π 4 ] (β ∈ [−0.5382, −0.3915] ). We have also analysed the transition from weak thermalization |Z+ to non-thermalization |X+ (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 β of every initial product state. We discover that all the strong thermalizing states we have found have energies, and thus β, close to zero.
Instead, for larger β > 0, we observe oscillations and weak thermalization as in the |Z+ initial state. An interesting case is the state of maximum β (the product state with minimal energy), which we can identify, by studying the energy landscape over the Bloch sphere, at θ ≈ 0.43, φ = π. The dynamics of this initial state shows also weak thermalization, with strong oscillations of all The integrable case is represented always by black dots. In the three cases, for the integrable limit we see the relaxation of the reduced density matrix to a state that is different from thermal. We notice also that the initial state |X+ shows in the non-integrable case a similar behaviour to the h = 0 limit.
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+ and |Y + , 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 β, the state of the system does not relax to the same state, at least during the long range of times we simulate.

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 FIG. 20: Distance between instantaneous reduced density matrix and thermal state for N = 3, and different strengths of the Hamiltonian parameter g. The initial states vary from |Y + to |X+ , as depicted on the Bloch sphere. For constant parallel field h = 0.5, we compare the case g = −1.05 (upper right pane) to g = −0.5 (lower left) and g = −1.5 (lower right).

the generalized thermal ensemble [3]).
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 val- ues of g, the behaviour approaches strong thermalization in most cases.
If we do the same now for the transition to nonthermalizing as observed between |Y + and |X+ , 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+ and |X+ . We observe that the oscillating behaviour of the initial state |Z+ 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+ is quite different, and in all cases we observe thermalization of some observables while others deviate from the thermal expectation value.