On the size of the space spanned by a nonequilibrium state in a quantum spin lattice system

We consider the time evolution of a state in an isolated quantum spin lattice system with energy cumulants proportional to the number of the sites $L^d$. We compute the distribution of the eigenvalues of the time averaged state over a time window $[t_0,t_0+t]$ in the limit of large $L$. This allows us to infer the size of a subspace that captures time evolution in $[t_0,t_0+t]$ with an accuracy $1-\epsilon$. We estimate the size to be $ \frac{\sqrt{2\mathfrak{e}_2}}{\pi}\mathrm{erf}^{-1}(1-\epsilon) L^{\frac{d}{2}}t$, where $\mathfrak{e}_2$ is the energy variance per site, and $\mathrm{erf}^{-1}$ is the inverse error function.


Introduction
Nonequilibrium dynamics in isolated many-body quantum systems have been being perceived as a captivating area where to look for the missing link between the physics of systems in equilibrium and the fundamental equation of quantum mechanics, the Schrödinger equation. On the one hand, the thermodynamic limit (large number of sites in spin lattice systems or large number of particles in particle systems) opens the door to critical phenomena, like spontaneous symmetry breaking and phase transitions. On the other hand, the limit of large time is characterised by a loss of information that resembles the statistical equivalence of microscopic configurations in a system in thermal equilibrium.
A weak form of relaxation can be defined by taking the time averages of the expectation values of the observables and sending the time to infinity [1]. Generally, this limit exists both when the number of degrees of freedom is finite and in the thermodynamic limit, but the non-commutativity of the limit of infinite time with the thermodynamic limit makes the definition ambiguous 1 .
The time average of a nonequilibrium state has been extensively studied when the thermodynamic limit is taken after the infinite time limit, which is the setting of the quantum recurrence theorem [3]. In the opposite order of limits, the limit of infinite time usually exists in a stronger sense, without average, and indeed most of the studies have been focussed on the infinite time limit of the expectation values [4]. There is however additional information that can be extracted from the time average when the number of degrees of freedom is large. Specifically, we can use it to estimate the size of the space spanned by the state in a given time window. In this paper we show that, for a relevant class of systems characterised by extensive energy cumulants, this estimate is almost independent of the system details.

Physical setting
We consider a quantum spin lattice system with L d sites, d being the dimension of the lattice. The system is prepared in a pure state |Ψ 0 that time evolves under a spin-lattice Hamiltonian H: |Ψ t = e −iHt |Ψ 0 . We do not specify other details, but we assume that the energy cumulants, denoted by L d e n , are proportional to the number of the sites L d e n = ∂ n t t=0 In particular, (1) is satisfied in generic spin lattice systems as long as interactions and correlations decay sufficiently fast to zero with the distance 2 . Let us consider the time averaged state: For given L, since the local space is finite 3 , the spectrum of H is discrete, and hence the infinite time limit exists; it is given by the so-called "diagonal ensemble" [1,6] lim t→∞ρ where |E form a basis diagonalizing H. We consider here the opposite limit of finite t and large L. We wonder how big it is the dimension D t of the space spanned by |Ψ 0 in the time window [t 0 , t 0 + t]. Strictly speaking, this is given by the rank ofρ t 0 ,t . The latter is however sensitive to infinitesimally small perturbations which do not really affect the dynamics of the physical observables. It is then more useful to approximate the state up to a given accuracy so as to reduce the dimension of the subspace, still capturing the relevant part of the dynamics. We do it at the level of the time averaged state, introducing a low-probability cutoff t , possibly dependent on the width of the time window, for the eigenvalues ofρ t 0 ,t . This leads to the following definition where θ H (x) is the Heaviside step function, and λ t is the corresponding cutoff in the eigenvalues. Note that the rank ofρ t 0 ,t can be obtained as a limit: t . We will come back to the choice of t later; we focus first on the solution of (4) for given t .

Loschmidt echo
Finding a solution to (4) is a hard problem, but there are crucial simplifications in the limit of large L. As it will be clear in the next section, these simplifications can be traced back to the behaviour of the overlap between the state at different times We remind the reader that the square of the absolute value of the overlap is also known as Loschmidt echo [7,8] in the specific case when the backward evolution is generated by a Hamiltonian which |Ψ 0 is eigenstate of. The series expansion about t 1 = t 2 of the logarithm of the overlap can be written as follows where, by assumption, each order of the expansion is proportional to the number of the sites, i.e., it is "extensive". In quantum spin lattice systems, extensivity is a non-perturbative feature, indeed one generally finds By definition f (t) is a nonnegative function with a zero at t = 0. Generally it remains finite in the limit of infinite time, and it can exhibit non-analytic behavior, which has been a subject of intensive investigations since 2012 [9]. (The interested reader can find some numerical data, showing the behaviour of the Loschmidt echo in one-and two-dimensional lattices, in Refs [10,11].) A simple but powerful property that we are going to assume and exploit is that normally f (t) has a single zero on the real line.
In other words, the overlap is exponentially small (in the number of the sites) everywhere but in the neighbourhoods of t = 0. The time window where it is not exponentially small has to shrink to zero in the thermodynamic limit L → ∞. In that region, the series expansion in the time (6) can also be interpreted as an asymptotic expansion in the number of the sites.

Entropies
We compute here the moments of the distribution of the eigenvalues ofρ t 0 ,t , namely tr[ρ α t 0 ,t ] for integer α. First of all, we note that they are independent of t 0 , asρ t 0 ,t = e −iHt 0ρ 0,t e iHt 0 is unitarily equivalent toρ 0,t and the moments are invariant under unitary transformations. From now on we will writeρ t instead ofρ 0,t . We start with the second moment tr[ where we have formally replaced the function f (t) with its series expansion. The contributions to the integral coming from regions where y increases with L are subleading 4 , and, in turn, only the term proportional to e 2 survives the limit. We then find We stress again that here t is a finite nonzero parameter (this expression does not capture the behaviour for t ∼ L − d 2 ). From (9) we infer the asymptotic behaviour of the second Rényi entropy We now compute a generic moment tr[ρ α t ]. We have The considerations made for α = 2 hold true also for α > 2: for any given integer α, we can neglect the cumulants higher than the second one at the price of introducing a relative error O(L −d ). We can then write the moments as follows 5 Thus, the asymptotic behaviour of the Rényi entropies S α [ρ t ] = 1 1−α log tr[ρ α t ] reads as As shown in Appendix A, the leading correction O(L − d 2 ) comes from a more careful integration over the variable that is not modulated by the gaussian.
A straightforward application of the replica trick gives the von Neumann entropy S vN

Distribution of eigenvalues
Since the support of the distribution of eigenvalues is bounded, the moments, and, in turn, the Rényi entropies, characterise the distribution completely (Hausdorff moment problem). The distribution can then be reconstructed using the approach proposed in Ref. [12]. To that aim, we define P ρ (λ) as the (nonnormalized) distribution of the eigenvalues λ j of the density matrix ρ: P ρ (λ) = j δ(λ − λ j ). It turns out that Φ ρ (λ) ≡ λP ρ (λ) can be written as Notwithstanding we computed only the leading order of the asymptotic expansion of the moments in the limit of large L, we expect the corrections to be subleading almost everywhere but close to λ = 0. Thus, we can use (12) to reconstruct the asymptotic distribution. Plugging (12) into (15) gives where Li1 /2 (x) is the polylogarithm of order 1 /2. The distribution of eigenvalues is then , with Γ(x) the gamma function. Note that, by definition, Pρ t (λ) is normalised to the dimension of the Hilbert space; the asymptotic result (17), on the other hand, despite capturing all the moments 1 0 dλPρ t (λ)λ n with integer n > 0, has a divergent integral. This is because the behavior for λ ∼ o(L − d 2 ) goes beyond the leading order of the asymptotic expansion.
The distribution Φρ t (λ), shown in figure 1, is instead correctly normalised Again, we expect the corrections to this asymptotic result to mainly affect the behaviour of Π(x) close to x = 0. In view of the universality of (18), we come to the conclusion that the large L behaviour of the distribution of the eigenvalues of a time averaged state is a fundamental feature of quantum spin lattice systems with extensive energy cumulants.

Effective rank of the time averaged state
Having computed the asymptotic behaviour of Φρ t (λ) and of Pρ t (λ) in the limit of large number of sites, we have access to the asymptotic solution of (4), provided that λ t ∼ O(L − d 2 ).
which are expected to be valid as long as t does not approach zero in the thermodynamic limit. Figures 1 and 2 provide a graphic representation of t and D ( t) t . Assuming t small, we can expand the inverse error function as follows so we find The next step is to establish a connection between the cutoff t and the error on the state. To that aim, let us introduce a tiny time scale δt such that t /δ is integer. The time average in [0, t] can be seen as the mean of the sample consisting of the time averaged states over the time windows [(n − 1)δt, nδt], with n ∈ [1, t /δt]. For asymptotically large L, it is reasonable to expect the errors produced by truncating the spectrum ofρ t to be randomly distributed over the sample. Under this assumption, since the variance of two independent random variables is additive, the error t of the mean of the sample has to scale as t ∼ δt δt /t, where δt is the truncation error in each time slice. This gives We are interested in observables that have a well-defined expectation value in the thermodynamic limit (for example, local observables), for which the variation of their expectation value in a tiny time window approaches zero as the width of the interval shrinks to zero. Thus, if δt is small enough, δt can be identified with the effective error on the time evolving state, and D ( δt ) t|δt is the effective dimension we are looking for. This observable-dependent step gives a meaning to "relevance": our truncation would be inadequate if we would consider observables that have different expectation values in |Ψ t and |Ψ t+τ L , with lim L→∞ τ L = 0. Note also that we are not allowed to choose a cutoff δt approaching zero in the thermodynamic limit because that would require the knowledge of the next orders of the asymptotic expansion, which is trickier (see Appendix A).
The projection on a space with size (22) is arguably the best approximation with error δt if we have only access to the time averaged state; but this is far from being optimal. It is more convenient to merge the reduced spaces in the time slices [(n − 1)δt, nδt], each of which is the span of D states. The size of the resulting space is bounded from above by the total number of elements, which is In the limit δt → 0, we reinterpret δt → as the truncation error on the state 6 ; we finally find the upper bound D ( ) t , which is tighter than (22). We can now state our main result: The size of the space that is approximately spanned by a nonequilibrium state in the time window [t 0 , t 0 + t] with error on the state is asymptotically bounded from above by 6 As in the previous discussion, we are restricting our attention to the system properties compatible with |Ψt Ψt| L→∞ ∼ lim δt→0 limL→∞ρ t,δt , which would have been an identity if the order of limits were reversed.

Numerical simulations of nonequilibrium dynamics
The result (23) is rather suggestive if reconsidered in the context of simulations of out-ofequilibrium many body quantum systems. First of all, the relevant space where the dynamics take place is proportional to the square root of the logarithm of the Hilbert space! In addition, within the assumptions of our calculation, a Lieb-Robinson velocity v LR generally exists [13,14], which bounds the speed at which information propagates throughout the lattice. Consequently, the dynamics of a compact subsystem of size d in the time window [0, t] display exponentially small finite-size effects, provided that L + 2v LR t: we can replace L by + 2v LR t + R in (23), committing an error that is exponentially small in R. In conclusion, in the limit of large time, the relevant size of the space does not grow faster than ∼ t d 2 +1 . An algorithm reducing the dynamics onto the relevant subspace would allow for investigations in much wider time windows than those accessible nowadays through state-of-the-art techniques (which, in spin chains, could be time-dependent density matrix renormalisation group (tDRMG) [15] and infinite time-evolving block decimation (iTEBD) [16] algorithms).

Quantum speed limit
Our findings are complementary to the studies on the minimum time τ required for arriving to an orthogonal state -"the quantum speed limit"; we mention here just the classical result by Mandelstam and Tamm [17] τ ≥ πL − d /2 /2 √ e 2 7 . Notwithstanding the similarity between this limit and our estimate, their meaning is quite different. In our approach a new orthogonal state starts counting in D ( ) t after developing a significant overlap, not necessarily equal to 1, with the time evolving state. In addition, the analogue of the quantum speed limit in the quantum systems considered is the time needed to have a so-called "dynamical phase transition" [9]. Generally, the latter time remains nonzero even in the thermodynamic limit, and, in a hypothetical case where it does not, the hypotheses behind the asymptotic expansion that we carried out would not be fulfilled.

Conclusion
We studied the nonequilibrium time evolution of a state under a Hamiltonian of a general quantum spin lattice system with energy cumulants proportional to the number of the sites. We have computed the leading order of the asymptotic expansion of the distribution of the eigenvalues of the time averaged state over a fixed time window in the limit of large number of sites. We used the asymptotic distribution to determine an upper bound for the dimension of the space visited by the state, and we have found that it is proportional to the square root of the logarithm of the Hilbert space. It would be interesting to generalise our results to the time evolution of critical states with super-extensive energy cumulants, like the ones considered in Ref. [5]. Finally, our bound does not distinguish chaotic systems from integrable ones (see also Appendix C), which are expected to time evolve with a lower complexity [19] (see also [20] and references therein for more recent investigations); we wonder whether a tighter bound could be established in integrable systems.

Acknowledgements
These notes have been prepared for the course "Quench dynamics and relaxation in isolated integrable quantum spin chains", held in IPhT Saclay from 25/01/2019 to 22/02/2019.

A Leading correction
In this appendix we work out the leading correction to the asymptotic behaviour of the moments tr[ρ α t ] when the number of the sites is large. As suggested by (8), that is a relative correction O(L − d 2 ) that comes from the constraint on the integration variables, which have to sum to zero. From the representation (11) it follows that the higher order cumulants start contributing at relative order O(L −d ), therefore our starting point is the second line of (12). The leading correction becomes visible if we integrate first in τ α . To that aim, we must move the first integral to the right, like we did in (8). The result is The leading term, which we already computed, comes from the 1 in the round bracket; the other term includes the leading correction. The domain of integration is rather complicated, but, for an asymptotically large number of sites, it can be extended to infinity. In addition, by reverting the sign of all τ j , we realise that the contribution from the term proportional to min(T α ) is equal to the contribution from the one proportional to − max(T α ). Thus we have This expression can be simplified further by defining the new variables and summing over the α − 1 possibilities for which max(T α ) = e − 1 /2 2 max(0, y j ). We finally obtain We have not found a closed form expression for the gaussian integral on the right hand side of the equation (already for α = 4 we end up with an integral of the error function). This makes it trickier to compute the leading correction in the distribution of eigenvalues, which we leave to future investigations.

B Alternative averages
In this appendix we consider nonuniform time averages where ℘ t is a probability distribution in [0, t]. The asymptotic behaviour of the moments of ρ t 0 ,t can be carried out as in the uniform case. Specifically, we have Let us change variables into where we rescaled (back) τ α because it does not appear in the gaussian anymore; we find Since the gaussian forces all the variables τ j with j ∈ 1, . . . , α − 1 to be O(1), we can extend their integration domain to infinity; in addition, at the leading order, they also disappear from the argument of ℘ t From this, we readily obtain the Rényi entropies and the von Neumann entropy We note that the von Neumann entropy is maximal when ℘ t (τ ) is uniform (℘ t (τ ) = 1 t ), suggesting (as expected) that the effective dimension D ( t) t is maximised by the uniform average. The distribution of eigenvalues can be computed with the method reported in the main text, and we find We point out that the change of scale in the eigenvalues brings the distribution into a universal form Finally, the dimension D ( t) t of the "weighted" space visited by the state is the solution to the following system (38)

C Numerical checks
In this appending we report some numerical checks of our findings in spin chains (d = 1).
Transverse field Ising chain. The transverse field Ising chain is described by the Hamiltonian where σ α ≡ · · ·⊗I −2 ⊗I −1 ⊗σ α ⊗I +1 ⊗I +2 ⊗· · · acts like the Pauli matrix σ α (α ∈ {x, y, z}) on site and like the identity elsewhere. A Jordan-Wigner transformation maps the spin chain into a chain of fermions, and the resulting Hamiltonian consists of two sectors where it acts like a quadratic form. This allows for exact diagonalization, and also for the exact solution of the dynamics when the initial state is a Slater determinant. For example, one can easily compute the function f (t) defined in (7) after a global quench of h : h i → h f (the system is prepared in the ground state of H(h i ) and then let to evolve under H(h f )) [9] f (t) = π 0 dk 2π log 1 + cos ∆ k 2 where We are therefore in a position to check our predictions for the first Rényi entropies ofρ t in the limit of large L. Figure 3 shows a comparison between the numerical evaluation of the entropies 8 (as logarithms of multidimensional integrals) and our asymptotic predictions. The agreement is excellent.
Exact diagonalization. Despite giving access to the first Rényi entropies, the exactly solvable model considered in the previous paragraph does not allow us to easily check the dimension of the relevant subspace. To overcome this problem, we have carried out a numerical analysis based on exact diagonalization algorithms in small spin chains with rather generic Hamiltonians. In principle our prediction is not expected to be accurate, as it was derived in the opposite limit of large L; the agreement between numerical data and prediction is nevertheless surprisingly good, as shown in figures 4 and 6.
We also check that our Ansatz t ∼ δt δt /t generates a subspace approximating the time evolving state with a given accuracy. To that aim, we have computed  Figure 4: The minimal number of eigenstates per unit √ L with probability larger or equal to 1 − t , with t = 0.15 / √ 1+100Jt in a small spin-1 2 chain with L = 6, 8, 10, 12 sites, obtained using exact diagonalization techniques. The initial state is the ground state of the ferromagnetic Ising Hamiltonian H 0 = −J (σ x σ x +1 + 2σ y ); time evolution is generated by H = J σ y σ y +1 + 0.5σ x σ x +1 + 1.5σ z σ z +1 + 0.25σ x + 0.3(−1) σ z . The dashed line is the prediction (19). Data seem to collapse to the prediction rather quickly.  Figure 6: The same as in figure 4 but for a different system. Here the initial state is fully polarised in the z direction and time evolution is generated by H = J σ x σ x +1 +2σ y σ y +1 + σ z σ z +1 , which describes an integrable system. which is the error made by projecting the state at the time t onto the subspace corresponding to T . As shown in figures 5 and 7 , the numerical data confirm the validity of the Ansatz. In all the examples that we considered, the maximal error is associated with the boundaries of the interval, while the error in the bulk of the time window is much smaller and approximately time independent.