{'fraction_non_alphanumeric': 0.04783940415592428, 'fraction_numerical': 0.010400828540578656, 'mean_word_length': 3.9810119635605314, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Unitary equilibrations: probability distribution of the Loschmidt echo

Closed quantum systems evolve unitarily and therefore cannot converge in a strong sense to an equilibrium state starting out from a generic pure state. Nevertheless for large system size one observes temporal typicality. Namely, for the overwhelming majority of the time instants, the statistics of observables is practically indistinguishable from an effective equilibrium one. In this paper we consider the Loschmidt echo (LE) to study this sort of unitary equilibration after a quench. We draw several conclusions on general grounds and on the basis of an exactly-solvable example of a quasi-free system. In particular we focus on the whole probability distribution of observing a given value of the LE after waiting a long time. Depending on the interplay between the initial state and the quench Hamiltonian, we find different regimes reflecting different equilibration dynamics. When the perturbation is small and the system is away from criticality the probability distribution is Gaussian. However close to criticality the distribution function approaches a double peaked,"batman-hood"shaped, universal form.


I. INTRODUCTION
A sudden change of the parameters governing the evolution of a closed quantum many-body system gives typically rise to a complex and fascinating dynamics. This so-called Hamiltonian quench is attracting an increasing amount of attention [1,2,3,4]. The reason for such an interest is, at least, twofold; in the first place this out of equilibrium phenomenon has been recently observed in cold atom systems [5,6]. Secondly, at a more conceptual level, the equilibration dynamics of a quenched quantum system plays a role in the very foundations of statistical mechanics [7,8,9,10,11,12]. New insights can be gained on the fundamental question about the emergence of a thermal behavior in closed quantum systems.
In this paper we study a prototypical dynamical quantity for a quantum quench: the Loschmidt echo (LE). This quantity is defined by the square modulus of the scalar product, of the time-evolved, (out of equilibrium) quantum state with the initial (equilibrium) one e.g., Hamiltonian ground-state.
In spite of the simplicity of its definition the Loschmidt echo L, or closely related quantities, convey a great deal of information in a variety of physical problems; for example L has been intensively studied in the context of Fermi edge singularities in x-ray spectra of metals [13], quantum chaos [14,15], decoherence [16,17,18], and more recently quantum criticality [19] and out-of-equilibrium fluctuations [20,21].
Typically (cfr. figure 1) the Loschmidt echo rapidly decays from its maximum value L = 1 at t = 0 and, after an initial transient starts oscillating erraticaly around the same well defined value. For finite size systems after a sufficiently long time a pattern of collapses and revivals is observed due to the almost-periodic nature of the underlying quantum dynamics. On the other hand for infinite volume systems the Hamiltonian spectrum generically becomes continuous and an asymptotic value L ∞ (coinciding with the average one) is eventually reached.
The main goal here is to investigate the statistical properties of the Loschmidt echo seen as a random variable over the ob-servation time interval [0, T ]. One of the key properties is that a small variance, by standard probability theory arguments, guarantees that the overwhelming majority of the time L(t) sticks very close to its average value [7,10,12]. This is the sense in which one can speak about equilibration dynamics and corresponding "equilibrium properties" of a finite system that is evolving unitarily and therefore cannot have attractive fixed points.
We shall show how the features of the probability distribution of the Loschmidt echo depend on a rich interplay between the initial state and quench Hamiltonian on the one hand and the system's size and observation time on the other. In particular we will focus on the potential role that the vicinity of quantum critical points may have on the features of the Loschmidt echo probability distribution function [4,17,19]. This latter analysis will be mostly carried over by exploiting exact results for quasi-free spin chain i.e., the quantum Ising model [19].
The paper is organized as follows: in Section II we give the general setting. Later we introduce a relaxation time T R and discuss the universality content of the L (t) before this timescale. In Section II C we define and study other relevant time scales, the time T 1 for necessary for observing the correct average, and revival times where large portion of L (t) are back in phase. In section II D we give explicit formulas for the moments of the LE assuming the non-resonant hypothesis. In section III we concentrate on a particular example and prove all the general results advocated so far for an exactly solvable case. Moreover we discover three universal behaviors for the whole LE probability distribution function. We draw some parallels with another natural quenched observable: the magnetization. Finally section V is devoted to conclusions and outlook.

II. GENERAL BEHAVIOR
Let us start by recalling a few elementary yet crucial facts. If H = n E n Π n is the system's Hamiltonian (Π n 's=spectral projections) the closed-system dynamics is described by the time-evolution superoperator U t = e −itH , H(X) = [H, X]. This superoperator is thought here of as a map of the space of trace-class operators X into itself ( X 1 := tr √ X † X < ∞). Closed quantum systems cannot equilibrate in the strong sense, as unitary evolutions U t do not have non-trivial i.e., non fixed, limit points in the norm topology for t → ∞ [32]. One may then wonder whether a weaker form of convergence can be achieved for t → ∞.
Let us then consider the expectation value of an observable A(t) := tr (U t (ρ 0 )A) and write the spectral resolution of the superoperator U t as a formal sum U t = E e −itE |E E|, here E (|E )denote the eigenvalues (eigenvector) of H. In finite dimensions the kernel of H is spanned by the Π n 's and gives rise to a time-independent contribution to A(t) i.e., A ∞ := n tr(Π n ρ 0 Π n A); the point is now to understand whether the remaining components involving the non-trivial time-dependent factors exp(−iEt) (E 0) admits a limit for t → ∞. In finite dimensions the E are a (finite) discrete set of differences of Hamiltonian eigenvalues e.g., E = E n − E m , and correspondingly A(t) − A ∞ is a quasi-periodic function: the long time limit of A(t) does not exist. On the other hand in the infinite dimensional case the spectrum of H can be continuous and, if the functionÂ(E) := ρ 0 , E E, A is sufficiently well behaved, using the Reimann-Lebesgue lemma, lim t→∞Â (E) exp(−iEt) = 0. Therefore in this case denotes the time-average over an infinite time interval. While this convergence cannot be achieved uniformly for the whole set of system's observables (in that it would imply strong convergence) it can be proven for specific families of A's e.g., local ones [22].
There is a third form of convergence that one can consider here: the convergence in probability. In the following we shall consider the above defined A(t) as a random variable over the real line of t endowed with the uniform measure dt/T with T → ∞. Suppose we have a sequence of {B L } L (think of L as the system size), we say that the B L 's converge to zero in probability if The meaning of this type of convergence should be clear: for large L the probability of observing a value of B L (t) different from zero is vanishingly small. In other words the fractions of t's for which B L (t) 0 is going to zero for L → ∞. As a matter of fact this is the type of convergence, with B L (t) = A L (t) − A L,∞ , that has been considered in [10,12]. The stochastic convergence (2) implies that the probability distributions of the random variables A L i.e., P L (α) := δ(α − A L (t)) are converging to the one of A ∞ := lim L→∞ A L,∞ i.e., lim L→∞ P L (α) = δ(α − A ∞ ). In this context we say that the initial state ρ 0 "relaxes" or "equilibrates" to ρ eq if it happens that A ∞ = tr(Aρ eq ). [33].
A typical strategy to demonstrate this kind of unitary equilibration is to prove (2) Figure 1: Typical behavior of the Loschmidt echo for the Ising model in transverse field. All curves refer to a size of L = 100. In the upper left panel the relaxation time T R is indicated. Using arguments as in section III D (see also notes [23,24]) one is able to show that, in the quantum Ising model at criticality, the first revival time is exactly var(B L ) goes to zero sufficiently fast for L → ∞. If this is the case one can use a basic probability theory result Pr{t Since, under the assumption ii) the RHS of this latter relation can be made arbitrarily small and given i), equation (2) holds true.
Notice that yet another way to formulate the kind of convergence (2) is by means of the concept of typicality [8,9]; the probability of observing a "non-typical" value i.e. one that deviates significantly from the mean one becomes negligible in the large L limit.

A. The Loschmidt echo
The time dependent quantity we are going to focus on in the rest of this paper is the Loschmidt echo (LE): where the state |ψ is possibly, but not necessarily, the ground state of the Hamiltonian H at a different coupling. In the sequel, statistical averages are always taken with respect to this state. In the following we will consider L (t) as a random variable with uniformly distributed t ≥ 0. Ideally we are interested not only in the first moment but in the whole probability distribution function. The probability of L to have value in Ω is given by For those x for which the probability density is well defined, it is given by .
The k-th moment of this probability distribution is given by µ k := x k P(x)dx = L k (t). Notice that the Loschmidt echo can be written as L(t) = ρ ψ , e −iHt (ρ ψ ) , where: ρ ψ := |ψ ψ|, H(X) := [H, X] and X, Y := tr(X † Y) denotes the Hilbert-Schmidt scalar product. From this it follows Performing the time average [34] one finds µ n = ρ ⊗ n ψ , P (n) (ρ ⊗ n ψ ) , where P (n) projects onto the kernel of H (n) . In particular the time average L = µ 1 is given by where ρ eq := P (1) (ρ ψ ). From the general discussion in Sect.

B. Short time regime and criticality
As already pointed out, typically the LE decays from its maximum value 1 at t = 0 and, after an initial transient, starts oscillating erraticaly around its mean value. In this section we will analyze the universality content of this initial transient and its dependence on the interplay between the initial state |ψ and the evolving Hamiltonian H.
We start by noticing that the LE Eq. (3) is the square modulus of a characteristic function χ (t) = e −itH which is the Fourier transform of the energy probability distribution: . Both χ and the LE can be expressed in terms of the cumulants of H: where · c stands for the connected average with respect to |ψ . The sums above starts from n = 1 because the zero order cumulant is zero: H 0 c = 0. Since H is a local operator, i.e. a sum of local "variables", we can expect in some circumstances, the central limit theorem (CLT) to apply. More specifically the version of the CLT we are going to consider here, is the following. In the thermodynamic limit, the probability distribution of the rescaled variable Y ≡ (H − H ) / H 2 c tends to a Gaussian (with variance 1 and mean zero) . In other words, all but the second connected moments of Y tend to zero when the volume goes to infinity.
When the CLT applies, for sufficiently large system sizes, the distribution of H will be of the form One can systematically compute corrections to this formula and order them as inverse powers of the system size L. Fourier transforming back we obtain the characteristic function and the LE It may seem that this expression for the LE could have be obtained right away by keeping only the first term in the expansion of the exponential in Eq. (6): This is simply a quadratic approximation, that does not rely on the CLT. Its validity requires t 2 σ 2 ≪ 1 that in turn implies t ≪ 1/ √ H 2 (as will be explained later this means roughly t ≪ L −d/2 where d is the space dimension). From figure 1 we see that typically L (t) decays from 1 and after an initial transient, starts oscillating around an average value L which will be computed below. Now, when the CLT applies, equation (7) can help us to define this transient or relaxation time T R given by e −T 2 R σ 2 = L. Roughly after this time one starts seeing oscillations in L (t). Since in general (see below) one has L ≃ e − f L d , and when the CLT applies σ 2 scales like the system volume, the relaxation time scales as The situation is different when one considers small variations of the parameters δh. In this limit the average LE is related to the well studied ground state fidelity F = ψ (1) |ψ (2) . More precisely one has L ≃ F 4 [18]. Close to critical points the behavior of the fidelity is dictated by the scaling dimension ∆ of the most relevant operator in H with respect to the critical state |ψ [25]. The precise scaling is the following: F ∼ 1 − const. × δh 2 L 2(d+ζ−∆) , where ζ is the dynamical critical exponent. Similarly one can show (see below) that σ 2 ∼ L 2(d−∆) . All in all this amounts to saying that the relaxation time for small variation δh around a critical point (roughly δh ≪ L −(d+ζ−∆) ) increases from O (1) to T R ∼ L ζ . In the thermodynamic limit instead, i.e. taking first the limit L → ∞, using standard scaling arguments, one can show that for δh small and h close to the critical point h c , the relaxation times diverges as T R ∼ |h − h c | −ζν , ν being the correlation length exponent. In figure 2 we plot the relaxation time for a concrete example that will be studied thoroughly in section III, the Ising model in transverse field. There one has ∆ = ζ = ν = 1, and so the singularities observed in figure 2 are of simple algebraic type: T R ∼ |h − h c | −1 around the critical points h c = ±1.
Let us now turn to discuss when we expect the CLT to work. Consider first the case when |ψ is the ground state of a gapped Hamiltonian and the connected energy correlators go to zero . In this case the connected averages of H scale as the volume: H n c ∼ L d . As a consequence the cumulants of the rescaled variable satisfy Y n c ∼ L −(nd−2d)/2 for n ≥ 2, which immediately implies the CLT in the sense given above.
Therefore we can have violation of the CLT only in the gapless case when the state |ψ is critical or when clustering fails. Let us then consider a critical state |ψ . Connected averages have a regular extensive part and a singular part which scales according to the most relevant component of H with scaling dimension ∆. When |ψ is the ground state of H at a different coupling, ∆ is the scaling dimension of the perturbation δH to the critical Hamiltonian. At leading order, we can write If ∆ < d/2 the cumulants of the rescaled variable Y don't go to zero but to universal constants [26,27] given by In this case the probability distribution of the energy is a, non-Gaussian, universal distribution. This kind of universal behavior has been observed for instance in [28] on an example where the scaling dimension is ∆ = 1/8. In the opposite situation where ∆ > d/2, all the cumulants of Y go to zero except for the first two, and the distribution function approaches a Gaussian in the large size limit. In the intermediate case ∆ = d/2 the cumulants of Y do not go zero but to a constant which is however not universal due to extensive contributions coming from the denominator. We recall that in d dimensional, zero temperature quantum mechanics, operators are classified into relevant, irrelevant, and marginal if their scaling dimension is respectively smaller, larger, or equal to d + ζ where ζ is the dynamical exponent. Hence we see that, even in the critical case, we observe deviation from the Gaussian behavior only if the perturbation δH is sufficiently relevant, specifically ∆ < d/2.
To finish let us remind the reader that the CLT also breaks down when clustering fails.
To summarize, when the CLT applies, the LE tends to a Gaussian and plotting the function L L t/ H 2 2 for different sizes L one should observe data collapse (see figure 3).
The same data collapse feature is observed when choosing different h (i) s, although the variance of this Gaussian is sensitive to that.

C. Equilibration and long time behavior
After having discussed the short time behavior of the LE related to the initial transient, let us now turn to its long time behavior.
We first re-write Eq. (3) in the eigenbasis of H = n E n |n n| where p n = | ψ|n | 2 . If the spectrum of H is non degenerate the superoperator P (1) acts as a dephasing in the Hamiltonian eigenbasis i.e, P (1) (X) = n n|X|n |n n|. In other words the time average of the exponentials in Eq. (8) gives simply δ n,m and equation (4) reduces to L = n p 2 n . As is well known this quantity is the purity of an equilibrium, dephased, state: ρ eq = n p n |n n|.
Time scales In the preceding section we already defined a relevant time scale, the relaxation time T R which is O (1) off-criticality while T R = O L ζ in the critical case and for sufficiently small variations δh ≪ L −(d+ζ−∆) .
In some situations it is useful to consider a finite observation time T . We will write L to indicate the corresponding average. It is natural to ask about the interplay between the observation time T and the linear size of the system L. In other words in general lim L→∞ lim T →∞ L L lim T →∞ lim L→∞ L L . Since taking larger system sizes has the effect of sending the revival times to infinity and the LE attunes its maximum value 1 at t = 0, typically the function lim L→∞ L L (t) has only one large peak at t = 0 whereas L L (t) has peaks at all the revival times. Correspondingly we expect lim L→∞ lim T →∞ L L > lim T →∞ lim L→∞ L L . This expectation has been confirmed for the case of the one dimensional quantum Ising model, see section III.
Another question which is relevant in the measurement process is how large must the observation time be to effec-tively measure L? That is, what is the condition to have L T 1 = L or more in general what is the smallest time T n such that one observes (L) n T n = (L) n ? Let us focus on T 1 . Looking at equation (8) one realizes that it suffices to have T ≫ ∆ −1 min , where ∆ min is the smallest gap in the whole spectrum i.e. ∆ min = min n,m (E n − E m ). We can address this question for the class of quasi-free Fermi systems in d spatial dimensions. In this case the energy has the form E n = k n k Λ k where k is a d-dimensional quasi-momentum-like label and we can assume the one particle energy Λ k to be positive. Then the gap is given by ∆ min = min β k k β k Λ k with β k = 0, ±1. By choosing β k = (−1) k , one obtains a ∆ min which is exponentially small in L in those (frequent) cases where Λ k is an analytic function of k [29]. However in quasi free systems the weights p n decrease exponentially with the number of excitations in n. In practice the highest weight is given for energy differences between the one and zero particle spectra: ∆ (1,0) = min k Λ k which is a constant of order 1 in the gapful case, while typically scales as L −1 for the critical case. The next largest amount of spectral weight is attained at a gap which is a difference between one particle energies ∆ (1,1) = min k,q Λ k − Λ q . It will be favorable to have k and q nearby in the region where Λ k is flat or almost flat. So we get . From this discussion we estimate that, at least in quasi-free systems, to have L T 1 ≃ L one must take T 1 = O (1) in the gapful case, while one has T 1 = O (L) at criticality. If one needs L T 1 ≃ L with a larger degree of precision than one must choose considerably larger time: 0 has a solution within the allowed set of k-vectors). Related time scales are revival times. We define a revival time to be that particular time for which a large portion of spectral weight p n p m has revived. More precisely T rev ω peak = 2π, where ω peak is a particular frequency E n − E m such that the weight p n p m is large. From the discussion above we expect T rev = O (1) when H has a gap above the ground state, while T rev = O (L) when H is critical. These expectations have been confirmed (see figure  1) on the hand of a solvable model that will be discussed in the next sections.

D. Moments of the Loschmidt echo
Having computed the time averaged LE we can now turn to higher moments. In doing this one has to distinguish cases where n = m from those where E n = E m in Eq. (8). So we write and the n-th moment is given by The computation of the average [X (t)] k can be done assuming a strong non-resonance condition. With this we mean the following. We say H satisfies a k-non-resonance condition if the only way to fulfill k l=1 E i l − E j l = 0 is to match the E i 's to the E j 's. Strong non resonance is k-non-resonance for any k. Note that this condition cannot be fulfilled when k becomes of the order of the Hilbert's space dimension. Now to compute [X (t)] k draw 2k points in two rows of length k. Imagine For [X (t)] 3  Both these diagrams give the same contribution (simply swap i's with j's) and the result is The number of terms in [X (t)] k , N (k) is the number of all permutations without fixed points and is given by However, among these N (k) terms, many of them give different contributions. Look for instance at [X (t)] 4 : Correctly one has 3 + 6 = N (4) = 9.
We collect here the first three moments while the cumulants are We can notice that each term in [X (t)] k has the same form of Later we will encounter situations where this function gives a good approximation to the Loschmidt echo probability distribution.

III. ISING MODEL IN TRANSVERSE FIELD
From now on we will give a detailed description of the Loschmidt echo for the case of an exactly solvable model. The model we consider is the Ising model in transverse field with Hamiltonian This model can be mapped to quasi-free fermions and so diagonalized exactly. At zero temperature we distinguish two phases: i) An ordered one in the longitudinal direction in for |h| < 1, and ii) A paramagnetic phase for |h| > 1. The points |h| = 1 are critical points where the system is described by a conformal invariant field theory with central charge c = 1/2.
Exploiting the fact that H decomposes into a direct sum of L/2 blocks 4 × 4, we are able to compute the complete dephased equilibrium state ρ eq = n p n |n n|. The result is where the multi-index α is α = (α 1 , . . . , α L ), Finally the weights are given by Using Eq. (10) together with (11) one can show (cfr. ref. [31]) that the dephased state has the following totally factorized form:

A. Short time regime
Let us first discuss the short-time, transient regime. Looking at the function ln L (t) one can readily see that all its nderivatives at t = 0 are the Riemann sums of a summable function irrespective of the h (i) s being critical. This means that all the derivatives of ln L (t) grow linearly with L, and together with equation (6) implies that the cumulants are linear even at criticality, i.e. H 2n c ∝ L. The same result could have been derived by noting that for |h| 1 the system is gapful and clustering. When |ψ (1) is critical, i.e. h (1) = 1, the scaling dimension of δH = H (2) − H (1) = −δh i σ z i is one and so according to the reasoning in section II B the cumulants of H grows linearly with L.
Accordingly the CLT applies. Since all the cumulants, including the variance, grow as L, plotting the function L L t/ √ L for different sizes L, one observes data collapse. This behavior is illustrated in figure 3. Clearly the plot reproduces a Gaussian with variance lim L→∞ H 2 c /L.

B. Long time, large sizes and the order of limits
Consider now a physical situation where an experimenter computes L L . We inserted the labels T and L to stress that both size and expectation time are finite, as is required in a true experiment. Here we want to study the interplay between T and L. Consider first the case where we send L to infinity, or more physically L is the largest scale of our system. In this situation the spectrum of H is practically continuous, and we can write the LE as In this approximation we sent all the revival times to infinity and so the function L (t) is no longer almost periodic, but rather tends to a precise limit as t → ∞. We can calculate the limit s (∞) and also the first corrections, as t → ∞. The procedure is outlined in the Appendix. The result is where s (∞) is the limiting value and A m/M are constants which depend on h (i) and are given in the appendix. The result (26) has already been found in [4], here we provide the explicit form of the asymptotic value, s (∞).
Since the function L (t) = e −Ls(t) has a limit at infinity, its time average is precisely this limit L = e −Ls(∞) . More precisely the distribution function becomes a delta function P (x) = δ x − e −Ls(∞) .
Consider now performing first the time average of Eq. (9). According to the discussion in section II C, this requires at least observation times as large as T ≫ L (if we are at criticality). The result in this case is (for the explicit computation see the following section) If now L is large, we can approximate the sum with the integral: L = e −Lg(h (1) ,h (2) ) . Calling δϑ k = ϑ (2) k − ϑ (1) k the two functions, g and s (∞) are given by We observed that the two averages e −Lg -obtained by first doing the time average and then taking large L -or e −Ls(∞)obtained by first considering L large and then doing the time average -are qualitatively very similar for most values of the parameters h (i) . The only region where there is an appreciable difference is when h (1) and h (2) correspond to different phases (either h (1) < 1 and h (2) > 1 or vice-versa).

C. Moments of the Loschmidt echo in presence of degeneracy
Having the explicit form of the LE Eq. (9) we can compute its time average and also other moments. Since the Ising model is mapped to a free Fermi system on a finite lattice, and given the form of the quasi-particle dispersion Λ k , its spectrum is non-degenerate. In other words the Ising Hamiltonian is 1-non-degenerate. However, for the same reason, it is not k-non-degenerate for k ≥ 2. This means that to compute moments higher than the first, we really need to use the explicit form Eq. (9) and cannot rely on the results of section II C. For the first moment this problem does not arise, and we can either use L = n p 2 n or do the time average of Eq. (9). Correctly the results coincide, and they rely on the fact that k (n k − m k ) Λ k = 0, implies n k = m k , i.e. that the spectrum is non-degenerate. The result is For later convenience we define α k = sin 2 ϑ (1) k − ϑ (2) k . To compute higher moments we first rewrite Eq. (9) as We write the n-th power of the LE as Now, when computing M products of Y (n) k terms, only the γ = 0 term will survive after taking the time average. In other words: 0,k M and so we have An explicit formula for g (n) 0,k is Noting that For example we have The variance of the LE is then given by The first moment and the variance are plotted in figure 4. One should note that close to the critical points h c = ±1 there appears a small region δh where the variance is large.
Equation (17) gives explicitly the variance in a case where the non-resonant hypothesis is violated. Since |α k | ≤ 1 generally the variance is given by the difference between two exponentially small quantities and so, a fortiori, is exponentially small in the system size L. However, looking at figure 4 one notes a small region of parameter close to the critical points, where the variance is large. As we will see this fact has important consequences.
It is now interesting to compare the result in equation (17) with what would have been obtained assuming a non-resonant condition. Clearly the second moment computed assuming non resonance has less terms than the correct one. Since for the LE all the contributions are positive, the non-resonant result ought to be smaller. In other words, for the variance we must have ∆L 2 ≥ ∆L 2 nr . Comparison with the non-resonant result To compute the variance assuming non-resonance we use Using this formula together with Eq. (11) we obtain We have verified that the inequality ∆L 2 ≥ ∆L 2 nr holds for all values of the coupling constants h (1) and h (2) . However the qualitative behavior of ∆L 2 nr is very similar to the true variance ∆L 2 .

D. The Loschmidt echo distribution function
We now turn to consider the whole probability distribution of the LE. As we have noted earlier, for any L finite being the spectrum discrete, L (t) is an almost periodic function. Actually L (t) belongs to a smaller class, since it is a trigonometric polynomial. In any case, most results we will present are valid for the larger class of almost periodic functions.
We now give the results for the whole LE probability distribution function. We have observed three kinds of universal behavior emerging in different, well defined regimes. i) Exponential behavior where the probability distribution is well approximated by ϑ (x) e −x/L /L (figure 5), ii) Gaussian behavior, (figure 6), and iii) A universal double peaked, "batmanhood" shaped function, (figure 7). More precisely we have the following scenario: • δh large. In this case, for L moderately large, the distribution is approximately exponential. The feature is more pronounced when h (1) , h (2) are in different phases, the limiting case being h (1) ≈ ±h (2) .
In other words, say that we fixed h (i) in order to have either a Gaussian or a double-peaked distribution. We can always find an L large enough such that the distribution becomes exponential in both cases. However, for the Gaussian case we must reach considerably larger sizes L ≫ |δh| −2 (compare figures 8 and 9). We have observed an exponential distribution in the region of parameters where the average LE is much smaller than one: L ≪ 1. Due to the bound L 2 < 2L 2 , one has ∆L < L 2 so that when L ≪ 1 even the variance is small. Since the LE is supported in [0, 1] and in particular L (t) must be positive we expect in the region L ≪ 1 a distribution with positive support, with a large peak very close to zero, and rapidly decaying tail. We have verified that an exponential distribution of the form ϑ (x) e −x/L /L gives a pretty good approximation in the region L ≪ 1. Note in any case, that the exponential form is always an approximation. In particular, for x → 0 the true distribution P (x) tends to zero for any value of the parameters. This happens since we have always L (t) > 0 strictly. And so generally 0 < L min ≤ L ≤ 1. This feature can be accounted for by adding a (small) ǫ term to the exponential: Let us now investigate the conditions under which the first moment is small and so we expect an approximately exponential behavior. Looking at equation (13) we see that L ≪ 1 holds when Lg h (1) , h (2) ≫ 1 where g is given by Eq. (14). Clearly the function g is zero (its minimum) when h (1) = h (2) , and is quadratic in the difference on the diagonal. Quite interestingly the function g is appreciably different from zero only when h (1) and h (2) correspond to different phases (i.e. either h (1) < 1 and h (2) > 1 or vice-versa) so that in these cases we observe exponential behavior even for moderately small lattices. When h (1) is close to h (2) , but away from critical points, the function g is quadratic in the difference δh so that L ≈ exp −const. × Lδh 2 . Hence to have L ≪ 1 and so to observe approximate exponential behavior we obtain the relation L ≫ δh −2 . At criticality and for δh small instead, we can use L ≈ F 4 where F is the fidelity which scales as (see section II B and ref. [25]) . In the quantum Ising model we are considering we have d = ζ = ν = ∆ = 1 and so the average behaves as L ≈ exp −const. × L 2 δh 2 . This means that for δh small around a critical point h c = ±1, the condition to have an exponential distribution becomes L ≫ δh −1 .
We now turn to consider the origin of the batman-hood shaped distribution function. As we have already noticed, the LE is a (finite) sum of cosines with given frequencies and amplitude. We can imagine a situation where only few frequencies contribute to the LE. In the limiting case, only two non-zero terms. That means that the LE can be approximated by where we can assume A, B positive. We have devoted some attention to the probability distribution generated by such a function. If ω A and ω B are rationally dependent, the function is periodic and the distribution function has square root singularities at all values of L where ∂ t L (t) = 0. However in our case the frequencies ω A/B are always rationally independent. In this case the vector x (t) = (ω A t, ω B t) wraps around the torus in a uniform way. We can then invoke ergodicity and transform the time average into a "phase space" average (in this case the phase space is x = (x 1 , x 2 )). Hence the probability distribution function is given by By using Eq. (18) this probability density can be written as The integral above can be expressed in terms of elliptic functions, but we won't need its explicit expression. Typically the function (19) is batman-hood shaped function (see fig. 7), with support in L − |A + B| , L + |A + B| and two peaks at χ = L ± |A − B|. The divergence at the peaks position is of logarithmic type, close to the peaks (i.e. χ = L ± |A − B| + ǫ) one has Note that we never observe the limiting case where B = 0 and L (t) becomes a periodic function. This means that even a very small spectral weight on B cannot be discarded. On the other hand the distribution function (19) seems to be quite stable against the presence of other oscillating terms with small spectral weight. This stability can be seen in fig. 7 where one clearly has at least three frequencies with reasonable spectral weight, but the probability density is still well approximated by a batman-hood. Spectral analysis To understand the behavior of P (L = x) we do a spectral analysis of L (t) to see which frequencies contribute most. In fact, for almost periodic functions there is a similar Fourier decomposition as for periodic functions. The Fourier expansion is given in this case byL disc (ω) = L (t) e iωt . Taking into account Eq. (8)L disc (ω) can be written aŝ We would like to know which frequencies have the largest weight. This is achieved by expanding the product in Eq. (9) [35] Now, since eachX k (t) is smaller than 1/2 in modulus, it is reasonable to approximate the LE with the first two terms of this expansion and we obtain In this approximation we only wrote the zero-frequency contribution, which corresponds to the mean, and the contribution coming from the one particle spectrum. The next term has also contributions coming from the two particle spectrum.
To be more precise, call E (n) the energy of the n-particle spectrum then a − E (0) b ∝ Λ k 1 + Λ k 2 (second order contribution with less spectral weight).
Note that we expect Eq. (21) to be approximately valid (with a different form for the amplitudes and the frequencies)  16) and (17). In the bottom panel one can notice that many frequencies contribute to the LE. The one particle contribution c (ω), Eq. (22), is given by the black dots while the red curve gives the true spectral decompositionL disc (ω) given by Eq. (20).
also for non integrable models in which a one-particle approximation works well. Now, if there is a regime where the amplitudes α k /2 are highly peaked around few quasi-momenta, in the limiting case only two, then the LE can be approximated as in Eq. (18) and we expect a double peaked distribution function. So we are led to study the (one-particle) amplitude function Generally c (ω) is a bell-shaped function, starting linearly from the band minimum E m , reaching a maximum value and then decreasing to zero at the band maximum E M . It is not difficult to show that [23], when δh is small and for roughly 1 − h (2) 10 −1 , c (ω) starts developing a peak, the width of which being proportional to 1 − h (2) . In the limiting case h (2) = 1, c (ω) has its maximum at E m = 0 and then decreases monotonically to E M = 2. So, for δh small and for 1 − h (2) 10 −1 , c (ω) is a peaked function. In order to have few frequencies fall within the peak, and so to have large spectral weight on few frequencies, we must additionally have L 1 − h (2) ≪ 1. This is easily seen analyzing the dispersion Λ k for h (2) close to the critical point [24]. All in all, the conditions to have a batman-hood distribution, are δh small and L 1 − h (2) ≪ 1. The feature is more pronounced when the h i are not precisely critical. In fact even though c (ω) is most  (21)). Lower panel: the red curve is the Fourier seriesL disc (ω) (projected spectral density), black curve shows the highest coefficient c (ω) given by Eq. (22), together with allowed frequencies (black dots).
peaked when h (2) = 1 (and the peak is at ω = 0), we have to remember that the allowed values of ω are ω n = Λ k n where k n = π (2n + 1) /L, and the smallest frequency is ω 1 = Λ π/L . If we perturb h (2) from 1 the peak of c (ω) shifts to the right, approaching ω 1 , so it is favorable to have h (2) 1. Not surprisingly, the conditions to have a batman-hood probability density, coincide with having a large variance (see figure 4 bottom panel). Generally fixing h (i) and increasing the size L, one eventually violates the quasi-critical condition L ≪ ξ. At this stage the double-peak feature tends to disappear and the distribution approaches an exponential one. This can be clearly seen in figure 8. From this figure one can have the impression that the "double-peak feature" is a prerequisite of short sizes, since in this case one has few frequencies anyway. As we have tried to explain instead, this feature survives for larger sizes, provided we shrink δh sufficiently ( figure 7).
Instead, when δh is small but h (i) are far from the critical point, then c (ω) is not peaked, and many frequencies have a large spectral weight (see figure 6). In this case the distribution becomes Gaussian.The emergence of a Gaussian distribution can be qualitatively understood in the following way. First write the LE according to its spectral decomposition L (t) = n A n e itω n where the amplitudes are precisely given by A n =L disc (ω n ) and are positive. When the frequencies are rationally independent the variables x n = tω n wrap uniformly around a large dimensional torus. Then one can consider each A n × e tx n as an independent random variable. The assumption δh is small but h (i) away from criticality corresponds to say that L (t) can be considered as a sum of many independent random variables, giving rise to a Gaussian distribution. When L ≫ |δh| −2 the conditions of independence breaks down and we recover an approximate exponential behavior.

IV. PROBABILITY DISTRIBUTION FUNCTION FOR THE MAGNETIZATION
In the same spirit we can compute the probability distribution of a local operator. The first candidate that comes to mind is the transverse magnetization. We computed the following time dependent observable m (t) = ψ (1) |e itH (2) σ z i e −itH (2) |ψ (1) . Using again equation (10) where the quasi-momenta range now in the whole Brillouine zone: k = π (2n + 1) /L, n = 0, 1, . . . , L − 1. Correctly, when h (1) = h (2) we recover the zero temperature equilibrium result σ z i = L −1 k cos (ϑ k ). From equation (23) we see that m (t) can be writtenexactly-as a constant term plus an an oscillating part with frequencies given by the single particle spectrum Λ k . The discussion on characteristic times becomes simplified as all time scales are uniquely determined by Λ (2) k . For example the time T 1 necessary to observe the correct mean: m T 1 = m must simply satisfy T 1 ≫ gap −1 which means T 1 ≫ L in the quasi-critical regime h (2) − 1 −1 ≫ L, while it suffices to have T 1 ≫ O (1) away from criticality. Given equation (23) it is not difficult to compute the mean and the variance, which are given by Some comments are in order here. First fixing h (1) , h (2) the variance (figure 10) goes to zero as L −1 and not exponentially fast as was the case for the LE. Second, the spectral weight associated to the frequency Λ (2) k is sin ϑ (2) k sin (δϑ k ). We observed that there are always many frequencies with large spectral weight. In other words the spectral weight function is never peaked.
These comments suggest us to expect a Gaussian behavior for the probability distribution function of the magnetization irrespective of the parameters approaching critical values. This has indeed been observed ( fig. 11).

V. CONCLUSIONS
The unitary character of the dynamics of a closed quantum system implies that whatever relevant notion of equilibration one might have has to be a subtle one. In this paper we investigated the unitary equilibration of a quantum system after a sudden change of its Hamiltonian parameter. To this aim we used a prototypical time-dependent quantity: the Loschmidt echo. We established how the global features of L depend on the physical properties of the initial state preparation and Here below for the reader's sake we summarize the main findings of the paper • Resorting to a cumulant expansion we characterized the "short" time behavior of L(t). Different regimes can be identified depending on the most relevant scaling dimension of the quench Hamiltonian. When the central limit theorem (CLT) applies one has a Gaussian decay over a time scale O(1) for gapped systems. For the critical case the time-scale becomes O L ζ in a small region |δh| ≪ L −(d+ζ−∆) . At critical points the CLT can be violated and L(t) takes a universal non-Gaussian form for sufficiently relevant perturbations.
• We discussed the general structure of the higher momenta of P i.e., µ k := L k (t) = x k P(x)dx using the so-called non-resonant hypothesis. We showed that all the µ k are bounded by those corresponding to a Poissonian distribution i.e., P(x) = ϑ (x) exp(−x/L)/L.
• Using exact results for the quantum Ising chain we rigorously analyzed the interplay between the chain length L and the averaging time T. In particular we showed how the limit lim T →∞ and the thermodynamical one i.e., lim L→∞ do not commute. While in finite systems the L(t) is an almost-periodic function, in the thermodynamical limit L ∞ := lim t→∞ L(t) exists and P(x) → δ(x − L ∞ ). We explicitly computed L ∞ and the way it is asymptotically approached for large t. We gave a general closed form for the exact µ k 's and compared with that obtained with the non-resonant hypothesis • For the quantum Ising chain we numerically investigated P(x). We identified three universal regimes a) An exponential one (P is Poissonian) when L is the largest scale of the system b) A Gaussian one for intermediate L and initial state and quench parameters close and off-critical c) A "Batman-hood" shape for P when the parameters are close to each other and close to criticality. This result holds in the quasi-critical region L|h (i) − 1| ≪ 1, (i = 1, 2).
• Finally, for the sake of the comparison of the Loschmidt echo with a prototypical observable, we computed the time-dependent magnetization after the quench and studied its long-time statistics. In this case only a Gaussian regime appears to be reachable.
We have shown that the Loschmidt echo encodes sophisticated information about the quantum equilibration dynamics. For finite system vastly different time-scales arise: short time relaxation is intertwined with a complex a pattern of collapses and revivals and eventually Poincare recurrences. Unveiling how these phenomena depend on spectral properties of the underlying Hamiltonians, is one of the key challenges in the way to understand emergent thermal behavior in closed quantum systems.
To compute the first correction to the limit note that α (ω) k smoothen the singularity at the band edge, so that for the leading correction we need only k = 1 in the expansion of the logarithm. The evaluation of the oscillating integral is done with a saddle point technique. The result is with constants given by (we assumed here h (2) > 0 so that the band minimum is E m = 2 1 − h (2) ) 1 + h (1) 2 h (2) 3/2 1 + h (2) .
The result (26) should be the same as the square modulus of Eq. (12) in [4]. However in [4] there appears only one frequency, corresponding to the lowest band edge. The discrepancy probably arises from a continuum approximation which discards the effect of the van Hove singularity present at the upper band edge. As we have seen both terms give similar contributions. In particular, even at criticality, the van Hove singularity at the upper band edge survives.
In figure 12 one can appreciate the validity of the approximation (26).
Since F n T is a bounded operator its expectation value divided by T goes to zero when T → ∞.
[35] A more precise approximation valid also in region where δh is large, is given by L (t) = k>0 1 + c k 0 +X k (t) ≃ L 1 + k>0 1 + c k 0 −1X k (t) . The two expressions coincide when δh is small.