General Features of the Relaxation Dynamics of Interacting Quantum Systems

We study numerically and analytically isolated interacting quantum systems that are taken out of equilibrium instantaneously (quenched). The probability of finding the initial state in time, the so-called fidelity, decays fastest for systems described by full random matrices, where simultaneous many-body interactions are implied. In the realm of realistic systems with two-body interactions, the dynamics is slower and depends on the interplay between the initial state and the Hamiltonian characterizing the system. The fastest fidelity decay in this case is Gaussian and can persist until saturation. A simple general picture, in which the fidelity plays a central role, is also achieved for the short-time dynamics of few-body observables. It holds for initial states that are eigenstates of the observables. We also discuss the need to reassess analytical expressions that were previously proposed to describe the evolution of the Shannon entropy. Our analyses are mainly developed for initial states that can be prepared in experiments with cold atoms in optical lattices.


I. INTRODUCTION
Despite the ubiquity of many-body quantum systems out of equilibrium, they are much less understood than quantum systems in equilibrium. To advance our understanding and to construct a general picture, it is necessary to single out what is universal about their behavior. Identifying factors that determine how fast these systems evolve in time [1][2][3][4][5] is also essential for the development of algorithms for quantum optimal control [6]. In these two contexts, the unitary evolution of isolated many-body quantum systems is of particular interest due, in parts, to the connection with current experiments in optical lattices [7][8][9][10][11][12][13][14][15][16][17]. The latter are quasi-isolated systems, where coherent evolutions can be studied for very long times.
The evolution of an isolated system can be initiated by changing instantaneously the parameters of a certain initial Hamiltonian which is brought into a new final Hamiltonian. This abrupt perturbation is referred to as a quench. The fidelity [18,19], which is defined as the overlap between the initial state and its evolved counterpart, is a way to characterize the system evolution. This quantity is also analogous to the characteristic function of the probability distribution of work [20][21][22] and is therefore likely to find applications in quantum thermodynamics, particularly in studies related with the quantification of the work done to take quantum systems out of equilibrium. The fidelity decays exponentially when the final Hamiltonian is chaotic [23][24][25][26][27][28][29][30]. In fact, this behavior is expected to hold even in integrable Hamiltonians provided the initial state be sufficiently delocalized in the energy eigenbasis [31][32][33].
Here, we extend the results obtained in Ref. [34] and show that the fidelity can have a faster than exponential behavior. The fidelity corresponds to the Fourier transform of the energy distribution of the initial state. The latter is also known as local density of states (LDOS) or strength function [35]. In the case of realistic systems with two-body interactions, the maximum LDOS is Gaussian. In this scenario, the fidelity decay is therefore also Gaussian and this behavior can persist until saturation. The slower exponential decay observed in previous studies occurs when the energy distribution of the initial state is restricted to a Breit-Wigner (Lorentzian) form or to a Gaussian shape that is not well filled. There are, however, situations where even the Gaussian decay can be surpassed. One, addressed here, happens when the system is described by full random matrices. This is not a very realistic approach, since full random matrices imply simultaneous interactions of many particles, but it serves to identify the ultimate lower bound for the fidelity decay in many-body quantum systems where the initial state has a single-peaked energy distribution.
We also discuss results for the evolution of the Shannon entropy. The fidelity gives the probability of finding the initial state in time, whereas the Shannon entropy captures the participation of other states. Analytical expressions were obtained showing that the Shannon entropy increases linearly in time in the limit of strong perturbation [29,32,33]. Even though this behavior is reproduced for the initial states considered here, the results do not match those previous analytical expressions. We speculate on the causes for the discrepancy and how it may be solved.
Another essential aspect of nonequilibrium dynamics, especially in connection with experiments, is the evolution of few-body observables. A complete description is a complex enterprise, since the evolution depends not only on the initial state and final Hamiltonian, but also on the individual properties of the various observables. However, when the initial state is also an eigenstate of the observables, their dynamics depends on the results for the fidelity and a simple general picture becomes available. We find distinct observables evolving according to very different Hamiltonians, but showing a very similar behavior.
The core sections of this paper are Secs. III, IV, and VI. Section III contains the main results about the relationship between fidelity decay and energy distribution of the initial state. Section IV extends this discussion to initial states that are accessible to experiments with cold atoms in optical lattices. Section VI analyzes the short-time dynamics of few-body observables. In the other sections, we cover the description of the model, their density of states, and the meaning of quench dynamics [Sec. II], as well as the results for the evolution of the Shannon entropy [Sec. V]. Concluding remarks are presented in Sec. VII.

II. MODEL, DENSITY OF STATES, AND QUENCH
A way to describe many-body quantum systems is to treat them statistically using full random matrices. This was Wigner's approach to describe heavy nuclei [36] and it was soon employed in the description of other complex systems, such as atoms, molecules, and quantum dots [37][38][39][40][41]. However, full random matrices do not capture the details of realistic quantum systems with few-body interactions, as the spin-1/2 systems considered here. Below we give a general overview of the differences between the two.

A. Full Random Matrices
Full random matrices are matrices filled with random numbers. Their only constraint is to satisfy the symmetries of the system they are trying to describe. The level spacing distribution of these matrices show level repulsion, which is a main feature of what is called quantum chaos [37][38][39][40][41]. The latter corresponds to properties of eigenvalues and eigenstates found in the quantum level that indicate whether the system in the classical level is chaotic or not. The term has in fact been extended to refer to those properties even in quantum systems without a classical limit. Level repulsion leads to the Wigner-Dyson shape of the distribution Π W D (s) of the spacings s between neighboring energy levels. The exact shape depends on the symmetries of the system. Ensembles of real and symmetric random matrices, the so-called Gaussian Orthogonal Ensembles (GOE's), imply time reversal invariance and lead to Π W D (s) = (πs/2) exp(−πs 2 /4).
The density of states P(E) (not to be confused with the local density of states P (E) analyzed in Sec. III) of full random matrices has a semicircular shape [40,42,43]. This is shown in Fig. 1 (a) for a full random matrix from a GOE.
FIG. 1: (Color online) Density of states for a full random matrix from a GOE (a) and for H ∆=0.5,λ=1 in the S z = 0 subspace of L = 16 (b). The GOE is normalized so that the length of the spectrum is 8J and its width is similar to that of H ∆=0.5,λ=1 ; D = 12 870.
The problem with full random matrices is that they imply the simultaneous interactions between many particles, while real systems involve few-body interactions, usually just two-body interactions, as the systems we study below. Early attempts to improve this picture led to the introduction of band random matrices [44] and two-body random ensembles [45][46][47][48]. The systems that we consider do not involve any randomness.

B. Spin-1/2 Model
The realistic spin-1/2 models that we investigate here describe real magnetic compounds [49][50][51], crystals of fluorapatite [52][53][54], and have also been simulated with optical lattices [11,15,17]. We focus on one-dimensional lattice systems with open boundaries and an even number L of sites. The Hamiltonian contains nearest-neighbor (NN) and possibly also next-nearest-neighbor (NNN) couplings, (1) It can be mapped onto systems of spinless fermions [55] or hardcore bosons [56]. In the equation above, = 1 and S x,y,z j = σ x,y,z /2 are the spin operators at site j; σ x,y,z j being the Pauli spin matrices. The coupling strength J, the anisotropy parameter ∆, and the ratio λ between NNN and NN exchanges are chosen positive, thus favoring antiferromagnetic order.
is the Ising interaction between NN (NNN) spins.
The Hamiltonian conserves total spin in the z direction, [ H, Other symmetries include parity, invariance under a global π rotation around the x axis when S z = 0, and conservation of total spin S 2 = ( L j=1 S j ) 2 when ∆ = 1. We work in the S z = 0 subspace, where the dimension of the Hamiltonian matrix is D = L L/2 . The noninteracting XX model (∆ = λ = 0) is trivially solved. The interacting XXZ case (∆ = 0, λ = 0) is solved with the Bethe ansatz [57]. The system undergoes a crossover to the chaotic regime as λ increases [33,41,58,59], the level spacing distribution gradually changing from a Poisson distribution, Π P (s) = exp(−s), in the case of the XXZ model [60] to the GOE Wigner-Dyson form.
We analyze the dynamics of the system for the following choices of parameters for the final Hamiltonian, Note that the value of λ leading to chaos depends on the system size. The larger the system, the smaller the parameter needs to be [61]. The density of states of the above Hamiltonians, independent of the regime of the system, is Gaussian, as shown in Fig. 1 (b) for H ∆=0.5,λ=1 (see other illustrations in [58]). This is typical of systems with two-body interactions [45,47,[62][63][64][65]. The majority of the states are close to the middle of the spectrum, where strong mixing occurs. The eigenstates (written in the mean-field basis) reach their highest level of delocalization in the center of the spectrum and are more localized close to the edges.
The width ω and the average energy E obtained from a Gaussian fit for the density of states of the Hamiltonians above are shown in Table I. The distributions get broader as the value of the anisotropy parameter and the strength of NNN coupling increase. They also shift their center away from zero and become more assymetric. The scenario we consider here is that of a quench. The initial state |Ψ(0) = |ini is as an eigenstate of an initial (unperturbed) Hamiltonian H I . The dynamics starts with the sudden change of some parameter(s) of this Hamiltonian in a time interval much shorter than any characteristic time scale of the model. This results in the final (perturbed) Hamiltonian H F with eigenvalues E α and eigenstates |ψ α = |ini . The unitary time evolution is given by The coefficients C ini α = ψ α |ini are the overlaps of the initial state with the eigenstates of H F . The evolution is computed numerically with full exact diagonalization for L ≤ 16 (D = 12 870) and with EX-POKIT [66,67] for larger system sizes. We examine up to L = 24 (D = 2 704 156). EXPOKIT is a software package based on Krylov subspace projection methods. Instead of diagonalizing the complete system Hamiltonian, the package computes directly the action of the matrix exponential e −i HFt on a vector of interest.

III. LDOS AND FIDELITY
The fidelity is one of our main quantities of interest. It corresponds to the probability of finding the system still in the initial state after time t. It is given by the overlap, F (t) is therefore equivalent to the Fourier transform in energy of the components |C ini α | 2 . The distribution P ini (E α ) of the components |C ini α | 2 in the eigenvalues E α is the so-called LDOS. Experimental measures of the density of states and LDOS is a subject of intense investigation, particularly in nuclear physics [62,[68][69][70]. From the connection between LDOS and the probability distribution of work, we can infer also the possibility of measuring the LDOS experimentally with Ramsey interferometric techniques [71,72].

A. Semicircular LDOS
For an initial state projected onto a full random matrix, P ini (E α ) agrees with the density of states and has again the semicircular shape. The envelope of the distribution is the function where 2E is the length of the spectrum and is the uncertainty in energy. An illustration is provided in Fig. 2 (a). The fidelity for this distribution is given by where J 1 is the Bessel function of the first kind. The behavior of the fidelity is in excellent agreement with the numerical results, as shown in Fig. 2 The maximum fidelity decay, when the LDOS has a single energy peak, is therefore given by Eq. (6). Apart from very short times, where F SC (t) ∼ 1 − σ 2 ini t 2 , even Eq. (6) is slower than the bound, F (t) ≥ cos 2 (σ ini t), derived in [1][2][3][4][5]. This latter result can be approached when the energy distribution of |ini involves more peaks well separated in energy. This is beyond the scope of this work, where the initial states are eigenstates of the unperturbed part of the final Hamiltonian.
As mention before, in realistic systems with few-body interactions, the density of states is Gaussian instead of semicircular. This has consequences to the LDOS, which cannot therefore excel the Gaussian shape.

B. Breit-Wigner LDOS
The energy of the initial state projected on the final Hamiltonian is In systems with two-(few-)body interactions and E ini close to the middle of the spectrum, as the strength of the instantaneous perturbation applied on H I increases from zero, LDOS broadens from a delta function to a Breit-Wigner form delineated by [32,33,35,62,73], where Γ ini is the width of the distribution. An example is given in Fig. 2 (c), where the quench considered is from the initial Hamiltonian H ∆=0.5,λ=0 to the final weakly chaotic Hamiltonian H ∆=0.5,λ=0.4 . The energy E ini of the initial state used is away from the edges of the spectrum. The above distribution leads to the exponential fidelity decay [24][25][26][27][28][29][30][31], as illustrated in Fig. 2 (d). For very short times, F BW (t) ∼ 1 − σ 2 ini t 2 , as expected from perturbation theory, but it soon switches to the exponential behavior.

C. Gaussian LDOS
As the perturbation to H I increases even further, the LDOS of initial states with E ini away from the edges of the spectrum eventually approaches a Gaussian form [28,29,32,33,35,62,73,74] of width Above |n corresponds to the eigenstates of H I and the basis in which the final Hamiltonian is written. The width σ ini depends only on the sum of the square of the off-diagonal elements of H F and can therefore be obtained before the diagonalization of this Hamiltonian.
In systems with two-body interactions, the Gaussian envelope of the LDOS, gives the maximum possible spreading of the initial state in the eigenvalues of the final Hamiltonian and is known as the energy shell [27-29, 32, 33, 35, 62, 73, 75]. The ergodic filling of the energy shell is used as a definition of chaotic states. When this happens, the components |C ini α | 2 of the state become random numbers following the Gaussian distribution. In this sense, a chaotic state may emerge even when one of the Hamiltonians involved is integrable, provided E ini is away from the borders of the spectrum.
The Gaussian distribution leads to a Gaussian fidelity decay controlled by σ ini , This Gaussian behavior was derived also via the central limit theorem [76]. Illustrations for P ini G (E) and F G (t) are given in Figs. 2 (e) and (f), respectively. The initial state chosen is an eigenstate of H ∆=0.5,λ=0 and it evolves according to the final Hamiltonian H ∆=0.5,λ=1 . Its energy is away from the edges of the spectrum. We stress that the fidelity decay, as seen in Fig. 2 (f), can be Gaussian until saturation (the saturation point is indicated with the dashed horizontal lines). This is in contrast with previous works, where despite the Gaussian LDOS, the expectation was for a Gaussian decay only at short-times, switching to exponential at longer times [26-29, 32, 33].
As E ini approaches the border of the spectrum, the initial state becomes more localized, the energy shell less filled, and the LDOS acquires a skewed Gaussian shape [75,77]. This is a consequence of the low density of states at the edges of the spectrum. Examples of this dependence on energy are provided in Figs. 3 and 4 in Sec. IV.

D. Relaxation Time
The equation for the fidelity involves diagonal and off-diagonal terms, The latter averages out in a system without too many degeneracies. After relaxation, the fidelity then saturates to its infinite time average F = |C ini α | 4 = IPR −1 ini , where the inverse participation ratio, IPR ini , measures the level of delocalization of the initial state in the energy eigenbasis. A large value of IPR ini indicates a delocalized state.
The saturation point is minimum when |ini is quenched to a full random matrix (or when the initial state is extracted from a full random matrix). The eigenstates of full random matrices are random vectors, thus |C ini α | 2 fluctuates around 1/D and for GOEs, IPR ini ∼ D/3 [62,63]. The relaxation time t R for an initial state evolved according to full random matrices can thus be obtained from In the case where both H I and H F are two-body-interaction Hamiltonians, the level of delocalization of the initial state depends on its energy. When E ini is close to the middle of the spectrum, IPR ini is large, although usually smaller than D/3, and it gets smaller as E ini approaches the borders. The minimum relaxation time for systems with two-body interactions and a single-peaked LDOS is therefore, The width and filling of the energy shell determine the lifetime of |ini .

IV. EXPERIMENTALLY ACCESSIBLE INITIAL STATES
In this section we extend the studies about fidelity decay to initial states that can be prepared experimentally with cold atoms in optical lattices. They are states where each lattice site has a spin either pointing up or pointing down in the z direction [58,[78][79][80]: The proposals for the preparation of domain walls in optical lattices require the application of a magnetic field gradient [81]. The Néel state [11,16,82,83] has been prepared with high precision using a bichromatic optical superlattice [16] and the evolution of quasi-local densities, currents, and coherences was then experimentally investigated. We recall that in the context of quench dynamics, the initial state is an eigenstate of H I . The eigenstates of the initial Hamiltonian also define the basis in which H F is written. For initial states where each excitation is confined to a single site, as above, H I corresponds to the Ising interaction of Hamiltonian (1). We refer to these states as site-basis vectors (they are also often called computational basis or natural basis). In this basis, the diagonal elements of the final Hamiltonian matrix depend on ∆ and λ, while the off-diagonal elements depend only on λ, since they follow from the flip-flop terms.
We study how the initial states above evolve according to the final Hamiltonians of Table I. This corresponds to a nonperturbative quench, where the off-diagonal elements of H F are much larger than the average level spacing. The shape of the energy distributions of these initial states is close to Gaussian, although the energy shell is not always well filled. Better fillings are associated with E ini closer to the middle of the spectrum. We discuss these distributions in detail in the next subsection before presenting the results for the fidelity.
For site-basis vectors, it is straightforward to calculate analytically n| H F |ini and, from it, the center E ini and the width σ ini of the energy shell. One sees that Eq. (10) reduces to where the connectivity M 1 (M 2 ) corresponds to the number of states directly coupled to |ini via the NN (NNN) flip-flop term. Notice that the total connectivity M of any site-basis vector is low, The values of E ini and σ ini for the three states above are given in Table II. Eini σini |DW , |PS , and |NS are chosen to magnify the effects of the anisotropy and of the NNN couplings. As ∆ increases, E ini is pushed to the edges of the spectrum, where the states get more localized. This dependence between E ini and ∆ is seen in Table II and in the panels for the energy distribution for |PS in Fig. 3 and for |NS in Fig. 4 (distributions for |DW are shown in Ref. [34,58]). For |DW and |PS , E ini is further pushed to the borders as λ increases, whereas for |NS , λ counterbalances the NN contributions and actually brings the state closer to the middle of the spectrum. Among the five cases in Fig. 3, E |PS is farthest from the center of the spectrum for the chaotic isotropic Hamiltonian [ Fig. 3 (d)], while for the Néel state this happens for the integrable isotropic H F [ Fig. 4 (a)]. Depending on λ, the energy of the initial state may also depend on L. In the integrable domain, a linear dependence on L occurs for |DW and |NS , while in the chaotic regime, it happens for |DW and |PS (see Table II).   Tables II, III); bin size = 0.05 J; L = 16.
The width of the energy shell does not depend on the anisotropy parameter, but it may be affected by the presence of NNN couplings. The shell broadens with λ for |DW and |PS [ Fig. 3], but not for |NS [ Fig. 4], since for this state M 2 = 0. The width is extensive on the system size for |NS (M 1 = L − 1) and |PS (M 1 = L/2, M 2 = L − 2), but not for |DW (M 1 = 1, M 2 = 2). The energy shell for the domain wall has the smallest σ ini among the three initial states. |NS has the largest M 1 and thus the largest σ ini when λ = 0, but it is surpassed by |PS when λ > 1/ √ 2. Despite being, in general, close to Gaussian, the LDOS for the initial states considered differ with respect to the filling of the energy shell. The latter depends on the interplay between ∆ and λ, and of course also on L. The values of the least square, used to quantify the deviation of the LDOS from the energy shell, are given in Table III. Small values indicate good filling of the shell. The exact values depend on the chosen bin size, which is somewhat arbitrary. It needs to be sufficiently small so that regions inside the shell where |C ini α | 2 is very small can be detected. Our choice was made to guarantee that for the same σ ini , the least square was smaller if IPR ini was larger. In Table III one finds also the values of IPR ini , which give information about how much spread the initial states are in the energy eigenbasis. Since ∆ pushes E ini to the edges of the spectrum, for a fixed σ ini , better filling occurs for smaller anisotropy. This is confirmed with the values of least square in Table III and  The dependence of the filing on λ is more subtle. For |NS , since the NNN couplings simply push E ini to the center of the spectrum, the behavior is monotonic: the filling improves and IPR |NS increases with λ be ∆ equal to 0.5 or 1. For |DW and |PS this behavior holds only for ∆ = 0.5. The chaotic anisotropic H ∆=0.5,λ=1 leads to the lowest least square value for the three |ini [cf. Fig. 3 (e), Fig. 4 (e), and Table III], |NS being the most delocalized one. When ∆ = 1 the improvement with λ occurs only as the parameter goes from 0 to 0.4 and the connectivity increases, while from λ = 0.4 to 1 the least square value increases. This happens because, despite the broadening of the LDOS, λ pushes E |DW and E |PS to edge of the spectrum, where the states are more localized. This unfavorable combination causes the worst filling of the shell for |PS to occur for the chaotic isotropic Hamiltonian. In this case the energy distribution is skewed and spiky [ Fig. 3 (d)] and IPR |PS has the lowest value [ Table III].
As a final remark, we note the unexpected relation between the width of the density of states [ Table I] and that of the energy shell. For L = 16, ω > σ ini for |DW under the six Hamiltonians, but this is not the case when the initial state is |PS or |NS . For the first, ω < σ ini for H ∆=0.5,λ=1.0 and for the latter this happens for all Hamiltonians, except the strongly chaotic ones.

B. Fidelity Decay
The fidelity decay reflects the results of the energy distribution of |ini . When the energy shell is well filled, the decay is Gaussian and this behavior may persist until saturation, whereas poor filling causes a mixture of Gaussian and exponential behavior, as illustrated in Fig. 5.  Table III The fidelity decays slowly for |DW , due to its low connectivity and narrow energy distribution. At short time the behavior is Gaussian and equal for systems with the same λ (same σ |DW ), but soon the curves for isotropic and anisotropic systems diverge, the first being slower than the latter, as expected from the filling of the shell. It is close to this point of separation that the exponential behavior takes over, but for the domain wall it does not remain until saturation. This state has a complicated dynamics at longer time, with the emergence of some plateaus indicating possible regions of pre-relaxation. Notice also that the infinite time average of the fidelity is very similar for the isotropic Hamiltonians, but the time for it to be reached depends on the strength of the NNN couplings, being shorter for larger λ. Among the Hamiltonians, H ∆=1,λ=1 leads to the fastest relaxation to equilibrium, since it combines largest σ |DW and largest saturation value.
The fidelity decay for |PS shares common features with |DW : at short time it is Gaussian and equal for Hamiltonians with the same λ, later it switches to an exponential behavior, and fastest relaxation to equilibrium happens for the chaotic isotropic system. However, contrary to |DW , the exponential decay of |PS persists until close to equilibration and IPR −1 |PS differs among the isotropic Hamiltonians. Furthermore, according to Table II, the fidelity decay rate increases with L for |PS , whereas σ |DW does not depend on the system size.
For the Néel state, where σ |NS is dissociated from ∆ and λ, the curves for F (t) fall on top of each other for the five Hamiltonians considered in the figure. The decay is Gaussian until saturation. As mentioned before, this is in contrast with previous studies where the exponential behavior superseded (or was expected to supersede) the Gaussian decay before saturation [25-29, 32, 33, 84-89]. We note that the slight acceleration of H ∆=0.5,λ=0 and H ∆=1,λ=0.4 must be caused by the presence of spikes far in the energy distribution. When these peaks are large and far in energy, they can add fast cosine decays to the Gaussian behavior.
The persistence of the Gaussian decay for σ ini t > 0 is not particular to the Néel state. We verified it for several sitebasis vectors. In fact, the majority of the site-basis vectors are much more delocalized than |NS and have comparable σ ini .
The Néel state emphasizes the role of the interplay between initial state and Hamiltonian. It is remarkable to find integrable and chaotic, isotropic and anisotropic systems, all showing the same dynamics. The difference appears only at the saturation point. The infinite-time average value decreases monotonically from ∆ = 1 to ∆ = 0.5 and from λ = 0 to λ = 1 (see IPR |NS in Table III). As a result H ∆=1,λ=0 is the first to reach equilibrium and H ∆=0.5,λ=1 is the last one. After relaxing, the fidelity fluctuates around IPR −1 |NS . The size of these fluctuations decreases exponentially with L [58].
An advantage of using site-basis vectors as initial states is the access that they give to exact analytical expressions for E ini and σ ini . In general, one needs exact full diagonalization to find these values, which limits the system sizes that can be studied. For the dynamics, on the other hand, there are alternative methods, such as Krylov subspace techniques or density matrix renormalization group, that can deal with larger L. Having access to σ ini without the need to resort to exact diagonalization allows us to compare the analytical expression in Eq. (12) with numerical results for F (t) for L > 16. In the bottom right panel of Fig. 5, we use EXPOKIT [66,67] and confirm the Gaussian fidelity decay for the Néel state up to saturation also for L = 24.

V. SHANNON ENTROPY
The Shannon (information) entropy is a delocalization measure that, just like IPR, depends on the basis. Written in the eigenstates of the final Hamiltonian, the Shannon entropy for a certain |ini corresponds to the diagonal entropy [90,91], which, after the quench, is a static quantity. In contrast, if the Shannon entropy is written in the eigenstates of H I , it will evolve in time. In this case, it quantifies the gradual spreading of the initial state and increased participation of the basis vectors in time.
Here, as in Sec. IV, we consider site-basis vectors |n , the initial state being one of them. The evolution of the Shannon entropy is given by where W n (t) = n e −i HF t ini 2 , and W ini (t) = F (t).
Obtaining an analytical expression for Sh(t) is non-trivial due to the dependence on the overlap between the evolved initial state and the other basis vectors. But we can estimate the short-time dynamics by expanding Eq. (16) and using the expression for fidelity F (t) given in Eq. (12). We obtain From the equation above, it is clear that ∆ does not affect the evolution of Sh when σ ini t < 1. States evolving under isotropic or anisotropic Hamiltonians show very similar behavior, as seen in Fig. 6 for |ini = |DW , |PS , and |NS . The role of the anisotropy becomes noticeable at longer times, particularly for the domain wall, where the evolution is slow and gives enough time for a clear separation of the curves before saturation.
In terms of regime, the entropy grows faster in the chaotic domain for |DW and |PS . In contrast, the initial evolution for the Néel state does not depend on λ, since M 2 = 0. The expression for the entropy simplifies to showing dependence only on the system size. The curves for |NS in Fig. 6 coincide. For short times, Fig. 6 indicates good agreement between the numerical results and the approximated Eq. (18). At later times, there is a visible transition from a quadratic to a linear behavior before saturation. This linear increase of the entropy is, however, not well described by the expressions obtained in Ref. [29] for initial states with a Breit-Wigner energy distribution and used also in [32,33] for initial states with a Gaussian LDOS. The derivation of the analytical expression in [29] is based on a cascade model that separates sets of states |n in classes. Each class is successively populated by states directly coupled with those from the previous class, leading to where Γ ini is the width of the Breit-Wigner. In Refs. [29,32,33], the connectivity of the initial state was large, so the last two terms were smaller than the first one. The increase of the entropy was then well captured by In Refs. [32,33], Γ ini was substituted by the width σ ini of the Gaussian LDOS. A semi-analytical expression was also proposed to describe the dynamics at both short and long times [29]. It corresponds to where N pc = exp(Sh) is obtained numerically by performing an average in a large time interval after relaxation. None of the three expressions can describe the results in Fig. 6. The problem is caused by the low connectivity of the initial states, where M ∝ L. The above expressions are valid for initial states with large connectivity, M ∝ D. In this latter case, the cascade model assumption of negligible probabilities of return from one class to previous classes holds, while it is violated for our |ini 's.
In Fig. 7 we compare our results for the Shannon entropy with Eqs. (20), (21), and (22). For the semi-analytical expression (22), we calculated F (t) and N pc numerically. The illustrations are for |PS [ Fig. 7 (a) and (c)] and |NS [ Fig. 7 (b) and (d)] evolving under H ∆=0.5,λ=1 and H ∆=1,λ=1 . In this case, |PS has the largest connectivity among the states studied, M = 3L/2 − 2, and |NS is the state that best fills the energy shell. The disagreement between the numerical results and the three expressions is evident.
One cannot discard, however, the possibility of extending the cascade model to include also low-connectivity-states. For this, it will be necessary to take into account that the connectivities may not be approximately constant for each class, especially for the first classes. For site-basis vectors, for example, we find great discrepancies. The Néel state is directly coupled to L − 1 states, while these states couple with L 2 − 5 states. These two connectivities are very different and much smaller than D. It is worth pointing out also that Eq. (22) agrees with Eq. (19) if N pc = L − 1. This suggests that a reassessment of the semi-analytical expression may need to deal with different values of N pc at different time intervals. A model capable of describing the evolution of the Shannon entropy until saturation even for states with low-connectivity is likely to shed light also on what to expect for the evolution of observables.

VI. FEW-BODY OBSERVABLES
The analysis of the evolution of few-body observables is, of course, more involved than the study of the fidelity decay. It depends on the overlaps between the evolved |ini and other basis vectors, as in the Shannon entropy, and also on the details of the observables A. However, a simple general picture, valid at short times, can be constructed for observables that commute with H I . In this case, the fidelity, and therefore σ ini , plays an important role in A(t).
The evolution of the observables is given by In the particular case where H I is the Ising interaction of Hamiltonian (1) and |n are site-basis vectors, as in Secs. IV and V, the off-diagonal elements of the final Hamiltonian, n| H F |ini , and σ ini are obtained from the flip-flop terms. Thus Eq. (24) does not depend on the anisotropy parameter. In other words, the short-time dynamics of few-body observables that commute with H I does not depend on the interaction. For the Néel state, where the offdiagonal elements exist only for the NN flip-flop term, not even λ is important, and completely different Hamiltonians (integrable and chaotic, isotropic and anisotropic) lead to a very similar initial relaxation of the observables.

A. Local magnetization
The on-site magnetization, S z j , commutes with H I . This is a simple and yet useful observable frequently measured experimentally [92]. It is straightforward to show that for |DW , |PS , and |NS the magnetization in the middle of the chain behaves, at short times, as In the integrable regime, S z,|NS L/2 (t) changes faster than for the other two states, since for |NS , the excitation on site L/2 has two neighboring sites to hop to, while for |DW and |PS it has only one site. When directly couplings between NNN are included and λ > 1/ √ 2, S z,|PS L/2 (t) becomes the fastest to evolve, since the excitation has now three sites to hop to, while |NS and |DW have only two. At the edges of the chain, border effects slow down the dynamics.
The above approximations agree well with our numerical results at short times (not shown). The curves with the same value of λ coincide for |DW and |PS , whereas for |NS , the curves fall on top of each other, independently of ∆ or λ.

B. Spin-spin correlations
The spin-spin correlation between sites i and j is given by In the z direction it commutes with H I . Using Eq. (24) and Table II, we find for neighboring sites in the middle of the chain that The expression for |PS above holds for mod(L, 4) = 0, when the spins on sites L/2, L/2 + 1 are parallel. When mod(L, 4) = 0, and the two middle spins are anti-parallel, the expression changes to C z,|PS L 2 , L 2 +1 (0)[1 − J 2 t 2 (1 + 3λ 2 )/2]. As expected and confirmed numerically (not shown), the magnitude of C z,|DW L 2 , L 2 +1 (t) decays slowly, especially in the integrable domain. Figure 8 compares the longitudinal correlation for |PS and |NS evolving under the same final Hamiltonians. Similarly to what was seen for fidelity, the initial decay of the magnitude of C z L 2 , L 2 +1 (t) for the Néel state is independent of the regime of H F , while for |PS it is faster in the chaotic domain. These distinct behaviors, anticipated from Eqs. (27) and (28), emphasize the significance of the initial state also for the dynamics of observables.
After a long time, C z L 2 , L 2 +1 (t) fluctuates around the equilibrium value C z L 2 , L 2 +1 , the fluctuations decreasing exponentially with system size [58]. The saturation value is closest to zero when E ini is closest to the center of the spectrum. This happens for |PS with the integrable Hamiltonians and for |NS with the strongly chaotic Hamiltonians, as can be seen in the insets of Fig. 8.
The insets of Fig. 8 give the scaling of C z L 2 , L 2 +1 with system size. For |NS , the correlations for the strongly chaotic Hamiltonians approach zero as L increases, while the results indicate that integrable and weakly chaotic Hamiltonians may retain memory in the thermodynamic limit. However, we cannot discard the possibility of an acceleration towards zero for L's larger than the ones considered here. The results for |PS are less conclusive. For this state, the direction of the spins in the middle of the chain depend on L. This causes the saturation value for small L to oscillate significantly from mod(L, 4) = 0 to mod(L, 4) = 0. To try to delineate a pattern, we show only the results for mod(L, 4) = 0. The correlations decrease with system size, but more points are necessary for an extrapolation to the thermodynamic limit.

C. Structure Factor
The structure factor is the Fourier transform of the spin-spin correlations, being therefore a nonlocal observable. In the z direction, commutes with H I . Above, κ = 2πp/L stands for momentum and p = 0, 1, 2 . . . , L is a positive integer. The evolution of the structure factor depends on κ. For example, for the Néel state, Eq. (24) becomes . It is only for κ = π that s z,|NS f (κ, 0) = 0, which allows for a large contribution from the width of the energy shell. This dependence on κ is reflected in Fig. 9. Since the magnitude of O(t 2 ) becomes very small as p → 1, higher order terms, with effects from ∆ and λ, become significant already at short times, which explains the early divergence of the curves for different Hamiltonians. In contrast, for κ = π the curves coincide up to the vicinity of the saturation point.
Contrary to the local observables C z i,j (t) and S z j (t), the evolution of the structure factor depends also on L. The decay is faster as L increases, although how fast it is depends again on κ. When κ = π, Eq. (30) simplifies to which shows the linear dependence of the term O(t 2 ) on L. For other κ's the dependence is much smaller. A similar analysis can be extended to |DW and |PS . In general, for these states, the evolution of s z,|ini f (κ, t) at short times depends on λ, as seen previously for other observables. The dependence on κ is again present, although it is not the same found for the Néel state. For instance, for |PS , the effects of σ |PS occur only when κ = π/2, where

D. Local spin current
The spin current is an observable of great interest for studies about quantum transport [93][94][95][96][97][98], which motivates having a closer look at it. Even though its evolution cannot be cast in the form of Eq. (24), the short-time dynamics shows again a simple dependence on λ that is comparable to what was found for the previous observables. In particular, the behavior for |NS is again very similar for all Hamiltonians considered.
The local spin current, I s,j , is associated with the conservation of total spin in the z direction, S z , and obeys the continuity equation [78], In the bulk, the local spin current agrees with the result from a periodic chain [93], −i[H, S z j ] = div(I s,i ) = (I s,j − I s,j−1 ), which leads to This observable does not commute with H I and I s,j (0) = 0. Also in contrast with the previous observables, n|I s,j |ini is imaginary, so the dominant term in the expansion of Eq. (23) is If the pair of spins on sites (j, j + 1) are parallel, as happens for |PS when mod(L, 4) = 0, then n|I s,j |ini = 0. This explains why, in Fig. 10, where L = 16, we selected I s, 9 . It is straightforward to show that for pairs of anti-parallel spins on sites (j, j + 1), The above expressions capture well the initial dynamics of the local spin current shown in Fig. 10. They make evident the lack of any influence of ∆ at short times, the dependence on λ for |PS , and its absence for |NS . For the Néel state, the proximity of the curves for very different Hamiltonians and even beyond the range of validity of Eq. (34) is again remarkable.

E. Entangled State
The simple picture developed for the short-time dynamics of C z i,j becomes less trivial when dealing with C x i,j . One reason is the disappearance of the role of the fidelity (and σ ini ), since C x i,j (0) = 0 when the initial state is a site-basis vector. The other issue are the contributions, already at O(t 2 ), from terms such as n| H 2 F |ini , which bring into play the effects of the anisotropy. Nevertheless, one can construct particular initial states, where behaviors close to that of Eq. (24) can be recovered also for C x i,j , and similar observables. Here, this is illustrated for an initial state corresponding to an entangled state, |ES which contains an EPR pair on sites L/2 and L/2 + 1, On each side of the EPR pair, there is a domain wall. To simplify the analysis, we fix system sizes such that mod(L, 4) = 0.   The energy and width of the energy distribution of |ES are Their dependence on ∆, λ, and L is similar to that for |DW , apart from an additional term ∝ λ in σ 2 |ES . The energy shell for |ES is, however, much better filled (compare Table IV and Table III). For the entangled state with H ∆=0.5,λ=1 , the least square is even better than that for the Néel state.
The fidelity decay of |ES is shown in Fig. 11 (a). The decay is slower than that of |PS and |NS (cf. Fig. 5) due to the domain wall structure that leads to smaller values of σ |ES . Yet, in the chaotic domain, the behavior is Gaussian until very close to saturation, as anticipated from the low values of least square.
Using the approximation · ↑↓ · | e −i HF t | ES ≈ · ↓↑ · | e −i HF t | ES , which holds at short times, we can expand C z L 2 , L 2 +1 (t) to O(t 2 ) as in Eq. (24). This gives For |ES , since |· ↑↓ · is directly coupled with |· ↓↑ · , C x L 2 , L 2 +1 (0) = 0 and we can find an expansion similar to that in Eq. (24) also for the correlation in the x direction. It gives At short times, both expressions above show good agreement with the numerical results in Figs. 11 (b) and (c), respectively. Since there is no dependence on ∆, the curves for equal λ are very similar. In fact, they practically coincide until close to saturation.

VII. CONCLUSION
We studied isolated quantum systems with interaction and quenched far from equilibrium. The initial state |ini was an eigenstate of an initial Hamiltonian H I and it evolved according to a final Hamiltonian H F . We analyzed numerically and analytically the behavior in time of fidelity, Shannon entropy in the basis of H I , and few-body observables. The focus was on initial states, system models, and observables that are accessible to current experiments with optical lattices. We showed that the system dynamics depends not only on the initial state or on the final Hamiltonian, but on the interplay between the two. Depending on the initial state, different H F 's may lead to very similar evolutions.
The fidelity is the Fourier transform of the energy distribution of the initial state. This distribution is referred to as LDOS. The envelope of the maximum possible LDOS corresponds to the energy shell. We investigated single-peaked LDOS of width σ ini . The ultimate lower bound for the fidelity decay emerges when H F 's are full random matrices and the LDOS is semicircular: F (t) = [J 1 (2σ ini t)] 2 /(σ 2 ini t 2 ), where J 1 is the Bessel function of first kind. If instead of many-body interactions, as implied by full random matrices, only two-body interactions are considered, then the maximum LDOS and, consequently, the fastest fidelity decay are Gaussian: F (t) = exp(−σ 2 ini t 2 ). For initial states that fill the energy shell, the Gaussian behavior persists for long times and can last up to saturation. In this latter case, the relaxation time is t R = ln(IPR ini )/σ ini , where IPR ini measures the level of delocalization of the initial state in the energy eigenbasis. This finding contrasts previous studies, where a switch to an exponential behavior was expected at long times [26-29, 32, 33].
For initial states corresponding to site-basis vectors, we derived an expression that captures the evolution of the Shannon entropy at short times. At longer times, analytical and semi-analytical expressions that had been successfully employed in the past [29,32,33] did not match our results. To reassess these expressions, one will need to take into account the connectivity of the initial states considered. Site-basis vectors are directly coupled with few states, the connectivity is of the order of the system size, whereas previously studied initial states had very high connectivity.
The analysis of few-body observables A was performed for local magnetization, spin-spin correlations, structure factor, and local spin current. Initial states corresponding to site-basis vectors are eigenstates of the first three observables in the z direction, that is [ H I , A] = 0. For short times and A(0) = 0, σ ini plays a central role in the evolution of such observables. The local spin current is not part of this general picture, but its behavior is comparable.
We showed that the Néel state, which is key in magnetism and has also been prepared in optical lattices [16], has a very interesting behavior. At short times, the evolution of fidelity, Shannon entropy, and observables where [ H I , A] = 0 is detached from the values of the anisotropy parameter and the regimes of the final Hamiltonians. The dynamics is the same, be the final Hamiltonian integrable, chaotic, isotropic, or anisotropic.
Our results for fidelity and observables contribute to the development of a general description for isolated quantum systems far from equilibrium. Studies about quantum control methods and quantum thermodynamics are also likely to benefit from our analysis of the relationship between LDOS and fidelity decay. The LDOS is related with the probability distribution of the work needed to take quantum systems out of equilibrium.