Nonequilibrium quench dynamics in quantum quasicrystals

We study the nonequilibrium dynamics of a quasiperiodic quantum Ising chain after a sudden change in the strength of the transverse field at zero temperature. In particular we consider the dynamics of the entanglement entropy and the relaxation of the magnetization. The entanglement entropy increases with time as a power-law, and the magnetization is found to exhibit stretched-exponential relaxation. These behaviors are explained in terms of anomalously diffusing quasiparticles, which are studied in a wave packet approach. The nonequilibrium magnetization is shown to have a dynamical phase transition.


Introduction
Recent experimental progress in ultracold atomic gases in optical lattices [1,2,3,6,5,4,7,8,9,10] has opened up fascinating new perspectives on research in the field of isolated quantum systems, both in equilibrium and out of equilibrium. In experiments the form of atomic interactions can be suddenly changed by tuning an applied magnetic field near a Feshbach resonance, which is known as a global quantum quench. On the theoretical side, one is interested in the time-evolution of different observables, such as the order parameter or some correlation function, after a quench. Fundamental questions concerning quantum quenches include (i) the functional form of the relaxation process in early times, and (ii) the properties of the stationary state of the system after a sufficiently long time.
Many results for quantum quenches have been obtained for homogeneous systems [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]; for example, the relaxation of correlation functions in space and in time is generally in an exponential form, which defines a quench-dependent correlation length and a relaxation (or decoherence) time. Many basic features of the relaxation process can be successfully explained by a quasiparticle picture [33,14,34]: after a global quench quasiparticles are created homogeneously in the sample and move ballistically with momentum dependent velocities. The behavior of observables in the stationary state is generally different in integrable and in non-integrable systems. For non-integrable models, thermalization is expected [12,13,14,15,16,17,18,19,20,21,22] and the distribution of an observable is given by a thermal Gibbs ensemble; however, in some specific examples this issue has turned out to be more complex [23,24,25,31]. By contrast, it was conjectured that stationary state averages for integrable models are described by a generalized Gibbs ensemble [12], in which each integral of motion is separately associated with an effective temperature.
Concerning quantum quenches in inhomogeneous systems, there have been only a few studies in specific cases; for example, entanglement entropy dynamics in random quantum chains [35,36,37] and in models of many-body localization [38,39]. In some of these cases the eigenstates are localized, which prevents the system from reaching a thermal stationary state.
A special type of inhomogeneity, interpolating between homogeneous and disordered systems, is a quasicrystal [40,41] or an aperiodic tiling [42]. Quasicrystals are known to have anomalous transport properties [43,44], which is due to the fact that in these systems the long-time motion of electrons is not ballistic, but an anomalous diffusion described by a power law. One may expect that the quasiparticles created during the quench have a similar dynamical behavior, which in turn affects the relaxation properties of quasicrystals.
Quasicrystals of ultracold atomic gases have been experimentally realized in optical lattices by superimposing two periodic optical waves with different incommensurate wavelengths. An optical lattice produced in this way realizes a Harper's quasiperiodic potential [47,48], for which the eigenstates are known to be either extended or localized depending on the strength of the potential. Different phases of Bose-Hubbard model with such a potential have been experimentally investigated [45,46]. There have also been theoretical studies concerning the relaxation process in the Harper potential [49,50].
In this paper we consider the nonequilibrium quench dynamics of the quantum Ising chain in one-dimensional quasicrystals. The quantum Ising chain in its homogeneous version is perhaps the most studied model for nonequilibrium relaxation [51,52,53,54,55,56,57,58,60,59,34,61,62,63,64,65,66]. Our study extends previous investigations in several respects and seeks to obtain new insights into quench dynamics in inhomogeneous systems. We focus on the Fibonacci lattice, for which many equilibrium properties of the quantum Ising model are known [67,68,69,70,71,72,73]; to our knowledge this is the first study of quantum quenches in such a lattice. Using free-fermionic techniques [75], we numerically calculate the time-dependence of the entanglement entropy as well as the relaxation of the local magnetization for large lattices. The numerical results are interpreted by a modified quasiparticle picture, in which the quasiparticles are represented by wave packets; we also obtain diffusive properties of the wave packets.
The structure of the paper is as follows. The quasiperiodic quantum Ising model and its equilibrium properties are described in section 2. The global quench process and some known results for homogeneous and random chains are presented in section 3. Our numerical results for the quasiperiodic chain are presented and interpreted in section 4. This paper is concluded with a discussion; some details of the free-fermionic calculation of the local magnetization are presented in the appendix.

The Model and its equilibrium properties
We consider the quantum (or transverse) Ising model defined by the Hamiltonian: where σ x i and σ z i are Pauli matrices at site i. The interactions, J i , are generally site dependent, which are parameterized as: where r > 0 is the amplitude of the inhomogeneity, and the integers f i are taken from a quasiperiodic sequence. Quasiperiodic lattices can be generated in different ways, such as by the cut-andproject method. Here we use the following algebraic definition for a one-dimensional quasiperiodic sequence: where [x] denotes the integer part of x, and ω > 1 is an irrational number. The Fibonacci sequence generated by the substitution rule: 0 → 01 and 1 → 0 starting with 0 corresponds to the formula in (3) with the golden mean ω = ( is the fraction of units 1 in the infinite sequence. Note that r = 1 represents the homogeneous lattice. The essential technique in the solution of H is the mapping to spinless free fermions [75,76]. First we express the spin operators σ x,y,z i in terms of fermion creation (annihilation) operators c † i (c i ) by using the Jordan-Wigner transformation [74]: Here and throughout the paper we denote the imaginary unit √ −1 by ı to avoid confusion with the integer index i. The Ising Hamiltonian in (1) can then be written in a quadratic form in fermion operators.
where N = L i=1 c † i c i is the number of fermions. The Hamiltonian (5) can be diagonalized through a canonical transformation [75], in which a new set of fermion operator η k is introduced by where the Φ k (i) and Ψ k (i) are real, and normalized by We then obtain the diagonal form of H: in terms of the new fermion creation (annihilation) operators η † k (η k ). The energies of free fermionic modes, ǫ k , and the components, Φ k (i) and Ψ k (i), can be obtained from the solutions of the eigenvalue problem: The spectrum of free-fermionic excitations, ǫ k in (8), plays a key role in equilibrium and non-equilibrium properties of the system. In equilibrium and in the thermodynamic limit the model has a quantum critical point at h = h c , the properties of which are controlled by the low-energy excitations. The value of h c is determined by the equation [76] ln h c = ln J, where the overbar denotes an average over all sites. With the parameterization given above, the critical point is given by h c = 1, independently of r. The lowest gap, ∆E = ǫ 1 , is zero for h < h c , and vanishes as ∆E ∼ (h c − h) ν , as h approaches h c . The singularity of the gap, measured by the gap-exponent ν = 1, does not depend on r; the same is true for the singularity of the specific heat: C v ∼ ln |h − h c |. Thus the transition belongs to the Onsager universality class [77], irrespectively of r. This means that the quasiperiodic modulation of the couplings represents an irrelevant perturbation at the critical point of the homogeneous model [78]. For h < h c the system is in the ordered phase, so that the local magnetization at site l is m l > 0. Upon approaching the critical point, the local magnetization goes to zero following a power law: the bulk magnetization m b decays as m b (h) ∼ (h c − h) 1/8 , which defines the critical exponent β b = 1/8, while the surface magnetization m 1 vanishes as m 1 (h) ∼ (h c − h) βs with β s = 1/2. For h > h c the system is in the disordered phase and the local magnetization vanishes in the thermodynamic limit.
While in equilibrium only the low-energy excitations are of importance, the complete energy spectrum contributes to nonequilibrium properties, which are investigated in this paper.

Nonequilibrium properties of homogeneous and random chains
We consider a quench process in which at time t = 0 the strength of the transverse field is changed suddenly from h 0 to another value, say h. The initial Hamiltonian with h 0 for t < 0 is denoted by H 0 , and its ground state is Ψ (0) 0 . For t > 0 the new Hamiltonian H with h governs the coherent time-evolution of the system; for example an observable, represented by the operatorÂ, has the time-evolution in the Heisenberg picture as:Â(t) = exp(ıtH)Â exp(−ıtH), and its expectation value for t > 0 is given by . Dynamics of the system out of equilibrium is governed by the complete spectrum of H and not only by the lowest excitations. Therefore, Hamiltonians with different spectral properties will have completely different nonequilibrium properties.
The form of the inhomogeneity in the couplings is generally crucial to the spectrum of a Hamiltonian. For example the spectrum of the homogeneous quantum Ising chain is absolutely continuous, thus all the eigenstates are extended. By contrast, the random chain has a singular point spectrum and the eigenstates are localized. The spectrum of quasiperiodic chains lies between the above mentioned two limiting cases [79,80]; for example, the spectrum of the Fibonacci chain defined in (1) is given by a Cantor set of zero Lebesgue measure, signaling that the spectrum is of a multifractal type, and it is called purely singular continuous [81] in the mathematical denotation. See [79] for precise mathematical definitions of different spectra.
Below we first briefly review nonequilibrium properties of the entanglement entropy and local magnetization after a quench in the homogeneous chain and in random chains.

Entanglement entropy
The entanglement entropy, S ℓ (t), of a block of the first ℓ sites in the chain is defined as S ℓ (t) = Tr ℓ [ρ ℓ (t) ln ρ ℓ (t)] in terms of the reduced density matrix: ρ ℓ (t) = Tr i>ℓ |Ψ 0 (t) Ψ 0 (t)|. Here |Ψ 0 (t) denotes the ground state of the complete system at time t > 0. The details of the calculation of S ℓ (t) in the free-fermion representation can be found in the appendix of [82].
For the homogeneous chain (corresponding to the case with r = 1 in (2)) in the limit L → ∞ and for ℓ ≫ 1 the results can be summarized as follows [33,36]: where v max is a maximum velocity. For a quench to a quantum critical point, the result in (10) is a consequence of conformal invariance [33]; for other cases, this behavior can be explained in the frame of a semiclassical (SC) theory [33,34]: entanglement between the subsystem and its environment arises when two quantum entangled quasiparticles, which are emitted at t = 0 and move ballistically with opposite velocities, arrive in the subsystem and in the environment simultaneously. The prefactors α = α(h 0 , h) and β = β(h 0 , h) have been exactly calculated [83] and these agree with the results obtained from the SC theory [34]. In [84] α(h 0 = 0, h) has been evaluated in a closed formula, which is a continuous function of h but at the critical point h = 1, its second derivative is logarithmically divergent.
In the random chain the excitations are localized and therefore the dynamical entanglement entropy approaches a finite limiting value. When the quench is performed to the random quantum critical point, the average entropy increases ultraslowly as log[log(t)] [36]. This behavior can be explained in terms of the strong disorder renormalization group [85,86,36,39].

Local magnetization
Another quantity we consider is the local magnetization, m l (t), at a position, l, of an open chain. Following Yang [87] this is defined for large L as the off-diagonal matrix-element: where Ψ (0) 1 is the first excited state of the initial Hamiltonian. Calculation of the magnetization in terms of free fermions is outlined in the appendix.
For the homogeneous chain the time-dependence of the local magnetization has been numerically calculated in [34,58]. For the quench performed within the ordered phases, h 0 < 1 and h < 1, the results in the limit L → ∞ and l ≫ 1 are given by: where the relaxation (decoherence) time τ and the correlation length ξ depend on the quench parameters h 0 and h. Exact expressions of these quantities have been derived recently [61,63,64]. In the small h 0 and h limit, accurate results can also been obtained from the SC theory [34,88]. In this framework the quasiparticles in terms of the σ operators are represented by ballistically moving kinks. Each time when a kink passes site l, the σ z l operator changes sign; thus kinks that pass a site an even number of times have no effect on the local magnetization. Summing up contributions of all kinks, we obtain the functional form in (12). If the quench is performed close to the critical point, the kinks have a finite width; this effect can be taken into account in a modified SC theory [34,65], which provides exact results.
For quenches involving the disordered phase with h 0 > 1 and/or h > 1, the results obtained numerically [34,58] or analytically by the form-factor approach [61,63,64] indicate that for bulk spins in large systems the first equation of (12) is modified as: where the prefactor, A(t) changes sign during the relaxation process, say defines a characteristic time-scale, which increases and becomes divergent as h → 1 + . This is a signal of a dynamical phase transition in the system. The order parameter can be defined as: which is positive (O > 0) in the oscillatory phase and O = 0 in the non-oscillatory phase.
In a disordered chain away from the random quantum critical point the bulk magnetization approaches a finite limiting value, which reflects the localized nature of the excitations. After a quench performed to the critical point, the average bulk magnetization has been found to vanish asymptotically in a very slow way [89],

Results for quasiperiodic chains
In this section we present our results for the quasiperiodic quantum Ising chain after a global quench, obtained by numerical calculations based on the free-fermion representation of the model. We concentrate on the Fibonacci chain with the parameter ω defined in (3) being the golden mean. We consider finite chains with a length fixed at a Fibonacci number F n . We have calculated the entanglement entropy and the local magnetization for system sizes up to L = F 17 = 1597. For the numerical calculation we solved hermitien and anti-hermitien eigenvalue problems, and calculated complex determinants using the LAPACK routine. For a given set of parameters (h 0 , h and r) the time-dependence of the entropy or the magnetization of a chain with L = 1597 was obtained in about one day of CPU time on a 2.5 GHz processor.
Below we present results for these two quantities separately.

Entanglement entropy
For a chain of total length F n with periodic boundary conditions, we have calculated the entanglement entropy S ℓ between a block of length ℓ = F n−2 and its environment which has a length of F n−1 . Various values of 0 < r < 1 for the inhomogeneity amplitude were considered. We start our numerical calculations from the fully ordered state with h 0 = 0 to a state with h > 0 both in the ordered and in the disordered phases, as well as at the critical point. The numerical results for S ℓ (t) − S ℓ (0) are shown in figure 1. For all cases considered, S ℓ (t) exhibits two time-regimes: in the late-time regime, the entropy is saturated to an L dependent value, similar to the behavior for the homogeneous chain; in the early-time regime, it increases with time as a power-law form: with some exponent σ < 1. Our numerical results show that the exponent σ depends on the value of the transverse field in the final state, while it does not vary (significantly) with the initial h 0 . The values of σ for r = 0.25, 0.5 and 0.75 are plotted in figure 6; for all cases considered, σ reaches its maximum at the critical point h = 1, and the increase with h in the ordered phase (h < 1) is much faster than the decrease in the disordered phase (h > 1). Furthermore, we have found that the exponent σ decreases with stronger inhomogeneity, that is with smaller value of r.
The power-law time-dependence of the entanglement entropy in (15) is a new feature of the quasiperiodic system: the increase in entropy is slower than in the homogeneous chain, but faster than in a random chain. This behavior can be explained in terms of quasiparticles that are emitted at time t = 0, and subsequently move classically by anomalous diffusion which has a power-law relationship between displacement and time, x ∼ t D , with a diffusion exponent 0 < D < 1. We note that in a homogeneous chain pairs of quasiparticles that contribute to the entanglement entropy move ballistically (i.e. x ∼ t) rather than moving by diffusion, which results in the linear growth of the entanglement entropy with time [33] (figure 2). The dynamics of the quasiparticles in our quasiperiodic lattice will be studied in more detail in section 4.3.

Local magnetization
The local magnetization, m l (t), is calculated for open chains of length L = F n . Generally m l (t) has a monotonic position dependence: m l1 (t) > m l2 (t) for l 1 < l 2 < L/2. We measured the magnetization at site l = F n−1 , which is considered as the bulk magnetization and denoted by m b (t). We have also studied the behavior of the surface magnetization, m 1 (t), for which some exact results are obtained.
We study the asymptotic behavior of the surface magnetization (given in (A.16)) for large t after a quench. If the quench is performed to the ordered phase, h < 1, the lowest excitation energy is ǫ 1 = 0 (i.e. cos(ǫ 1 t) = 1); consequently P 1,2k−1 (t) in (A.7) has a time independent part. This results in a non-oscillating contribution to the surface magnetization: m 1 = lim t→∞ t 0 m 1 (t ′ )dt ′ which is given by: Schematic illustration of the light cones of quasiparticles for a homogeneous quantum Ising chain (a) and for a chain with an aperiodic modulation of the couplings (the thin/thick lines between sites represent weak/strong couplings according to a Fibonacci sequence) (b). The quasiparticle excitations emitted at time t = 0 move ballistically in the homogeneous lattice, while their motion is anomalous diffusive with x ∼ t D (D < 1) in the quasiperiodic lattice. Pairs of quasiparticles moving to the left or right from a given point are entangled; they will contribute to the entanglement entropy between a region A (the region with orange sites) and the rest of the chain, region B, if they arrive simultaneously in A and B. and defines its stationary value. Recall that Φ 1 (1) = m 1 (h, t = 0), i.e. it is equal to the equilibrium surface magnetization [90,91], which is finite for h < 1, and zero in the disordered phase. Similarly, Φ (0) 1 (1) > 0 for h 0 < 1 and zero otherwise. From this it follows that the stationary nonequilibrium surface magnetization is m 1 > 0, if both h < 1 and h 0 < 1. Otherwise the stationary surface magnetization vanishes. If the quench starts from the fully ordered initial state h 0 = 0, then Φ (0) 1 (j) = δ 1,j and m 1 = Φ 2 1 (1); thus we obtain the simple relation: which is generally valid between the stationary value of the nonequilibrium surface magnetization and its equilibrium value. From (17) it follows that the critical exponent β ne s for the nonequilibrium surface magnetization and the critical exponent β s for the equilibrium surface magnetization are related as: β ne s = 2β s . According to (17) and [92], for the Fibonacci chain close to the critical point h → h c = 1, we have We numerically calculated the time-dependence of the bulk magnetization after a quench from the fully ordered initial state, h 0 = 0, to different values of h. For fixed values of the inhomogeneity r = 0.25, 0.5, 0.75, the results for the double logarithm of |m b (t)| are shown in figure 3(a-c) as functions of ln t. In each case one can observe a linear dependence, which implies that the magnetization has asymptotically which corresponds to equation (13) for a homogeneous system, with µ = 1. Before analyzing the decay exponent µ, we first study the behavior of the prefactor A(t).
Like in a homogeneous chain as discussed in section 3.2, there is a dynamical phase transition between a non-oscillating phase for h < h * (r), where the order-parameter O defined in (14) is zero, and an oscillating phase for h > h * (r), where O > 0. In the oscillating phase, the characteristic time-scale defined as the period time, t per (h, r), becomes divergent as h → h * (r) + . An example for this behavior is illustrated in figure 3, panel (d), in which ln |m b (t)| as a function of t is shown in the window 50 < t < 100 for different values of h at r = 0.5; as seen in this figure, the curves for h = 0.86, 1.0 and 1.25 oscillate, whereas the oscillations vanish for h = 0.81 and h = 0.84. We identify the dynamical phase transition point as h * = 0.850 (5). In this quasiperiodic model the dynamical phase transition does not coincide with the equilibrium phase transition, since h * (r) < 1 for r < 1. Estimates of h * (r) versus r are shown in figure 4; the data are well approximated by a power-law h * (r) ∼ r α with

Interpretation by wave packet dynamics
As is known from previous studies on the homogeneous chain, dynamical features of the entanglement entropy and the local magnetization can be well described by the dynamics of quasiparticles. To understand the dynamical properties of the quasiparticles emitted after a quantum quench in the quasiperiodic lattice, we regard the quasiparticles as wave packets and study their dynamics using a method that has been applied to studies of transport properties of quasicrystals [44,94].
We construct a wave packet connecting sites k and l at time t in the form: which is localized at t = 0 since W l,k (0) = δ l,k (cf. equation (7)). For a Hamiltonian with eigenfunctions φ q (l) and eigenvalues ǫ q , a wave packet can be obtained by: W l,k (t) = q cos(ǫ q t)φ q (l)φ q (k), which corresponds to the first term in (19). We note that (19) is just a linear combination of the four time-dependent factors in (A.7), which describe the time dependence of the fermion operators. The width of the wave packet starting from site k after time t is given by: The spreading of a wave packet in a perfect crystal with absolutely continuous energy spectrum is known to be ballistic, i.e. the width increases linearly in time. A heuristic argument is the following [95]: the energy scale ∆ǫ defined by the typical variation of the energy levels is proportional to the inverse of the time that a wave packet needs to spread over the chain. In the case of the absolutely continuous spectrum, we have ∆ǫ ∼ L −1 , which gives d ∼ L ∼ t. In case of a singular continuous spectrum as for our quasiperiodic lattice, there are many energy scales ∆ǫ ∼ L −1/α with a number of exponents α. One then expects that for large t the wave packet in the infinite quasiperiodic lattice shows anomalous diffusion in the form d(k, t) ∼ t D(k) with a diffusion exponent D(k), which may depend on the starting position. Here we determine the value of D(k) numerically. After a global quench, quasiparticles are emitted everywhere in lattices, therefore d(k, t) should be averaged over different initial positions, In our numerical calculations chains of length L = F 17 = 1597 with periodic boundary conditions were considered. First we have confirmed that the wave packet constructed in our method moves ballistically in the homogeneous chain (with r = 1), corresponding to D = 1. In the quasiperiodic chains the motion is indeed anomalous diffusive with D < 1, which is seen in figure 5 where  The variation of D with h at a fixed r is shown in figure 6, compared with the exponent σ for the entanglement entropy and the exponent µ for the local magnetization. Here one can observe that the agreement between these three exponents is very good for h < h * (r), i.e. in the non-oscillating phase, but the exponent for the magnetization deviates in the oscillating phase (h > h * (r)). The discrepancy in the oscillating phase implies that the semiclassical picture breaks down in the oscillating phase, where the quasiparticles cannot be well described by the moving kinks in the magnetization.

Discussion
In this paper we have studied the nonequilibrium dynamics of quasiperiodic quantum Ising chains after a global quench. In a quench process, the complete spectrum of the Hamiltonian is relevant for the the time evolution of various observables. For the quasiperiodic quantum Ising chain the spectrum is in a very special form, which is given by a Cantor set of zero Lebesgue measure, i.e. purely singular continuous. We have calculated numerically two quantities: the dynamical entanglement entropy and the relaxation of the local magnetization. The entanglement entropy is found to increase in time as a power-law (see (15)), whereas the bulk magnetization decays in a stretched exponential way (see (18)). Both behaviors can be explained in a quasiparticle picture, in which the quasiparticles move by anomalous diffusion in the quasiperiodic lattice. The diffusion exponent has been calculated by a wave packet approach, and good agreement has been found with the exponents that we obtained for the entropy and for the magnetization. We note that the anomalous dynamics found in the global quench process is similar to the transport properties of quasicrystals.
Relaxation of the bulk magnetization is found to present a nonequilibrium dynamical phase transition. The non-oscillating phase, in which the magnetization is always positive, and the oscillating phase, in which the sign of the magnetization varies periodically in time, is separated by a dynamical phase transition point, at which the time-scale of oscillations diverges. This singularity point, due to collective dynamical effects, is different from the equilibrium critical point.
A similar nonequilibrium dynamical behavior is expected to hold for other quasiperiodic or aperiodic quantum models as long as the spectrum of the Hamiltonian is also purely singular continuous; there is a large class of such models, for example the Thue-Morse quantum Ising chain. If, however the spectrum of the Hamiltonian of the model is in a different type, such as the Harper potential which has extended or localized states, the nonequilibrium dynamics is expected to be different than the case we consider in this paper. To calculate the local magnetization in (11), we need to first calculate the time dependence of the spin operator σ x l (t) at site l in the Heisenberg picture. We introduce at each site two Majorana fermion operators,ǎ 2l−1 andǎ 2l , defined in terms of the free fermion operators η † k and η k (given in (6)) aš
As a special case, the surface magnetization is expressed as: