Hybrid infinite time-evolving block decimation algorithm for long-range multi-dimensional quantum many-body systems

In recent years, the infinite time-evolution block decimation (iTEBD) method has been demonstrated to be one of the most efficient and powerful numerical schemes for time-evolution in one-dimensional quantum many-body systems. However, a major shortcoming of the method, along with other state-of-the-art algorithms for many-body dynamics, has been their restriction to one spatial dimension. We present an algorithm based on a \textit{hybrid} extension of iTEBD where finite blocks of a chain are first locally time-evolved before an iTEBD-like method combines these processes globally. This in turn permits simulating the dynamics of many-body systems in the thermodynamic limit in $d\geq1$ dimensions including in the presence of long-range interactions. Our work paves the way for simulating the dynamics of many-body phenomena that occur exclusively in higher dimensions, and whose numerical treatments have hitherto been limited to exact diagonalization of small systems, which fundamentally limits a proper investigation of dynamical criticality. We expect the algorithm presented here to be of significant importance to validating and guiding investigations in state-of-the-art ion-trap and ultracold-atom experiments.


I. INTRODUCTION
Over the past century, a lot of research in condensed matter physics has focused on exploring, understanding, and engineering exotic phenomena emerging from the interactions of many particles whose individual behavior falls short of the richness contained in that of their collective gestalt. 1 Straightforward examples of the emergent behavior of interacting microscopic particles are the distinct phases of matter due to spontaneous breaking of global symmetries, topological phase transitions and topological order, laminar flow, turbulence, and countless other phenomena the appearance of which necessarily requires the collective interactions of the microscopic constituents. 2 In recent years, the study of dynamical criticality in quantum many-body systems has received huge interest. 3,4 In particular, two concepts of dynamical phase transitions have been at the forefront of these investigations, where both are mostly concerned with dynamics in the wake of a quench. The first relies on the dynamical behavior of a local order parameter. Criticality can then be determined by the infinite-time value of this order parameter, such that when it takes a finite value this means that the system has settled into a (global) symmetry-broken steady state, whereas a zero order parameter indicates that the system has settled into a symmetric steady state. 5,6 However, criticality can also be determined by the intermediate or prethermal dynamics of the order parameter, rather than its value in the steady state. This is particularly intuitive in the case of models with no finite-temperature phase transition where a quenched system always ends up in a symmetric steady state. Nevertheless, how the order parameter evolves in time can exhibit either one of two dynamical behaviors (phases): an asymptotic decay to zero where the order parameter never changes sign, or an envelope that oscillates about zero while asymptotically decaying to it. These phases have been characterized for quenches within the ordered phase and from the ordered phase to the symmetric phase, respectively, in the nearest-neighbor quantum Ising, 7,8 Bose-Hubbard, 9 and XXZ 10 chains. Recently, this notion of critical behavior in the intermediate dynamics of the order parameter has been extended to long-range interacting quantum spin chains. 11 The second concept of dynamical phase transitions builds on interpreting the overlap of the time-evolved wave function of a many-body system with its initial state as a dynamical analog of the thermal partition function, with complexified evolution time standing in for inverse temperature. 12,13 Therefore, just like critical temperatures bring about thermal phase transitions in equilibrium, critical times define so-called dynamical quantum phase transitions (DQPT) in out-of-equilibrium manybody dynamics. The theory of DQPT has recently witnessed rapid development, with studies of its properties covering various integrable and nonintegrable systems.
Nevertheless, these various phenomena require great effort to generate whether experimentally or numerically. In generic quantum many-body systems the effective Hilbert space grows exponentially with system size, rendering exact diagonalization methods simply insufficient, particularly when one is interested in critical behavior, and more so when that criticality is in the out-of-equilibrium realm. Sophisticated numerical simulations have taken on a key role in studying such systems, in particular the celebrated density matrix renormalization group (DMRG) algorithm, 14,15 and its implemen-tation using matrix product states 16 (MPS) has allowed the exploration of various equilibrium phases of matter. Additionally, its extension to simulating the time evolution of quantum many-body systems has revealed a vast range of exotic dynamical properties in strongly interacting systems. [17][18][19][20][21] At the same time, the development of high-precision experimental techniques has made possible the realization of highly tunable quantum systems in the laboratory. Such engineered systems enable the study of not only the static properties of the underlying model, but also its dynamical properties, such as for example in the wake of quench, the parameters of which can now be very well controlled in ion-trap 22 and cold-atom 23 experiments. Moreover, state-of-the-art setups today can realize quantum many-body systems with various forms and ranges of interaction profiles and in different lattice geometries. [24][25][26][27][28] Such quantum simulators 29-32 promise to provide deep insights into the physics of both microscopic and macroscopic systems, including the study of critical dynamical phenomena in generic quantum manybody models.
To probe, verify, and calibrate a quantum simulator, the development of methods to reliably predict its properties is vital. Currently, the static properties of longrange interacting one-and two-dimensional systems in the thermodynamic limit are well studied using extended DMRG techniques. 33,34 In the past few years, time evolution in long-range quantum spin chains has witnessed considerable success through the time-dependent variational principle (TDVP) in uniform MPS. 19,35,36 The dynamics of models on finite two-dimensional lattices has also been well studied using several time-evolution schemes. 20,21 However, uncovering the dynamical properties of higher-dimensional models would be the useful next step in light of the tremendous advancement in experimental setups that can now readily and reliably realize these systems. Nevertheless, numerical techniques have been restricted to just one spatial dimension when it comes to generic quantum many-body models in the thermodynamic limit. For example, numerical studies on both of the above concepts of dynamical phase transitions have been restricted to one-dimensional (1D) many-body Hamiltonians, with the exception of integrable models, such mean-field and exactly solvable free-fermionic systems, and small finite two-dimensional lattices studied in exact diagonalization where ascertaining dynamical critical behavior is necessarily inadequate so far away from the thermodynamic limit.
In this work, we introduce a new MPS algorithm that can simulate the dynamics of effectively two-dimensional (2D) systems by time-evolving mappings thereof to 1D chains with long-range fixed-length interaction profiles. This method is a combination of two techniques. The first involves local time evolution of several sites on the 1D chain using an algorithm for time-evolving a finite system. The second technique is based on iTEBD-like evolution that evolves different bulks within the 1D Hamilto-nian using Suzuki-Trotter expansion. Very importantly, this method allows simulating the time evolution of longrange 1D models without assuming translational symmetry up to the size of a unit cell. By mapping 2D systems of infinite length in just a single direction to a 1D chain containing a unit cell with the appropriate interaction profile and corresponding size, our method is able to simulate the dynamics of such 2D systems and adequately determine their critical properties.
The rest of the paper is organized as follows: In Sec. II we provide an overview of matrix product states and operators, followed in Sec. III by a description of the Suzuki-Trotter decomposition and Krylov-subspace expansion, both of which form the cornerstones of our hybrid timeevolution algorithm. In Sec. IV the main features of our method -the hybrid infinite time-evolution block decimation (h-iTEBD) algorithm -are presented along with an error analysis. We explain in Sec. V how one can map two-dimensional models to one-dimensional chains of fixed-length long-range interactions that are amenable for implementation using our method. Further benchmarking is presented in Sec. VI in the form of a comparison of prominent DQPT results calculated using our method and other well-established techniques in 1D, and we additionally present results on DQPT in the quantum Ising model on a triangular lattice. We conclude in Sec. VII.

II. MATRIX PRODUCT STATES AND MATRIX PRODUCT OPERATORS
Matrix product states (MPS) are the representation of a state on a lattice with L sites through a product of matrices. 37,38 Let the local state at the i th site be σ i , which can take any state within the local Hilbert space at site i. A component of the probability amplitude associated with a local state σ i is represented by a matrix M σi i with a bond dimension m. The probability amplitude of a global state |σ 1 , . . . , σ L can then be represented as a product of all the M σi i matrices. With this representation, a state |Ψ can be represented with L matrices with dimensions of m × m as follows: A product of two matrices is invariant upon substituting identity in the middle. Therefore, one can impose a condition such as through some local unitary transformation. A wavefunction can then be represented with A (B) matrices as is what remains after left (right) orthogonalizing all of the M matrices from the left (right). The representation in (4) can further be canonicalized using Vidal's Λ Γ notation: 39 where the original A σm m matrix corresponds to Λ σm m Γ σm m . Its infinite extension can be achieved by infinitely repeating the unit cell of length L.
An MPS |Ψ with a bond dimension m can be compressed to an MPS |Ψ with a lower bond dimension m < m. 37,40 For M i to be an optimal matrix that best approximates M i with lower dimension, it must satisfy Similarly to MPS, operators on a system can be represented in a matrix product form. 21,37,38,41,42 This form is called a matrix product operator (MPO) and is formulated as follows. In general, the component of an operator on a chain can be written as follows: If an operator can be expressed as a product of local operators, then most of the W matrices are identity. In such a case, MPO has a more efficient form. When a HamiltonianĤ consists of the sum of local terms, it can be written as a sum of parts localized to the left (Ĥ Left ) and right (Ĥ Right ) of a boundary ∂i, which resides between the i th and i + 1 st sites. Then this Hamiltonian can be decomposed aŝ Consequently, this is the equivalent of a product of the matrices of local operators. Decomposition ofĤ at ∂i in the i th iteration looks likê whereÔ ,Ô ,Ô ,Ô are local operators that recover the summation term in (11). This decomposition can be done until the unit cell size is reached. The MPO of an infinite system is then just an infinite repetition of L local matrices.

III. OPERATOR EXPONENTIATION: SUZUKI-TROTTER DECOMPOSITION AND KRYLOV SUBSPACE EXPANSION
An operator exponentiation of a sum of noncommuting operator is a nontrivial challenge when it comes to time-evolving systems with long-range interactions. One way to perform it is by approximating the exponent of the sum to a product of exponentiated summands. For a sum of operators with small norm, this is possible with a controlled error bound in Suzuki-Trotter decomposition. 17,37,39,43,44 It decomposes e −i∆t(Â+B) as This can be expanded to one higher order, without requiring any additional computational resources, as This allows exponentiation to be performed individually in the subspace of the summands with an error of the order ∆t 3 . Ideally, the operatorsÂ andB would be sums of commuting summands, i.e., [Â,B] = 0. However, this is generically not so. In the case where [Â,B] = 0, it is optimal that the value of the commutator be as small as possible.
After the decomposition, the exponentiated summands must be calculated. However, exponentiation using exact diagonalization is numerically very expensive if not Full Time Evolution by ∆t where |Ψ(t 0 ) is a Krylov vector at time t 0 . Once the approximated operator is found then exp(−i∆tM) -or, in general, any function ofM -can be calculated by the method of one's choice.
The Lanczos method 20,47-49 is one of several algorithms that one can use to proceed from this point. It tridiagonalizes an operatorM into a small operator by iteratively finding a suitable Krylov basis starting from |Ψ(t 0 ) . If there is such a transformation, thenM can be written as follows: Acting withM on |Ψ(t 0 ) while enforcing the orthornormalization condition V i |V j = δ i,j yields a recursive relation for α i , β i , and the Lanczos vectors |V i as follows: One can set a small parameter and stop the iteration when the value of β i goes below that threshold. The tridiagonalized operator with a reduced rank can then be exponentiated using the method of one's choice.
Once the subspace is spanned by the Krylov vectors, we can write where c i is the (i, 0) th component of exp(−i∆tM) in the subspace.
In general, the error from the approximation has the bound 49-51 Therefore, the error scales as O (∆t n ). This error bound is only valid if {|V i } are orthonormalized to infinite precision. However, the orthonormality of {|V i } expressed as MPS breaks down upon compression. This error can be reduced by introducing correction terms through adopting semiorthogonal 52 or full reorthonormalization methods. 48 In the case of our algorithm, such measures are unnecessary because h-iTEBD does not require Krylov subspace dimension of more than 3. This is because the error from the 2 nd -order Suzuki-Trotter decomposition scales as ∆t 3 such that the incurred accuracy from a Krylov subspace expansion larger than 3 does not contribute significantly to improving the precision. For three Krylov vectors, the error from broken orthonormalization is effectively negligible.

IV. ALGORITHM
The hybrid infinite time-evolving block decimation (h-iTEBD) algorithm presented here evolves a state in time under a Hamiltonian with translational invariance down to a length scale of L sites. Inspired by the iTEBD algorithm, h-iTEBD achieves this by breaking the Hamiltonian into two parts: an intra-unit cell part where the interactions are contained within the unit cell, and an inter-unit cell part where the interactions cross the unit cell boundaries. Each part can then be evolved locally by some local time-evolution algorithm, which in our case is Krylov subspace expansion. Those separate local evolutions of inter-and intra-unit cell parts are then combined via Suzuki-Trotter iterations, as shown in Fig. 1. Although Krylov subspace expansion is used for the local time evolution, other methods such as TDVP 19,35,53 can equally well be used, and this is one of our future directions. Fig. 1 shows the schematic diagram of the time-evolution routine using a tensor network diagram.
We are interested in time-evolving a Hamiltonian that is translationally symmetric with unit cell size of N sites. The Hamiltonian can be expressed as a sum of the com-ponentĤ intra,i , which has only terms with interactions contained within the i th unit cell, and a second com-ponentĤ inter,i , which contains all the interactions that cross the unit cell boundary between the i th and i + 1 st unit cells:Ĥ = l Ĥ intra,l +Ĥ inter,l .
Using Suzuki-Trotter decomposition, the time-evolution operatorÛ = exp −i∆tĤ can be approximated aŝ U acting on a state |Ψ can then be calculated with Krylov subspace expansion. The inter-and intra-unit cell Hamiltonians have freedoms on distributing the terms that qualify to be in either one of them. The total error can be minimized by choosing their distribution such that Ĥ inter,m ,Ĥ intra,m and Ĥ inter,m+1 ,Ĥ intra,m are smallest. For the calculations performed in this paper, we have chosen the distribution such that the overlapping terms inĤ inter,m have the same contribution as the ones inĤ intra,m . There are three contributors to the error in this method. In the earlier time steps, these are, from largest to smallest, the error from the Suzuki-Trotter decomposition, the error from the local time evolution, and the truncation error. Shown in Fig. 2 is the early-time error of the method analyzed using the Heisenberg model by following the error analysis in Zaletel et al. 21 The Néel state is used for the initial state and the fidelity F(t) = 1− Ψ(t)|Ψ(t) reference is computed using 4 th -order iTEBD as a reference exact state.
We have used h-iTEBD with Krylov subspace dimension of 3 (h-iTEBD-K3) and 7 (h-iTEBD-K7). As expected, the early-time behavior differs significantly from the regime where the error from Krylov subspace expansion is comparable to the subsequent Suzuki-Trotter step (h-iTEBD-K3), to that where the error is negligible (h-iTEBD-K7). In the former regime, the error scales as t 3 , which indicates that the most dominant error comes from the local time evolution (Krylov subspace expansion). In contrast, if the subspace dimension is sufficiently large (h-iTEBD-K7), the error scales as approximately t 2 , as expected from the 2 nd -order Suzuki-Trotter decomposition.
The implementation of h-iTEBD is available as part of the Matrix Product Toolkit, 54

V. MAPPING FROM TWO DIMENSIONS
We now illustrate one of the major capabilities of our algorithm, which is the simulation of the dynamics of nonintegrable 2D quantum many-body models upon an appropriate mapping to a 1D chain with long-range fixedlength interactions that is translationally invariant down to a unit cell. We shall now describe this mapping.
For simplicity and without loss of generality, let us consider a square lattice of size N sites × ∞ sites. If the lattice is translationally symmetric over the direction in which the length is infinite (horizontal according to Fig. 3), then this lattice can be mapped onto an infinite one-dimensional chain with a unit cell of N sites.
As shown in Fig. 3, the mapping can be done by numbering the sites along the direction in which the length is finite (vertical). Starting from site 0 at the top, the site number increases going downwards. Upon reaching the bottom, we continue on to the top of the next vertical slice to the right, and so on and so forth. Obviously, sites i and i + N have the same interaction profile. The resulting one-dimensional chain is translationally symmetric down to a unit cell with N sites, with a nontrivial long-range fixed-length interaction profile. Consequently, h-iTEBD can be used to simulate the time evolution on this type of lattice.

VI. DYNAMICAL QUANTUM PHASE TRANSITIONS
We also test our algorithm by calculating the Loschmidt-echo return rate 12 where |Ψ 0 is an initial state quenched by a Hamiltonian H at t = 0. This return rate is a dynamical analog of the thermal free energy in which complexified time it stands for inverse temperature. Much the same way as nonanalyticities in the thermal free energy indicate critical temperatures at which thermal phase transitions occur, nonanalyticities in (22) indicate critical times at which dynamical quantum phase transitions (DQPT) occur. 12,13 The theory of DQPT has advanced significantly in the last few years, and experimental observations of the phenomenon have been possible in ion-trap 22 and ultracoldatom 23 setups. Initially, it was demonstrated in the seminal work of Ref. 12 in the case of the nearest-neighbor transversefield Ising chain that if a quench in the transverse-field strength crosses the equilibrium quantum critical point, cusps will appear in the return rate. However, quenching across the equilibrium quantum critical point was then shown to be in general neither a necessary nor a sufficient condition to obtain nonanalytic behavior in the return rate. 56,57 Indeed, the phenomenon of DQPT proved to be significantly richer with anomalous cusps appearing for arbitrarily small quenches within the ordered phase [58][59][60][61][62][63] when local excitations formed the lowest-lying quasiparticles in the spectrum of the quench Hamiltonian. 64,65 The study of DQPT in one spatial dimension has resulted in a very good understanding of the phenomenon in both integrable and nonintegrable systems. In higher dimensions the investigations have been mostly restricted to integrable models. 66,67 However, we have recently used h-iTEBD to study DQPT in the nonintegrable twodimensional transverse-field Ising model on a square lattice by mapping it onto a chain as explaied in Sec. V. 63 In addition to presenting DQPT results on this model on a triangular lattice, we also include in this Section benchmarking results of h-iTEBD in relation to classic examples from the field.
First, we present the benchmarking results in Fig. 4, which are direct comparisons of the Loschmidt-echo return rate obtained from h-iTEBD-K7 and TDVP for the nearest-neighbor XXZ chain in Fig. 4(a), and the axial next-nearest-neighbor Ising chain (ANNNI) in Fig. 4(b). The Néel state is quenched with the XXZ Hamiltonian with J > 0, which complement those of the same model on the square lattice in Ref. 63. In both of these models, the system is implemented on an infinite cylinder. The Hamiltonian is constrained to have translation invariance, insomuch that each point on the lattice interacts with its six nearest neighbors. To use h-iTEBD, the model is projected onto an infinite one-dimensional lattice, as described in Sec. V, that is translationally invariant down to a unit cell of size N sites, where the 2D lattice is of size N sites × ∞ sites. The projection breaks the two-dimensional translational invariance of the original Hamiltonian, but this does not remove the two-dimensional characteristics of the model. For example, the magnitude of the critical transverse-field strength where the equilibrium quantum phase transition occurs is h e c ≈ 4.53J (see Appendix A) as calculated iDMRG. It is much larger than the equilibrium critical point of the one-dimensional model (1J) due to the higher connectivity which makes the system robust against fluctuations. Thus the projected quasi-twodimensional model has equilibrium characteristics beyond its one-dimensional counterpart.
Furthermore, as shown in Fig. 5, the magnetization and the Loschmidt-echo return rate exhibit behavior that is fundamentally different from that in one-dimensional short-range models. In particular, the top panel of Fig. 5 shows the dynamical critical point for quenches from h i = 0 to be h d c ≈ 2.25J, which separates between magnetization profiles that do not cross zero (for quenches to h f < h d c ) and those that do (quenches to h f > h d c ). We see that for small quenches below h d c the magnetization does not show exponential decay but rather it relaxes to a finite nonzero value, characteristic of models with a finite-temperature phase transition. This can only be possible in spatial dimension d > 1 for models with nearest-neighbor interactions and Z 2 symmetry. Furthermore, we calculate the return rate in the wake of the quench from h i = 0 to h f = 2.1J < h d c = 2.25J, shown in the bottom panel of Fig. 5. The return rate shows clear anomalous 58 behavior in that the corresponding magnetization profile makes no zero crossings yet the return rate exhibits five nonanalytic cusps in the same time interval. 64 As discussed in detail in Refs. 63 and 64, these anomalous cusps are due to the lowest-lying excitations being local in nature. In the nearest-neighbor quantum Ising chain, this is never the case as the lowest excitations are freely propagating domain walls. As first shown in Ref. 64 and proven analytically in Ref. 65, the presence of local excitations as the lowest-lying quasiparticles of the quench Hamiltonian is a necessary condition for anomalous cusps to arise. Hence, this clearly establishes that we are effectively working in 2D based on the equilirium and dynamical characteristics of this model.

VII. CONCLUSION
We have introduced an MPS algorithm for timeevolving systems with generic fixed-length long-range interactions. The algorithm is a hybrid of a local timeevolution scheme and a global one. An error analysis and a benchmarking of the algorithm has been carried out using Krylov-subspace expansion as the local timeevolution scheme and Suzuki-Trotter decomposition as the global one.
The analysis on scaling of the error with respect to the quasi-exact evolution result by 4 th -order iTEBD shows behavior that agrees with the theoretical estimate. The benchmarking of the algorithm with respect to TDVP on integrable and nonintegrable models also shows excellent agreement. It is shown that the algorithm is also applicable to short-ranged systems, as expected.
To go beyond local models, we applied the method to study DQPT in a quantum Ising triangular lattice placed on the two-dimensional surface of an infinite-length cylinder. The mapping of the lattice sites onto an infinite chain, gives rise to a Hamiltonian with nontrivial longrange interactions with translational symmetry down to a unit cell of size N sites.
The simulated dynamics of the quantum Ising triangular model reveal, in both the order parameter and Loschmidt return rate, properties that are fundamentally different to dynamical behavior in short-range onedimensional systems. In particular, we see that the magnetization in the wake of a small quench within the ferromagnetic phase decays to a finite nonzero value indicative of the presence of a finite-temperature phase transition. Indeed, such behavior cannot occur in short-range chains as the lower critical dimension for Z 2 symmetry is d = 1. Moreover, the return rate displays anomalous cusps, defined as cusps that occur without the order parameter ever making zero crossings, which is in stark contrast to the nearest-neighbor quantum Ising chain described in the seminal work of Ref. 12 on DQPT, where cusps are one-to-one connected to zeros of the magnetization. The fact that the algorithm can calculate these critical phenomena in out-of-equilibrium illustrates its applicability and effectiveness in simulating the dynamics of models in dimensions d > 1. It is worth mentioning that the results we have presented here on the quantum Ising triangular lattice are in full qualitative agreement with those of the quantum Ising square lattice of Ref. 63.
The future prospects lead in two directions. Firstly, considering the use of different local time-evolution algorithms, such as TDVP, may yield much more precise results. This may actually resolve the necessity for truncating long-range interactions beyond the unit-cell length, and may allow the time evolution with truly long-range Hamiltonians. Secondly, this technique can be readily applied to many of the paradigmatic two-dimensional lattice models of current interest in state-of-the-art experiments in ion-trap and cold-atom laboratories. Determining the dynamical critical point, or at least calculating the short-time dynamics of such models can be fruitful on benchmarking and testing experimentally realizable strongly correlated systems. The equilibrium quantum phase transition is at h ≈ 4.53J, vastly different from its 1D counterpart due to much higher connectivity on the triangular lattice. Fig. 6 is the ground state of the quasi-two-dimensional ferromagnetic Ising model on a triangular lattice given by the Hamiltonian (25) with 6 sites per circumference. The phase diagram for its antiferromagnetic counterpart has been calculated in infinite DMRG (iDMRG) using a mapping similar to that in Sec. V in Ref. 69. An MPS ansatz of the ground state is calculated using iDMRG. The critical point of the quantum phase transition for lattice is predicted to happen at h c ≈ 4.53J. This critical value is larger than that of the model on a square lattice and much larger than its one-dimensional counterpart. This is due to the increased connectivity in the triangular lattice, which would then require a large critical field for fluctuations to be strong enough to break order in the system.