Critical quench dynamics of random quantum spin chains: Ultra-slow relaxation from initial order and delayed ordering from initial disorder

By means of free fermionic techniques combined with multiple precision arithmetic we study the time evolution of the average magnetization, $\overline{m}(t)$, of the random transverse-field Ising chain after global quenches. We observe different relaxation behaviors for quenches starting from different initial states to the critical point. Starting from a fully ordered initial state, the relaxation is logarithmically slow described by $\overline{m}(t) \sim \ln^a t$, and in a finite sample of length $L$ the average magnetization saturates at a size-dependent plateau $\overline{m}_p(L) \sim L^{-b}$; here the two exponents satisfy the relation $b/a=\psi=1/2$. Starting from a fully disordered initial state, the magnetization stays at zero for a period of time until $t=t_d$ with $\ln t_d \sim L^{\psi}$ and then starts to increase until it saturates to an asymptotic value $\overline{m}_p(L) \sim L^{-b'}$, with $b'\approx 1.5$. For both quenching protocols, finite-size scaling is satisfied in terms of the scaled variable $\ln t/L^{\psi}$. Furthermore, the distribution of long-time limiting values of the magnetization shows that the typical and the average values scale differently and the average is governed by rare events. The non-equilibrium dynamical behavior of the magnetization is explained through semi-classical theory.

In a system with spatial inhomogeneities, such as defects [70,71,72,73,74,75,76], quasi-periodic or aperiodic interactions [77,78,79], the non-equilibrium relaxation process becomes qualitatively different from the homogeneous case. In the semiclassical picture the quasi-particles are also inhomogeneously created and they move in a complex diffusive way. With randomness the eigenstates of the Hamiltonian can be localized, which leads to Anderson localization [80] in non-interacting systems or many-body localization [81] if the particles are interacting; in these cases the quasiparticles generally move only finite distances away from their place of creation or follow some ultra-slow diffusive motion. As a result the non-equilibrium long-time stationary state of disordered systems retains memory of the initial state and therefore thermalization does not take place [82].
Concerning the functional form of the relaxation process after a quench in random quantum systems, there have been detailed studies about the time-dependence of the entanglement entropy [83,84,85,86,87]. If the system consists of non-interacting fermions -such as the critical XX-spin chain with bond disorder or the critical random transverse-field Ising chain -the dynamical entanglement entropy grows ultraslowly in time as and saturates in a finite system at a value where ℓ denotes the size of a block in a bipartite system and can be chosen to be proportional to the size of the system L [84,87]. These scaling forms can be explained by a strong disorder renormalization-group (SDRG) approach [88]. Recently, the SDRG method, which was designed as a ground state approach, has been generalized to take into account excited states [89,90,91]; this generalized RG method is often abbreviated as RSRG-X [89]. By this generalized SDRG method the ratio of the prefactors in (1) and (2) is predicted as b/a = ψ ne , where ψ ne = 1/2 is a critical exponent in the non-equilibrium process and describes the relation between time-scale and length-scale as For interacting fermion models due to many-body localization the time-dependence of the dynamical entropy is S(t) ∼ ln ω t with ω ≥ 1, while the saturation value follows the volume law, S(ℓ) ∼ ℓ [86].
In the present paper we study the relaxation of the order parameter of the random transverse-field Ising chain after a quench to the critical final state and into a ferromagnetic state. We consider relaxation processes from an initial ferromagnetic state and from a fully paramagnetic state by a sudden change of the strength of the transverse field. To circumvent numerical instability as observed in previous calculations for large systems using eigenvalue solver routines [84,87], we use multiple precision arithmetic to study the time-evolution through direct matrix multiplications.
The rest of the paper is organized as follows. The model and the method for calculating the local magnetization are described in section 2 . The numerical results for relaxation of the magnetization following different quench protocols are presented and discussed in section 3. A summary is given in section 4. Details of the time evolution of Majorana fermion operators are given in the appendix.

The model and the method
The model we consider is the random transverse-field Ising chain of length L defined by the Hamiltonian: in terms of the Pauli matrices σ x,z i at site i. In this paper we will consider open chains with free boundary conditions. The couplings, J i , and the transverse fields, h i , are position dependent random numbers taken from the uniform distributions in the intervals [0, 1] and [0, 1]h, respectively. The strength of the random transverse field, h, is time-dependent; for t < 0 its value is h 0 and for t > 0 it changes suddenly to h ( = h 0 ). The initial and the final Hamiltonians are denoted by H 0 and H, respectively. The system for t < 0 is prepared in the ground state |Ψ (0) 0 of the initial Hamiltonian; after the quench for t > 0 the system evolves according to the final Hamiltonian H and the state of the system is described by |Ψ 0 (t) = exp(−ıHt)|Ψ (0) 0 , which is generally not an eigenstate of H. Here and throughout the paper we denote the imaginary unit √ −1 by ı to avoid confusion with the integer index i and set = 1. To calculate the time-dependent expectation value A(t) of an observableÂ, we work in the Heisenberg picture usingÂ H (t) = exp(ıHt)Â exp(−ıHt) and evaluate A(t) = Ψ  The standard way to deal with the Hamiltonian of the transverse field Ising chain is the mapping to spinless free fermions [92,93]. The spin operators σ x,y,z i are expressed in terms of fermion creation (annihilation) operators c † i (c i ) by using the Jordan-Wigner transformation [94]: where a ± j = (σ x j ± ıσ y j )/2. The Ising Hamiltonian in (4) can then be written in a quadratic form in fermion operators: In this paper we are interested in the relaxation of the local order parameter (magnetization), m l (t), at a position l. Following the method by Yang [95], for a given sample of large length L we define m l (t) as the off-diagonal matrix element of the longitudinal magnetization operator: m l (t) = Ψ is the first excited state of H 0 . In the free-fermion representation we introduce two Majorana fermion operators,ǎ 2i−1 andǎ 2i at site i with the definition: which obey the commutation relations: The spin operators are then expressed in terms of the Majorana operators as: The calculation of the magnetization, m l (t), involves the evaluation of the determinant of a 2l × 2l antisymmetric matrix C with the elements being the correlation functions , and C ji = −C ij . For details see the appendix of Ref. [78]. For a random system any physical quantity requires an averaging over disorder realizations. We denote the disorder average by an overbar; accordingly, the average local magnetization is expressed as

In equilibrium
In equilibrium (i.e. for t < 0), the critical behavior of the local magnetization m l (0) has been analytically studied by the SDRG method [96], and the random quantum critical point (at h 0 = h c = 1) is found to be controlled by an infinite-randomness fixed point [97,88]. In the thermodynamic limit (L → ∞) one needs to discriminate between bulk (l/L = O(1)) and surface (l = O(1)) points, and the corresponding average magnetization is denoted by m and m s , respectively. In the ordered phase, The surface magnetization follows similar behavior, however with a surface exponent β s = 1. At the critical point of a finite system the scaling of the average magnetization is governed by rare samples (or rare regions), which have the magnetization of order of O(1), whereas in typical samples the magnetization behaves as m typ ∼ exp(−A √ L). The fraction of rare events scales as P rare ∼ L −x , and the same holds for the average magnetization: m ∼ L −x ; here the exponent x = β/ν is related to the critical exponent ν = 2 of the average correlation function. Similarly, the scaling behavior of the surface magnetization involves the exponent x s = β s /ν. Concerning equilibrium dynamical scaling at the infinite-randomness fixed point, it is extremely space-time anisotropic so that the typical length, ξ, and the typical time, τ , is related as ln τ ∼ ξ ψ , (10) with an exponent ψ = 1/2.

Out-of-equilibrium
Some out-of-equilibrium properties of the random transverse-field Ising chain have been predicted by the RSRG-X approach [89,90,91]. For a quench to the critical state, i.e. h = h c in the final Hamiltonian, the relation between the typical length scale and the typical time scale is in the same form as (10) in the equilibrium case and even the corresponding exponents are the same ψ = ψ ne = 1/2 (cf. (3) and (10)). However, available numerical results [84,87] for the relaxation of the entanglement entropy are not in complete agreement with the RSRG-X prediction.
Here we focus on the relaxation of the magnetization in the random chain, which to our knowledge has not been previously studied. The crucial point of the method is the calculation of the time evolution of the Majorana operators introduced in (6), which can be given in a closed form, provided the Hamiltonian is expressed in a quadratic form in terms of fermion operators, like the form in (5). The calculation for a general quadratic Hamiltonian is described in the appendix. Definingǎ(t) as a vector with componentsǎ m (t) for m = 1, 2, . . . , 2L, we can express the time evolution of the Majorana operators aš where the 2L × 2L matrix P(t) is given by To evaluate the matrix exponential in (12), one can use spectral decomposition of M by diagonalizing the large 2L × 2L matrix. For disordered systems, because of some extremely small eigenvalues standard eigenvalue solvers would fail to converge for some large-size samples, leading to significant numerical errors. This problem was observed in our preparatory work. Therefore, we reformulated our numerical procedure to avoid using any eigenvalue solver routine and instead solve the time-evolution problem by matrix multiplication using multiple precision arithmetic. In our numerical procedure we first evaluated the matrix exponential at a unit time step, t st = 1, using the Taylor expansion: and calculated the sum of the first hundred terms with multiple precision arithmetic; the absolute value of the last term of this sum is less than 2 100 /100! ≈ 1.3 × 10 −128 since the eigenvalues of M are in the range [−2, 2] [98]. The truncation error in (14) is therefore sufficiently small even for octuple precision. For larger time steps we used the identity: P(t 1 )P(t 2 ) = P(t 1 + t 2 ) and iterated it with t 1 = t 2 , starting from t 1 = t st = 1 until a given time step t n = 2 n . We have checked the accumulation of errors by the condition that the matrix P should be orthogonal. While calculating the dot products of the arrays of P, the off-diagonal products have to be zero. The sum of the absolute value of the numerically calculated off-diagonal dot products are found to be less than 10 −30 , even after n = 250 iterations, which represents the longest time in our calculation.

Relaxation of the magnetization
In this chapter we present our numerical results for the relaxation of the magnetization of the random transverse-field Ising chain from two different initial states. In the first part the initial state is fully ordered with h 0 = 0 for the parameter of the transverse field; in the second part we consider a fully disordered initial state with 1/h 0 = 0. For each system sizes about 10, 000 samples were considered to obtain the disorder average. We have used doubledouble or quadruple precision arithmetic for our numerical study and also in some cases checked the accuracy by comparing the results in the large-time limit with those obtained with octuple precision.

Relaxation from a fully ordered initial state
We first consider quenches starting from a fully ordered initial state with h 0 = 0.
3.1.1. Ferromagnetic final state Before we study quenches to the critical point, we first consider the situation in which the system is quenched to a state within the ferromagnetic state with h < 1. In figure 1(a) we present the time dependence of the average magnetization after a quench to h = 0.5 for different chain lengths from L = 16 to L = 128. As seen in this figure, there is first a decay, followed by a plateau extending to time t d , after which (t > t d ) the magnetization decays sharply before it saturates at a second plateau; the delay time t d is L-dependent and scales approximately as ln(ln t d ) ≈ ln L + const, as shown in the inset of figure 1(a). Similar characteristics is observed for other values of h < h c (shown in figure 1(b)) when the final state is within the ferromagnetic phase.
To explain the behavior observed in figure 1 we recall that according to a semiclassical theory [36,18,41,42] the local magnetization at a given site decreases in time if an odd number of quasi-particles, created by the quench, pass through the site. In a random chain the motion of the quasi-particles is limited within a finite localization length ξ loc (h), which depends on h and grows as the critical point is approached; thus the decay of the magnetization stops at some point of time when the travel distance of the quasi-particles is reached, which explains the presence of the first plateau and why this plateau is smeared out when h approaches the critical point ( figure 1(b)). The second decay of the magnetization is related to the lowest excitation. In the ferromagnetic phase, the lowest energy gap vanishes exponentially with the system size as ǫ 0 (L) ∼ exp(−L/ξ), ξ being the correlation length, thus the corresponding time scale t d ∼ 1/ǫ 0 grows exponentially with L, as shown in the inset of figure 1(a). At the critical point, the lowest excitation vanishes as ǫ 0 (L) ∼ exp(−AL 1/2 ) (see (10)); thus, different relaxation behavior is expected for a quench to the critical point and this is discussed in the following subsection.

Critical final state
Relaxation of the average magnetization from the fully ordered initial state to the critical point is shown in figure 2 for various chain lengths.
Here the first plateau observed in a ferromagnetic final state is missing, and instead there is a transient region showing a linear relation between ln m(t) and ln(ln t) with a slope a ≈ 0.14, corresponding to a logarithmically show decay: collapse. Here we comment on the logarithmic decay described in (15). A similar logarithmic decay is also present in the XX-chain with bond disorder, where a larger exponent a = 2 is found [90,87]. The critical states of the random-bond XX-chain and the random transverse-field Ising chain are both governed by an infinite-randomness fixed point in SDRG language. The slower decay (corresponding to a smaller exponent a) for the random transverse-field Ising chain is due to effective spin clusters formed during renormalization. In the semi-classical picture [42], relaxation of the order parameter is contributed by quasi-particles created during the quench; a spin is flipped if the trajectory of a quasi-particle crosses the given site at a later time. For the relaxation of a larger cluster more quasi-particles are needed, which takes more time.
The magnetic moment of an effective spin cluster of the random transverse-field Ising chain scales as µ ∼ L d f , with a fractal dimension d f > 0. The process of formation of spin clusters is missing in the SDRG procedure for the XX-chain, where only random singlets are formed, corresponding to d f = 0. The different SDRG procedures for the Ising chain and the XX-chain lead to different functional forms of the decay. In a random system the distributions of some observables are also of interest. Here we focus on the distributions of the large-time limiting values m p of the magnetization. In figure 3(a) we plot the distribution of the logarithmic values ln(m p ) for t exp(exp(4.5)) for different system sizes. The distribution of the logarithmic variable is broad and becomes broader with increasing L. A good scaling collapse can be obtained in terms of the scaled variable y = ln(m p )L −α , thus determined by the behavior of the distribution function in the limit of y → 0 − , which is assumed to be in a power-law form: with χ = 0 in our case shown in figure 3(b). The average of m p is then given by thus, with α = b and χ = 0 we recover the relation m p ∼ L −b , as given in (16).

Relaxation from a fully disordered initial state
Now we turn to the case in which the system is quenched from a fully disordered initial state (h 0 = ∞) to the critical point (h = 1). Looking at the time-dependent average magnetization for various system sizes plotted in figure 4(a), there is a period of time 0 < t < t d (L) just after the quench where the magnetization stays negligible and practically corresponds to the initial magnetization before the quench. After this delay time t d , which is L-dependent, the average magnetization starts to increase rapidly and in the large-time limit it approaches a plateau, the value of which has a power-law L-dependence: with b ′ = 1.46, as measured in the inset of figure 4(a). Using the scaled variables mL b ′ and ln t/L ψ , we achieve a good data collapse for the time-dependent magnetization in the scaling plot shown in figure 4(b). The increase of the average magnetization in the disordered chain after the quench is unexpected since semiclassical theory predicts a monotonous decrease of the dynamical magnetization in a homogeneous system. To understand the mechanism behind the increase of the average magnetization, we have studied the distribution of the large-time limiting values m p (L), which is presented in figure 5(a). The distribution is extremely broad even on the logarithmic  scale and it gets broader with increasing system size; in terms of a scaling combination y = ln m p L −α ′ with α ′ = 0.5, the inset of figure 5(a) shows a scaling plot with good data collapse. We can then conclude that the typical value scales with L as m typ p ∼ exp(−C L α ′ ) and is much smaller than the average given in (20). Thus also in this quench process the average value is determined by rare events, in which the largetime limiting values of the magnetization are m p (L) = O(1). Similar to the situation for a quench from an ordered state as discussed in section 3.1, we have found a powerlaw formP (y) ∼ (−y) χ ′ for y → 0 − with an estimated exponent χ ′ = 2.0 (shown in figure 5(b)). Thus, we have found the size-dependence of the average magnetization in (19) as m p (L) ∼ L −α ′ (1+χ ′ ) = L −1.5 , which agrees with (20) within the error bars.
The scaling behavior of the large-time limiting value of the dynamical order parameter can be better understood if the boundary site (l = 1) of an open chain is considered. The large-time limiting value of the surface magnetization, m s p , can be exactly calculated (see equation (16) in Ref. [78]): Here: For Φ where Φ 1 (1) = m s eq is the equilibrium value of the surface magnetization of the chain [100,101], evaluated in the final state, i.e. with h = h c . It is known that the typical value of m s eq at the critical point scales as exp(−CL 1/2 ) [101, 102] and the same scaling form holds for m s p ; thus the proper scaling combination for m s p is y = ln(m s p )L −1/2 , which in the same form for bulk spins discussed above. Concerning the scaling behavior of the average value m s p , it is dominated by rare events. As described in Ref. [101], there is a close analogy between the random transversefield Ising chain and a one-dimensional random walk: to a given sample with a set of couplings J i , and transverse fields h i , i = 1, 2, . . . L, one can assign a onedimensional random walk which starts at the origin and takes consecutive steps, the length of the i-th step being ln(h i /J i ). For a rare realization of m s eq , the associated random walk has a surviving character, i.e. it stays at positive position for all steps. Concerning m s p , here for a rare realization both m s eq and the product lm−1 j=1 hj Jj should be of order of O(1); in the language of random walks, this means that the walk is surviving and returns after l m steps. This particular event has the probability: P(l m ) ∼ l −3/2 m (L − l m ) −1/2 and its average value scale as: lm P(l m )/L ∼ L −3/2 , which means that m s p (L) ∼ L −3/2 . The numerically observed scaling behavior for the bulk magnetization in (20) is similar to this form.

Conclusions
In this paper we have studied numerically the time dependent magnetization of the random transverse-field Ising chain after a quench. In order to obtain accurate numerical results we have mainly performed matrix product operations with multiple precision arithmetics in our calculations, instead of any eigenvalue solver routine. In this way we obtained our finite-size results up to length L = 256 free from numerical instability. We note that some problems with numerical instability using eigenvalue solvers occurred when studying the entanglement entropy in large systems like L = 256 [84,87].
For a quench from a fully ordered initial state to a state within the ferromagnetic phase, the average magnetization is shown to relax first to a plateau value, followed by a second relaxation after a delay time t d ∼ exp(cL) to a second plateau value. This second relaxation is attributed to quasi-localized modes which are present in finite systems in the ferromagnetic phase. Both plateau values are h-dependent, and finite because in the free-fermionic picture the emitted quasi-particles are localized and travel only a finite distance, which reduces the initial order parameter by a finite fraction only. This behavior is similar to that observed in the Aubry-André model when the quench is performed to the localized phase [79].
When the system is quenched from a fully ordered initial state to the critical finale state, the relaxation of the average magnetization is logarithmically slow (see (15)), and the large-time limiting value decays as a power of L (see (16)). Between the time scale (t) and the length scale (L) we have found the same relation as in equilibrium, ln t ∼ √ L, as predicted by the RSRG-X approach [89,90,91]. We have also studied the distribution of large-time limiting values of the magnetization: the typical value decreases exponentially with the system size as exp(−CL α ), more rapidly than the average value, which decays as a power-law and is determined by rare realizations. The relaxation process after the quench can be explained by the diffusion of quasiparticles in a random environment; the travel distance (ℓ) of the quasi-particles within time t scales like ln(t) ∼ √ ℓ. In a quench process from a fully disordered initial state to the critical point, one observes a delay time, t d ∼ exp(CL 1/2 ), in which the average magnetization is negligible, and for t > t d a rapid increase of the average magnetization toward an asymptotic large-time limiting value, which has a power-law L-dependence (see (20)). The distribution of the large-time limiting magnetizations shows that the typical and the average values scale differently, and the average is determined by rare events. Qualitatively and even quantitatively similar behavior is found for the time dependence of the surface magnetization, which has exact scaling arguments. It has been explicitly shown that there are rare samples, in which after a delay time a stationary surface magnetization of order of O(1) develops. These rare samples will then dominate the average value. Such a phenomenon is similar to phase ordering dynamics in classical systems [103], but not expected in closed homogeneous quantum systems.  with the initial condition P mn (0) = δ mn . The solution as a shorthand matrix notation is given in (12), in which the exponential can be evaluated by using the spectral decomposition of M. Indeed, the eigenvalues of M are given by the energies of the freefermionic modes of H in (A.1), which can be seen by forming M 2 and comparing it with the results in Ref. [92] . Details of the calculation using the spectral decomposition of M are presented in the appendix of Ref. [78].