Quantum quench dynamics in the transverse-field Ising model: A numerical expansion in linked rectangular clusters

We study quantum quenches in the transverse-field Ising model defined on different lattice geometries such as chains, two- and three-leg ladders, and two-dimensional square lattices. Starting from fully polarized initial states, we consider the dynamics of the transverse and the longitudinal magnetization for quenches to weak, strong, and critical values of the transverse field. To this end, we rely on an efficient combination of numerical linked cluster expansions (NLCEs) and a forward propagation of pure states in real time. As a main result, we demonstrate that NLCEs comprising solely rectangular clusters provide a promising approach to study the real-time dynamics of two-dimensional quantum many-body systems directly in the thermodynamic limit. By comparing to existing data from the literature, we unveil that NLCEs yield converged results on time scales which are competitive to other state-of-the-art numerical methods.


Introduction
Understanding the dynamics of isolated quantum many-body systems out of equilibrium is an active area of research of modern theoretical and experimental physics [1][2][3]. A popular nonequilibrium protocol in this context is a so-called quantum quench [4]. In such quench protocols, the system's Hamiltonian H depends on some parameter λ, and the system is prepared in an eigenstate |ψ(0) of H, e.g., the groundstate, for an initial value λ i . Next, the value of λ is suddenly changed, λ i → λ f , such that |ψ(0) is no eigenstate of H(λ f ), and the system exhibits nontrivial dynamics. For an isolated quantum system undergoing unitary time evolution, it is then intriguing to study if and in which way the system relaxes back to equilibrium. Central questions are, for instance, how the (short-or long-time) dynamics can be described in terms of "universal" principles [1,[5][6][7][8][9][10][11], what are the relevant time scales of relaxation [12][13][14], and whether or not the long-time values of physical observables agree with the prediction of, e.g., a microcanonical or canonical ensemble (i.e. thermalization) [15][16][17].
One possible mechanism to explain the emergence of thermalization in isolated quantum systems is given by the eigenstate thermalization hypothesis (ETH) [18][19][20]. While the validity of the ETH has been numerically tested for a variety of models and observables (see, e.g., [21][22][23][24][25][26][27]), there are also classes of systems which violate the ETH and fail to thermalize. One such class is given by integrable models, where the extensive number of conservation laws prevents the applicability of standard statistical ensembles [28]. Instead, it has been proposed that integrable models equilibrate towards a generalized Gibbs ensemble (GGE), which maximizes the entropy with respect to the conserved charges [29][30][31]. In addition, it is now widely believed that some strongly disordered systems can undergo a transition to a many-body localized (MBL) phase, where the ETH is violated as well [32,33]. Moreover, there has been plenty of interest recently in models which are, in a sense, intermediate cases between "fully ETH" and "fully MBL". This includes, e.g., models featuring "quantum scars" where rare ETH-violating states are embedded in an otherwise thermal spectrum [34][35][36][37][38], as well as models which exhibit a strong fragmentation of the Hilbert space due to additional contraints [39,40].
From a numerical point of view, studying the nonequilibrium dynamics of isolated quantum many-body systems is a challenging task. This is not least caused by the fact that for an interacting quantum system, the Hilbert space grows exponentially in the number of constituents. Nevertheless, thanks to the continuous increase of computational resources and the development of sophisticated numerical methods including, e.g., dynamical mean field theory [41], Krylov subspace techniques [42,43], dynamical quantum typicality [44], or classical representations in phase space [45], significant progress has been made. Especially for onedimensional systems, the time-dependent density-matrix renormalization group, as well as related methods based on matrix-product states (MPS), provide a powerful tool to study dynamical properties for system sizes practically in the thermodynamic limit [46,47]. However, since these methods rely on an efficient compression of moderately entangled wavefunctions, the reachable time scales in simulations are eventually limited due to the inevitable buildup of entanglement during the unitary time evolution.
The growth of entanglement becomes even more severe in spatial dimensions larger than one. Despite recent advances involving MPS-based or tensor-network algorithms [48][49][50][51][52][53], as well as the advent of innovative machine-learning approaches [54][55][56], the time scales numerically attainable for two-dimensional quantum many-body systems are still comparatively short. While the dynamics of such two-dimensional systems can nowadays be accessed in experiments with quantum simulators [57][58][59], the development of efficient numerical techniques is paramount.
In this paper, we scrutinize the nonequilibrium dynamics for quantum quenches in the Ising model with transverse magnetic field. While this model is exactly solvable in the case of a chain and has been studied in numerous instances, our main focus is on nonintegrable geometries such as two-and three-leg ladders and, in particular, two-dimensional square lattices. To this end, we rely on an efficient combination of numerical linked cluster expansions (NLCEs) and the iterative forward propagation of pure states in real time via Chebyshev polynomials. Depending on the model geometry, the initial state, and the strength of the quench, the nonequilibrium dynamics is found to display a variety of different behaviors ranging from rapid equilibration, over slower monotonous relaxation, to persistent (weakly damped) oscillations. Most importantly, from a methodological point of view, we demonstrate that NLCEs comprising solely rectangular clusters provide a promising approach to study the real-time dynamics of two-dimensional quantum many-body systems directly in the thermodynamic limit. By comparing to existing data from the literature, we unveil that NLCEs yield converged results on time scales which are competitive to other state-of-the-art numerical methods. This paper is structured as follows. In Sec. 2, we introduce the models, observables, and quench protocols which are studied. In Sec. 3, we then discuss the employed numerical methods, while our results are presented in Sec. 4. We summarize and conclude in Sec. 5.

Models, observables, and quench protocols
We study the Ising model with ferromagnetic nearest-neighbor interactions and transverse magnetic field, described by the Hamiltonian where the first sum on the right hand side runs over all pairs of nearest neighbors and m, L is the total number of sites, J > 0 sets the energy scale, g > 0 denotes the strength of the transverse field, and σ x,z are Pauli matrices at site . Note that the Hamiltonian (1) is symmetric under the global spin-flip operation σ z → −σ z . In this paper, the transverse-field Ising model (1) is considered for different lattice geometries such as chains (L = L x ), two-and three-leg ladders (L = L x × 2, L = L x × 3), and two-dimensional square lattices (L = L x × L y ). While we generally intend to obtain results in the thermodynamic limit L → ∞ (see Sec. 3.1 for our numerical approach), we consider finite system sizes as well. In the case L < ∞, one has to distinguish between open boundary conditions (OBC) and periodic boundary conditions (PBC), where for chains and ladders the latter only applies in the x direction.
On the one hand, in the case of a chain, H is a paradigmatic example of an integrable model and can be solved exactly by subsequent Jordan-Wigner, Fourier, and Bogolioubov transforms [60], see also Appendix A. For g < 1, H is in a ferromagnetic phase with a twofold degenerate groundstate. At the critical point g = 1, H undergoes a quantum phase transition towards a paramagnetic phase with unique groundstate for g > 1. On the other hand, for a two-dimensional square lattice, H is nonintegrable [24,25,61], and the quantum phase transition between an ordered phase and an unordered phase occurs at the larger transverse field g = g c ≈ 3.044 [62]. For intermediate cases, such as "semi-infinite" multi-leg ladders on a cylinder geometry, the exact value of g c can vary due to finite-size corrections [50].
In this paper, we consider quench protocols starting from fully polarized initial states |ψ(0) . Namely, we either study quenches starting from |ψ(0) = |↑ , where all spins are initially aligned along the z axis, or quenches starting from the state |ψ(0) = |→ , where all spins point in the x direction. Note that written in the common eigenbasis of the local σ z , |→ is a uniform superposition of all 2 L basis states. Moreover, while the state |↑ is an eigenstate of H for vanishing field g = 0, the state |→ is the groundstate of H for g → ∞.
Given the states |↑ and |→ , we study the nonequilibrium dynamics resulting from quantum quenches to weak (g < g c ), strong (g > g c ), or critical values (g = g c ) of the transverse field, i.e., depending on the initial state these are quenches either within the same equilibrium phase, or to or across the critical point. Due to the quench, the fully polarized states |↑ and |→ are no eigenstates of H anymore and evolve unitarily in time ( = 1), Consequently, the expectation values of observables acquire a dependence on time as well. In particular, we here consider the dynamics of the transverse and the longitudinal magnetization,

Numerical approach
We now discuss the numerical methods which are employed in this paper. Throughout this section, we exemplarily focus on the transverse magnetization X(t) . The calculations for Z(t) are carried out analogously.
The main idea of NLCEs is that the per-site value of an extensive quantity in the thermodynamic limit can be obtained as a sum over contributions from all linked clusters which can be embedded on the lattice [77], where the sum runs over all connected clusters c with multiplicities L c and weights W c (t).
The notion of a connected cluster here refers to a finite number of lattice sites where every site of the cluster has to be connected to at least one other cluster site. Given a cluster c, its weight W c (t) is obtained by an inclusion-exclusion principle. That is, the quantity of interest (here the dynamics of the magnetization X) is evaluated on the cluster c (with OBC) and, subsequently, the weights W s (t) of all subclusters s of c have to be subtracted [77], While NLCEs yield results in the thermodynamic limit (such that a finite-size scaling becomes unnecessary), it is instead crucial to check the convergence of the series. To this end, the sum in Eq. (6) is usually organized in terms of expansion orders [77]. For instance, one could group together all clusters which comprise a certain number of lattice sites. Then, an expansion up to order C refers to the fact that all clusters with up to C lattice sites are considered in Eq. (6). Moreover, the NLCE is said to be converged if the outcome of Eq. (6) does not depend on the value of C.
At this point, it is important to note that in actual simulations, the maximum order C that can be reached is limited by two factors: (i) the exponential growth of the Hilbert-space dimension with increasing cluster size, and (ii) the necessity to identify the (possibly very large number of) topologically distinct clusters and to calculate their weights. Since a larger expansion order typically leads to a convergence of Eq. (6) up to longer times [75] (or down to lower temperatures for thermodynamic quantities [66]), it is desirable to include clusters as large as possible. In this paper, we therefore aim to mitigate the limitations (i) and (ii) by two complementary approaches. First, instead of using exact diagonalization to evaluate X(t) (c) , we here employ an efficient forward propagation of pure states (see Sec. 3.2), which is feasible for significantly larger Hilbert-space dimensions. Secondly, in order to reduce the enormous combinatorial costs to generate (and evaluate) all clusters with a given number of sites, we rely on the fact that the sum in Eq. (6) can also converge for different types of expansions, as long as clusters and subclusters can be defined in a self-consistent manner [77]. In particular, we here restrict ourselves to only those clusters which have a rectangular shape [68]. Due to this restriction, the number of topologically distinct clusters is significantly reduced and the calculation of weights W c (t) becomes rather simple since all subclusters are rectangles as well, (c) For the two-dimensional square lattice, we restrict ourselves to clusters with a rectangular shape. A cluster c = (x, y) with x > y is equivalent to its 90 • -rotated counterpart c = (y, x), such that we only need to consider c which enters Eq. (6) with L c = 2. For square-shaped clusters with x = y, we have L c = 1.
see Fig. 1. Specifically, given a cluster c = (x, y) of width x and height y, its weight W (x,y) (t) follows as [78] W Note that in practice, we only need to consider clusters with x ≥ y. On the one hand, for a cluster c = (x, y) with x > y, the cluster c = (y, x), i.e., c rotated by 90 degrees, is topologically equivalent, such that c simply enters the sum (6) with a multiplicity L c = 2. On the other hand, for square-shaped clusters with x = y, we have L c = 1. Note further that for the case of a chain, Eq. (6) reduces to a single difference between X(t) (c) evaluated on the largest and the second-largest cluster [75]. The same holds true for two-leg (or three-leg) ladders if only clusters with y = 2 (or y = 3) are considered.

Pure-state propagation
Evaluating the unitary time evolution of the initial states |ψ(0) according to Eq. (4) in principle requires the exact diagonalization (ED) of the Hamiltonian H. In order to access system (and cluster) sizes beyond the range of ED, we here subdivide the evolution up to time t into a product of discrete time steps, where δt = t/Q. If the time step δt is chosen sufficiently small, then there exist various approaches to accurately approximate the action of the exponential exp(−iHδt) such as, e.g., Trotter decompositions [79], Krylov subspace techniques [42], or Runge-Kutta schemes [80,81].
In this paper, we rely on an expansion of the time-evolution operator in terms of Chebyshev polynomials, for a comprehensive overview see [82][83][84][85].
Since the Chebyshev polynomials are defined on the interval [−1, 1], the spectrum of the original Hamiltonian H has to be rescaled [85], where a and b are suitably chosen parameters. In practice, we use the fact that the (absolute of the) extremal eigenvalue of H can be bounded from above according to [83] max where E max (E min ) is the largest (smallest) eigenvalue of H, and N ,m denotes the number of nearest-neighbor pairs , m , i.e., the number of bonds of the lattice. By choosing a ≥ E, it is guaranteed that the spectrum of H lies within [−1, 1]. As a consequence, we can set b = 0. Note that while this choice of a and b is not necessarily optimal, it proves to be sufficient [83] (see also Appendix B).
Within the Chebyshev-polynomial formalism, the time evolution of a state |ψ(t) can then be approximated as an expansion up to order M [85], where the expansion coefficients c 0 , c 1 , . . . , c M , are given by with J k (aδt) being the k-th order Bessel function of the first kind evaluated at aδt. [Note that the notation in Eqs. (12) and (13) assumes b = 0.] Moreover, the vectors |v k are recursively generated according to with |v 1 = H |v 0 and |v 0 = |ψ(t) . Given a time step δt (and the parameter a), the expansion order M has to be chosen large enough to ensure negligible numerical errors. In this paper, we typically have δtJ = 0.02 and M = 15, which turns out to yield very accurate results (see Appendix B). As becomes apparent from Eqs. (12) and (14), the time evolution of the pure state |ψ(t) requires the evaluation of matrix-vector products. Since H is a sparse matrix, these matrixvector multiplications can be implemented comparatively time and memory efficient. In particular, we here calculate the matrix elements of H on the fly and use parallelization to reduce the runtime. Thus, the memory requirements are essentially given by the size of the state |ψ(t) and the auxiliary states |v k−1 , |v k , and |v k+1 . As a consequence, it is possible to treat system (or cluster) sizes significantly larger compared to ED (here up to 28 lattice sites with a Hilbert-space dimension of d ≈ 10 8 ). Since the transverse-field Ising model (1) does not conserve the total magnetizations X or Z, the corresponding quantum numbers cannot be used to block-diagonalize H. Moreover, the clusters entering the NLCE are defined with open boundary conditions such that translational invariance cannot be exploited.

Results
We now present our numerical results for the quench dynamics of X(t) and Z(t) in chains, ladders, and two-dimensional lattices. Our main focus is to analyze the convergence properties of the NLCE by comparing to direct simulations of finite systems with periodic boundary conditions and to existing data from the literature.

Chains
The transverse-field Ising chain is a paradigmatic example of an exactly solvable model and analytical solutions have been known for a long time [60,[86][87][88] (see also Appendix A). Since quantum quenches in the Ising chain have been studied extensively before (see, e.g., Refs. [89][90][91][92][93][94][95][96]), the present section should be mainly understood as a consistency check for our numerical methods and a preparation for the study of ladders and two-dimensional lattices in Secs. 4.2 and 4.3. (It might be fair to say, however, that explicit visualizations of the analytical solutions, e.g., for the full time-dependent relaxation process of X(t) for specific initial states and transverse fields g, are less often available in the literature.) In Figs. 2 (a)-(c), the dynamics of the transverse magnetization X(t) is shown for quenches starting from the initial state |ψ(0) = |↑ and different values of the transverse field g = 0.5, 1, 2. (Recall that the quantum critical point is g = 1 for the chain geometry.) Numerical data obtained by NLCE for expansion orders C = 24, 25 are compared to (i) a simulation for a finite chain with L x = 25 and PBC, and (ii) the exact, analytically known result [see Eq. (16) in Appendix A]. Starting from its initial value X(0) = 0, we observe that X(t) quickly increases and exhibits a peak at short times, before equilibrating towards a constant long-time value. This stationary value is reached already for times tJ ≈ 2. While this overall behavior of X(t) is very similar for all values of g considered, the long-time value X(t → ∞) is found to vary with g. In particular, it is known that this long-time value can be described in terms of a suitable GGE [28].
Generally, we find that the NLCE results in Figs. 2 (a)-(c) are well converged on the time scales depicted, i.e., the curves for expansion orders C = 24 and C = 25 agree convincingly with each other. (Note that for g = 2 and tJ ≈ 10, a first indication for a breakdown of convergence can be seen.) In addition, the curves for the finite chain with L x = 25 in Figs. 2 (a)-(c) also coincide nicely with the NLCE data for L → ∞. Thus, finite-site effects appear to be less relevant in these cases. Moreover, our numerical results for X(t) agree perfectly with the analytical solution.
Next, in Figs. 2 (d)-(f), we consider quenches starting from the state |ψ(0) = |→ . Despite the obvious difference that X(t) now starts at a maximum, X(0) = 1, the general picture is very similar compared to the previous case of |ψ(0) = |↑ . Namely, X(t) exhibits a rapid decay and equilibrates rather quickly towards its long-time value. Especially for g = 1 [ Fig. 2 (e)], however, we now observe pronounced finite-size effects, i.e., the curve for L x = 25 deviates from the analytical solution for times tJ 5 and exhibits oscillations. In contrast, the NLCE results for C = 24, 25 remain converged up to at least tJ = 10. This is a remarkable result since the largest cluster in the NLCE also only has 25 lattice sites, i.e., the computational complexities of the NLCE and the simulation of the finite system are essentially the same.
Depending on the details of the quench, we thus find that performing a NLCE can yield a numerical advantage over the direct simulation of finite systems. On the one hand, if finitesize effects are weak, the results for finite chains can be very similar to the actual L → ∞ dynamics (and also remain meaningful on longer time scales where the NLCE breaks down). On the other hand, the presence of strong finite-size effects [e.g. at the quantum critical point, cf. Fig. 2 (e)] appears to favor the usage of NLCEs which yield the dynamics directly in the thermodynamic limit. This is a first result of the present paper. As will be discussed in more detail in the upcoming sections, a similar parameter-dependent advantage (or disadvantage) of performing a NLCE occurs for ladder geometries and two-dimensional lattices as well.

Ladders
Let us now turn to the results for two-and three-leg ladders, which can be seen as intermediate cases between the chain geometry (cf. Sec. 4.1) and the two-dimensional square lattice (cf. Sec. 4.3). Since exact solutions for the dynamics of ladders are absent, we cannot compare our numerical data to analytical results. (For additional remarks on the transition from integrability to nonintegrability, see also Appendix C.) In Fig. 3, we consider quenches starting from the state |ψ(0) = |↑ in two-leg ladders with different transverse fields g. Here, the data is obtained by NLCE for expansion orders C = 24 and C = 26, i.e., the largest clusters involved are of size 12 × 2 or 13 × 2. As shown in Fig.  3 (a), the dynamics of the longitudinal magnetization Z(t) displays a strong dependence on the value of g. On the one hand, for g = 2, Z(t) rapidly decays, exhibits a minimum at tJ ≈ 1, and equilibrates to zero for tJ 3. On the other hand, for g = 1, the decay of Z(t) towards zero is distinctly slower and much more monotonous. Moreover, for g = 0.5 (i.e. a quench within the same equilibrium phase), the decay of Z(t) is almost indiscernible on the time scale shown, and we additionally observe that Z(t) exhibits small oscillations for this value of g. The corresponding dynamics of the transverse magnetization X(t) is shown in Fig. 3 (b). While X(t) quickly equilibrates towards a stationary value for g = 2, X(t) displays essentially undamped oscillations for g = 0.5 and does not equilibrate on the time scale shown here. Let us comment on the convergence properties of the NLCE data in Fig. 3. Both for Z(t) and X(t) , we observe that the NLCE remains converged for longer times if the value of g is smaller. Specifically, we find that the series breaks down at tJ ≈ 4 for g = 2, at tJ ≈ 8 for g = 1, while no breakdown can be seen for g = 0.5. Comparing these NLCE data to direct simulations of ladders with periodic boundary conditions and L x = 12, a good agreement is found on short to intermediate time scales (or even longer for g = 0.5). In particular, the simulation for the finite ladder turns out to be advantageous for a strong transverse field g = 2, since it captures the stationary value of Z(t) and X(t) for a longer time than the NLCE. Similar to our previous results for chains, however, it becomes clear from Fig. 3 (a) that the usage of NLCEs is in turn beneficial for g = 1, where finite-size effect appear to be stronger and the NLCE captures the monotonous decay of Z(t) up to longer times compared to the finite-system data.
To proceed, Fig. 4 shows results for quantum quenches starting from the initial state |ψ(0) = |→ , with data for two-leg ladders in Fig. 4 (a) and data for three-leg ladders in Fig. 4 (b). Since Z(t) = 0 due to the spin-flip symmetry of H, we only have to consider X(t) in this case. We find that X(t) generally behaves very similar for the two different ladder geometries. Specifically, X(t) rapidly decays towards an (approximately constant) stationary value which is naturally higher for a higher value of g. Note however, that for L y = 2 and g = 0.5, as well as for L y = 3 and g = 2, X(t) still exhibits some residual fluctuations, i.e., perfect equilibration is absent. Concerning the convergence properties of the NLCE, we find that, analogously to the previous case of |ψ(0) = |↑ , the NLCE remains converged significantly longer for g = 0.5 compared to g = 2. Especially the early breakdown of convergence for L y = 3 and g = 2 in Fig. 4  Numerical data obtained by NLCE for different expansion orders C are compared to direct simulations of finite ladders with PBC.
As a side remark to conclude the study of ladder geometries, let us note that Ref. [97] has recently discussed the possibility of quantum scars in transverse-field Ising ladders. Specifically, Ref. [97] has considered small values of g and "density-wave" initial states of the form |ψ(0) ∼ |↑↓↑↓ · · · . These initial states were found to exhibit a large overlap with rare, weakly entangled eigenstates, leading to quasi-periodic revivals in the dynamics. As detailed in Appendix C, the fully polarized states |↑ and |→ studied in the present paper, in contrast, do not exhibit such a significant overlap with the weakly entangled eigenstates. These special eigenstates therefore do not play a distinguished role for the quench dynamics presented in Figs. 3 and 4.

Two-dimensional square lattice
We now come to the last part of this paper, i.e., the quantum quench dynamics in the two-dimensional transverse-field Ising model. Note that dynamical properties of this model [50,52,56,58,98,99], as well as the emergence of thermalization [24,25,61], have been studied before by a variety of approaches. By comparing our results to existing data from the literature, let us demonstrate in this section that numerical linked cluster expansions based on rectangular clusters only, combined with an efficient forward propagation of pure states, provide a competitive alternative to other state-of-the-art numerical approaches.
As a first step, it is instructive to compare our results to earlier NLCE data from Ref. [72]. This comparison is shown in Figs. 5 (a) and (b), where the dynamics of the transverse magnetization X(t) is studied for quenches from |→ and |↑ with g = 1. (Recall that g c ≈ 3.044 for the two-dimensional lattice.) Importantly, Ref. [72] has considered all topologically distinct cluster geometries in the expansion and has used ED to evaluate the respective weights. Due to the computational bottlenecks of NLCEs discussed in Sec. 3.1, Ref. [72] was consequently limited to rather small clusters with up to 10 lattice sites. In Fig. 5, we find that our NLCE with solely rectangular clusters nicely reproduces the data from Ref. [72]. In particular, while the results of Ref. [72] are converged for times tJ < 1, the rectangular NLCE up to expansion order C = 28 (i.e. the largest clusters are of size 7 × 4, 14 × 2, 28 × 1) yields converged results on time scales which are approximately twice as long. This demonstration, that a NLCE restricted to rectangular cluster geometries can be better than a NLCE comprising all (pos- sibly nonrectangular) clusters, is an important result of the present paper. In this context, let us add that the inclusion of rectangles with different length ratios appears to be crucial to achieve a good convergence. For instance, we have checked that an expansion using solely square-shaped clusters (1 × 1, 2 × 2, . . . , 5 × 5) performs very poorly instead (not shown here). Moreover, let us note that the short-time dynamics in Figs. 5 (a) and (b) can be accessed also by the direct simulation of a 5 × 5 lattice with PBC. Next, let us study quenches starting from the state |ψ(0) = |↑ such that Z(0) = 1 and X(0) = 0, and consider a strong transverse field g = 2.63g c ≈ 8, i.e., a quench across the quantum critical point. Again, we consider clusters with up to 28 lattice sites in the NLCE. In Fig. 6 (a), we find that Z(t) displays pronounced oscillations with an amplitude that is weakly damped over time. Correspondingly, the transverse magnetization X(t) in Fig. 6 (b) exhibits damped oscillations as well (with a frequency that is twice as large). It is instructive to compare these NLCE data for the thermodynamic limit to a simulation of a 5 × 5 lattice with PBC. Specifically, one observes that for such a finite system and times tJ 1, the oscillations of Z(t) and X(t) die away rather quickly. This is in contrast to the NLCE results for L → ∞ which capture the persistent oscillations on a longer time scale. In addition, we compare our NLCE results for Z(t) in Fig. 6 (a) to recent data digitized from Ref. [56], which are computed by an artificial neural-network (ANN) approach for a 8 × 8 lattice. While the NLCE and ANN data agree nicely with each other for times tJ < 1, the NLCE remains converged also on longer time scales. In particular, the ANN data from Ref. [56] up to times tJ 1 can be reproduced even by the smaller 5 × 5 lattice. Thus, for the parameter regime considered in Fig. 6, it appears that the NLCE can be better than the direct simulation of finite systems with PBC as well as the ANN approach from Ref. [56]. This is another important result of the present paper.
Finally, we also consider quenches starting from the state |ψ(0) = |→ . The values of the transverse field are chosen as g = 0.1g c , 1g c , 2g c , which again allows us to compare to ANN data from Ref. [56], as well as to data from Ref. [52] based on infinite projected entangled pair states (iPEPS). For all values of g shown in Figs. 7 (a)-(c), we find a convincing agreement between the data from Refs. [52,56] and our NLCE results up to expansion order C = 28, with convergence times that are rather similar for all three methods. In order to put the convergence times into perspective, it is again helpful to compare the NLCE data to a simulation of a finite 5 × 5 lattice with PBC. While finite-size effects appear to be less important for g = 0.1g c and g = 2g c , we observe pronounced finite-size effects for g = g c already at short times tJ ≈ 0.5 due to, e.g., the divergence of the relevant length scales at the quantum critical point. Importantly, the NLCE results for g = g c in Fig. 7 (b) remain converged up to times tJ ≈ 1.5. One explanation for the advantage of NLCEs at the quantum critical point might be given by the fact that the expansion involves a variety of clusters with different ratios of width and height such that one can capture the dynamics on longer time and length scales. This is another central result of this paper.

Conclusion
To summarize, we have studied the nonequilibrium dynamics of the transverse and the longitudinal magnetization resulting from quantum quenches with fully polarized initial states in the transverse-field Ising model defined on different lattice geometries. To this end, we have relied on an efficient combination of numerical linked cluster expansions and a forward propagation of pure states via Chebyshev polynomials. Depending on the geometry and the parameter regime under consideration, the quench dynamics has been found to display a variety of different behaviors ranging from quick equilibration, over slower monotonous relaxation, to persistent (weakly damped) oscillations. As a main result, we have demonstrated that NLCEs comprising solely rectangular clusters provide a promising approach to study the dynamics of two-dimensional quantum many-body systems directly in the thermodynamic limit. While the organization of the NLCE becomes straightforward due to the simple cluster geometry, the memory efficient pure-state propagation made it possible to include clusters with up to 28 lattice sites. Especially, for quenches to the quantum critical point, where finite-size effects are typically strong, we have shown that NLCEs can yield converged results on time scales which compare favorably to direct simulations of finite systems with periodic boundary conditions (also in the case of chains or ladders). By Data obtained by NLCE for expansion orders C = 24, 27, 28 are compared to the simulation of a 5 × 5 lattice with PBC. Additionally, we show iPEPS data digitized from Ref. [52] and ANN data for a 8 × 8 lattice digitized from Ref. [56].
comparing to existing data from the literature, we have demonstrated that the reachable time scales are also competitive to other state-of-the-art numerical methods. While NLCEs with rectangular clusters have been used before to obtain thermodynamic quantities [100] or entanglement entropies [68], the present paper unveils that such NLCEs also provide a powerful tool to study the real-time dynamics of quantum many-body systems (although truly long times are still out of reach). A natural direction of future research is to further explore the capabilities and the convergence properties of NLCEs regarding the real-time dynamics of quantum many-body systems. In this context, it might be promising to consider other building blocks for the expansion such as, e.g., clusters that consist of multiple corner-sharing 2 × 2 squares [66]. Moreover, it will be interesting to study other two-dimensional lattice geometries such as triangular or Kagome lattices with nonrectangular cluster shapes. Finally, the usage of supercomputing will be helpful to include higher expansion orders and to extend the convergence of the NLCE to even longer time scales [76].
Note added: After this paper was submitted, we became aware of the related work [101] which appeared in the same arXiv posting as our manuscript. While Ref. [101] also presents NLCE calculations for the dynamics of two-dimensional systems using an expansion in rectangles, its focus is on the application of NLCEs to disordered systems or inhomogeneous initial states. In addition, while Ref. [101] employs ED to evaluate the contributions of the clusters, the present paper highlights the usefulness of efficient pure-state propagation methods to reach expansion orders beyond the range of ED and to extend the convergence times of the NLCE.

A Exact solution for the integrable chain
In the case of a chain geometry, the transverse-field Ising model (1) is a paradigmatic example of an integrable model and can be diagonalized by means of subsequent Jordan-Wigner, Fourier, and Bogolioubov transforms [60], Since quantum quenches in the transverse-field Ising chain have been studied extensively before, and since the focus of this paper is on the numerical analysis of nonintegrable geometries, we here refrain from providing more details and refer to the large body of existing literature instead [89][90][91][92][93][94][95][96]. Given the notation of H in Eqs. (1) and (15), as well as an initial state |ψ(0) which is chosen as the groundstate of H for some transverse field g , the dynamics of the transverse magnetization X(t) for a quench g → g is then given by [86][87][88]90], where we have used the abbreviations and E k and k are defined like their unprimed counterparts, but with g → g . In order to obtain the results shown in Fig. 2 of the main text, we have numerically evaluated the integral in Eq. (16) either for g = 0 (|ψ(0) = |↑ ) or for g → ∞ (|ψ(0) = |→ ).

B Accuracy of the pure-state propagation
While we have already demonstrated that our numerical results agree very well with existing data, let us nevertheless discuss the accuracy of the Chebyshev-polynomial expansion which is used to evaluate the time evolution of the pure states |↑ and |→ . To this end, Fig. 8, shows the dynamics of the transverse magnetization X(t) for a cluster of size L x × L y = 7 × 3 (with OBC), initial state |ψ(0) = |→ , and transverse field g = g c ≈ 3.044. First, in Fig. 8 (a), we set the discrete time step to δtJ = 0.02 and depict results for different expansion orders M = 5, 10, 15, 20 (curves). On the one hand, for small M = 5, we observe clearly unphysical results (e.g. X(t) > 1), which can also be explained by the fact that the norm ψ(t)|ψ(t) (symbols) is not conserved over time for this choice of M . On the other hand, for M = 10, 15, 20, all curves for X(t) are perfectly on top of each other, i.e., convergence with respect to M has been reached, and ψ(t)|ψ(t) = 1.

C Eigenstate entanglement and spectral decomposition of initial states
Let us discuss some properties of the fully polarized initial states |↑ and |→ . To this end, we first study the entanglement (von Neumann) entropy S |n of the eigenstates |n of H, where ρ A is the reduced density matrix on a subsystem A, obtained by tracing over the degrees of freedom in the complement B. In Figs. 9 (a) and (b), S |n is shown for a chain and a two-leg ladder respectively, numerically obtained by exact diagonalization for L = 14 sites, transverse field g = 0.5, and periodic boundary conditions. In both cases, we have chosen A as one half of the system, i.e., the first 7 lattice sites in case of the chain, or the first three rungs and one site of the fourth rung in case of the ladder. On the one hand, for the integrable chain geometry in Fig. 9 (a), we find that S |n is comparatively small at the edges (consistent with the area-law entanglement scaling of groundstates [102]), while weakly and strongly entangled states coexist in the bulk of the spectrum (see also Refs. [103,104]). On the other hand, for the two-leg ladder in Fig. 9 (b), the fluctuations of S |n in the center of the spectrum are clearly smaller, i.e., the eigenstates are typically stronger entangled. This behavior of S |n can be interpreted as an indication of the transition from integrability to nonintegrability [104], by going from chains to ladders. In addition, we can identify a small number of eigenstates |n with energy close to E = 0 in Fig. 9 (b), which exhibit a distinctly lower value of S |n . This appears to be consistent with the recent proposal of quantum scars in transverse-field Ising ladders in Ref. [97]. Next, it is useful to study S |n in combination with the overlap between the initial states |ψ(0) = |↑ , |→ and the eigenstates |n , where E n is the eigenvalue of H belonging to |n . As shown in Figs. 9 (c) and (d), this spectral distribution is narrow and peaked at the groundstate in the case of |↑ , while P |ψ is much broader for |→ , both for the chain and the ladder. Thus, a quench to g = 0.5 with |ψ(0) = |↑ , results in a dynamics which is strongly dominated by the groundstate with a significantly smaller admixture of excited states. Note that exactly for such a situation, i.e., a quantum many-body system with one macroscopically populated eigenstate, an analytical prediction for the temporal relaxation has been recently obtained in Ref. [105]. While this is beyond the scope of the present manuscript, it appears that quantum quenches in transversefield Ising chains or ladders can be promising candidates to test such predictions. Finally, as already pointed out in the main text, we note that the fully polarized initial states |↑ and |→ do not exhibit a distinguished overlap with the rare, weakly entangled eigenstates discussed in Fig. 9 (b). These potential quantum scars therefore do not play an essential role for the resulting quench dynamics.