Spectral gaps of two- and three-dimensional many-body quantum systems in the thermodynamic limit

We present an expression for the spectral gap, opening up new possibilities for performing and accelerating spectral calculations of quantum many-body systems. We develop and demonstrate one such possibility in the context of tensor network simulations. Our approach requires only minor modifications of the widely used simple update method and is computationally lightweight relative to other approaches. We validate it by computing spectral gaps of the 2D and 3D transverse-field Ising models and find strong agreement with previously reported perturbation theory results.


I. INTRODUCTION
Understanding the behavior of quantum many-body systems in 2D and 3D lies at the heart of contemporary physics.One of the key enigmas in this domain is the computation of the spectral gap, i.e., the energy difference between the ground state and the first excited state, as the spectral gap serves as a key parameter characterizing the intrinsic properties of various quantum phenomena, including superconductivity, superfluidity, spin liquids [1], and quantum annealing [2].
Here, we introduce a numerical method designed to determine the spectral gap in quantum systems.This approach is computationally lightweight and seamlessly integrates with widely used imaginary time propagators, enabling its broad applicability across diverse quantum systems.The computational efficiency of the method positions it as a versatile tool to unravel the essential characteristics of quantum systems.
To validate this numerical approach, we apply it to the well-known 2D transverse-field Ising model.The obtained results exhibit excellent agreement with existing analytical solutions derived through perturbation theory, establishing the reliability of the method.We then expand our calculations to encompass the 1D Haldane chain, as well as the 3D transverse-field Ising model.This latter challenge was previously considered beyond the reach of existing computational techniques, and as a consequence, estimates of the spectral gap in 3D have, until now, been confined to series expansions [3,4].We numerically benchmark the results of our method against these results, showcasing the broad applicability and computational prowess of this method in addressing intricate quantum systems.

II. THEORETICAL APPROACH A. Main theorem
We begin by stating the main analytical result.
We note that Theorem 1 is a modification of the theorem in Ref. [5] and that previously, other works have targeted the development of similar, albeit approximate, spectral gap expressions (see, e.g., Refs.[6,7]).It has also been long recognized that the energies of the low lying excited states can be extracted via the multiexponential fitting of a correlation function (see, e.g., Refs.[8, Sec.VII.C] and [9]).However, such a fit is numerically unstable.Importantly, Eq. ( 4) does not suffer from this drawback.
In the thermodynamic limit, quantum systems (including the 2D and 3D transverse-field Ising models studied below) often have a continuum spectrum in addition to discrete energies.Equation ( 4) is also applicable in such systems; E 1 is either the first excited state, if the system has at least two discrete energy levels, or the infimum of the continuous spectrum, if there is only one bound state.
A practical application of Eq. ( 4) to estimate the spectral gap requires a numerical scheme for imaginary time evolution to calculate the ground state.There are numerous computational techniques that could be considered for this purpose, e.g., tensor networks, quantum Monte Carlo [10], and even quantum computing-based methods [11][12][13][14][15][16], and exploration of each would constitute valuable future work.
Here, we specifically consider imaginary time evolution within tensor networks.Tensor networks, as detailed in [17][18][19][20], are an apt wave-function ansatz for capturing the entanglement structure of ground states.This approach is particularly effective in 2D [21] and in 1D using time-evolving block decimation [22].
Here, we apply the simple update method [21,35] for the imaginary time evolution of iPEPS and observables calculations [36,37].We calculate the spectral gap via Eq.( 4) by tracing an expectation value during the imaginary time evolution.

B. Transverse-field Ising model in 2D and 3D
The Hamiltonian has the following form: where i and j label sites on the square or cubic lattice, ⟨ij⟩ restricts summation over the nearest neighbor pairs, and σ x i and σ z i are the conventional Pauli matrices acting on the site i.We focus on the case of ferromagnetic coupling (J > 0).At zero temperature the system can be in two different phases: symmetry-broken ferromagnetic phase at small g/J and disordered paramagnetic phase at large g/J.In 2D, the phase transition occurs at g/J = 3.04438(2) [38], while in 3D, the estimates are less precise, g/J ≈ 5.29 [39].We will use the following parametrization: in the ferromagnetic phase, we set J = 1 and measure the gap in units of J.In the paramagnetic phase, we instead set g = 1.
The excitations in the model have the following structure: in the ferromagnetic phase, the vacuum is either all spins are aligned either up or down.An excitation in this phase is created by flipping the spin at a single site.In 2D, flipping one spin creates four domain walls.As a result, in the case of g = 0, the spectral gap is ∆ = 8J.The series expansion up to O(g 20 ) for the spectral gap in the ferromagnetic phase in 2D can be found in Ref. [40].
In the noninteracting case where J = 0, all spins are aligned along the positive x direction in the ground state.The elementary excitation corresponds to the spin flip on one site, with the excitation energy ∆ = 2g, or simply ∆ = 2 due to the chosen convention in the paramagnetic phase.The series expansion up to O(J 13 ) for the spectral gap in the paramagnetic phase can be found in Ref. [41] and Ref. [3] in 2D and 3D, respectively.
It is important to note that in both the ferromagnetic and paramagnetic phases, the excitation can be obtained from the ground state by applying the local operator σ y i , since this operator reverses the magnetization along both the z and x axes.More generally, the operator i σ y i is expected to exhibit a nonzero overlap between the ground state and the first excited state, regardless of the value of g/J.Furthermore, since σ y i contains the purely imaginary matrix elements and we take the real valued initial function |ϕ(0)⟩, the operator O = i σ y i satisfies the condition (3), thereby validating application of Eq. ( 4) for the spectral gap.

C. Gapped topological system: 1D Haldane chain
The 1D Haldane model can be viewed as the onedimensional antiferromagnetic Heisenberg chain with spin S = 1.The model Hamiltonian is It was conjectured by Haldane that this system is in the gapped phase [42,43].This was later confirmed by extensive numerical calculations [44][45][46][47] and also by the solution of the Affleck-Kennedy-Lieb-Tasaki (AKLT) model [48], which is in the same phase.Furthermore, the gapped phase is nontrivial-this is the symmetryprotected topological order, which is protected simultaneously by several symmetries [49].In the tensor network calculations these topological properties of the model expose themselves as degeneracies in the entanglement spectra [49], which we also observe in our calculations.
For the AKLT model the structure of excitations is approximately known, since these are determined by the spin triplets on the lattice bonds.The excitations (for the AKLT model) can be created by the operator O = i S y i S z i+1 (note that the additional higher-energy excitations are also created, but they vanish quickly in imaginary-time evolution).We employ the same O for the gap extraction in the Haldane chain, since the Haldane chain and AKLT models are continuously connected.

D. Spectral gap calculations via tensor networks
There are several different ways to calculate the gap and excitation spectra with tensor networks.The easiest approach is to add a penalizing ground state projector to the Hamiltonian and then re-run the groundstate calculation [50,51].This approach works well for 1D systems, where the DMRG optimization can be efficiently employed.Another strategy is based on the tensor-network approach for the excitation ansatz [52][53][54][55][56][57].This approach allows to calculate not only the gap and lowest excited state, but also the dispersion relation, certain two-particle bound states [58], the scattering matrix [59], and topological excitations [55,60].This approach can also be applied to 2D systems on cylinders [61], helixes [62] or in infinitely extended 2D systems with the iPEPS excitation ansatz [34,[63][64][65].There is also an approach based on the Lanzcos algorithm adaptation to the MPS [66,67].Note that the excitation energies can also appear in the projected Hamiltonian problems in the DMRG sweeps [68] and may be connected with the ground state transfer matrix spectrum [69,70].For certain class of models it was also proposed to estimate the gap by studying phase transitions upon the introduction into Hamiltonian of additional terms commuting with the Hamiltonian [71].For many-body localized systems there are tailored approaches to study not only the lowest excited states but also eigenstates in the middle of the spectra [72,73].
The approaches just reviewed are very effective for quasi-1D systems and for certain 2D systems within the variational iPEPS methodology.Unfortunately, for 2D systems the computations are rather demanding, since they require either complex summations of correlation functions with the corner transfer matrix renormalization group (CTMRG) [34,64] or automatic differentiation through CTMRG [65].These approaches are also currently very difficult to generalize to 3D and complex graph structures.Here, we show that the spectral gap can be easily extracted via the simplest algorithm of tensor network optimization, which is applicable to arbitrary lattice or network structures.We illustrate this with the spectral gap estimation for the 3D transverse Ising model.
Our algorithm for calculating the spectral gap runs as follows: We employ the infinite projected entangled pair states of the bond dimension D and optimize the randomly initialized iPEPS wave function (within a predefined periodic unit cell) with the imaginary time evolution (see Fig. 1).We discretize the imaginary time propagation of a step dτ and apply either the Trotterized gates representing exp(−Hdτ ) or a PEPO-like approximation of the exp(−Hdτ ), constructed via the W II method [76], to the iPEPS.Following the application of these gates/PEPO, the bond dimension D of the iPEPS wave function typically increases and must be truncated back to its original value.This truncation is performed using the simple update method, a cost-effective and widely used approach in iPEPS optimization, especially for gapped nontopological systems [21,35].
To determine the gap, it is necessary to find the operator averages of the commutator [H, O] of the Hamiltonian H and the observable O creating excitations.For example, we employ O = i σ y i for the Ising model, and for the Haldane model.The precise computation depends on the exact CTMRG contraction [77][78][79][80] at every step of the imaginary time evolution, rendering the simple update method ineffective due to CTMRG's significantly higher computational cost compared to the simple update's truncation procedures.Fortunately, our focus is on gapped systems, where it was recently demonstrated that the simple update method, with sufficiently large bond dimensions (D), can accurately determine observables in these systems [36].This finding led to the development of the gPEPS strategy, applying the simple update method to optimize the iPEPS wave function and to compute its observables [31,36] (see also a similar discussion in [37]).This approach has proven successful for challenging 2D problems, some 3D systems, and for simulating the IBM kicked Ising experiment classically [81,82].
Let us provide some additional comments on different computational schemes within the simple update iPEPS methods: (1) PEPO scheme: the optimization using the PEPO application, where PEPO approximates the exponent of the Hamiltonian exp(−Hdτ ), and the iPEPS is truncated using the superorthogonalisation canonical form [37,74,75].The unit cell consists of one site.
(2) MPO scheme: the optimization is performed by the consequent applications of horizontal and vertical MPO to the iPEPS wave function, where the MPOs encode the exponent of the Hamiltonian.The truncation is performed again using the superorthogonalisation.This approach is faster than the PEPO scheme.The unit cell consists again of one site.
(3) Gate scheme: the optimization is performed with the Trotterized gates.This is the most widely used scheme in practical calculations.Usually, superorthogonalization is not performed in the truncation procedure, but it can be used as an additional step between gate applications, which allows for the further stabilization of the algorithm.The unit cell consists of the two different tensors placed in the checkerboard pattern.This method is the most unstable in terms of the gap estimation.In particular, even if the commutator decay is approximately exponential in the approach, its numerical derivative at dτ ∼ 10 −3 shows sizable fluctuations, which are absent in the two previous schemes.The gap extraction in the approach requires in this case an additional data flattening.
Below, we employ the simple update (SU) to obtain the expectation value of the commutator [H, O], where O = i σ y i is the observable with nonzero overlap between the first excitation and ground state.The commutator expectation value is determined after each imaginary time step dτ .We then plot the logarithmic quantity C(τ ) = ln |⟨[H, O]⟩(τ )| [see also Eq. ( 4)] and determine the interval of τ with a linear dependence.Within this interval, we interpolate the function C(τ ) = C 0 − τ ∆, which contains the estimate of the spectral gap ∆.
Note that the proposed approach may encounter difficulties when addressing topological excitations, such as anyons in 2D systems or kinks/domain walls in 1D systems.This challenge arises because the single topological excitations generally cannot emerge from a translationally invariant tensor-network wave function through the action of a local operator.Specifically, it is not feasible to obtain the domain wall energy in the 1D Ising model inside the ferromagnetic phase.Consequently, we conclude that addressing topological excitations necessitates more complex approaches.Such methods are akin to those discussed in Ref. [55] for 1D systems.
Below, we present results obtained using three different schemes of the simple update iPEPS optimization: PEPO, MPO, and Trotter gates.Specifically, for the 2D system, we illustrate the iPEPS wave function, the calculation of observables, and the application of PEPO/MPO in Fig. 1.Regarding the application of Trotter gates within the simple update approach, we direct readers to Refs.[21,35] for further details.

III. RESULTS
As an initial iPEPS wave function |ϕ(0)⟩, we used random product states with bond dimension D = 1, noting that such initial states lead to better convergence.Figure 2(a) displays the dynamics of the commutator expectation value over imaginary time τ , plotted on a logarithmic scale.Computations using all three methods (MPO, PEPO, and Trotterized gates) are shown.It is clear that all methods result in the exponential decay of the commutator (evidenced as a linear region on the graph), with nearly identical slopes, indicative of the excitation gap.However, the stability and accuracy of these methods generally vary.For instance, with the Trotter gates, for large D or, surprisingly, small dτ we observe instabilities.Such instabilities, which can occur at arbitrary τ , lead to significant variations in the numerical derivatives.Consequently, for methods based on Trotter gates, one needs to use a linear regression to average out the numerical instabilities to extract the spectral gap, which makes it difficult to estimate the error from the time discretization.In contrast, the numerical derivatives obtained with PEPO/MPO methods are always smooth and stable [Fig.2(b)], enabling direct estimation of the gap.
According to Fig. 2(b), the PEPO and MPO methods estimate the spectral gap ∆ = 1.074 (almost insensitive to the step size dτ , with the calculations carried down to dτ = 0.002).As discussed below, the perturbation series expansion yields ∆ = 1.083 [see Fig. 3(a)].For Trotter gates, ∆ ≈ 1.073 [see Fig. 2(b)], but only for the specific dτ = 0.2; for dτ < 0.2, the results become sensitive to the step size due to significant numerical instabilities, as mentioned in Sec.II D. Below, we utilize the MPO method, since it is faster than PEPO and more reliable than the Trotter gates optimization.We have also calculated the average values with the much more accurate CTMRG approach.The results for the gap agree between the different approaches for the observable calculation, even though the average values of the observables themselves may be underestimated for the SU method of observables calculation.
Next, we benchmark this method of the gap estimation against the series expansions [41].In Fig. 3(a), we compare the gap values computed via the simple update imaginary time evolution (using MPOs) of the iPEPS wave function with series expansions in J (g = 1).In the region of J far from the critical transition point J c ≈ 0.329, the results from the series expansions and the simple update show a strong agreement.However, in the critical region, the predictions diverge, with the simple update underestimating the transition point, while the series expansions overestimate it.The iPEPS approach becomes unreliable in this region due to the simple update's reliance on mean-field environments, which are ineffective amidst long-range correlations.To achieve more accurate results in the near-critical regime, more complex tensor-network methods should be employed: either gap extraction from the full update evolution [83] (complemented with CTMRG for calculating averages) or employing some type of variational evolution [84].Additionally, we compare these results with D = 5 calculations, noting differences predominantly near the critical regime.
It is important to note that in the critical region, gap estimation using this approach faces challenges: first, the commutator decay tends to be polynomial rather than exponential over the extended periods of τ , and second, the simple update often inaccurately determines the precise position of the transition point, leading to a shifted gap prediction.In Fig. 3(b) we benchmark our numerical results in the ferromagnetic phase against the series expansions [40].Unfortunately, the predictions from the series expansions are reliable only for g < 1.5.Within this range, our results agree well with the series expansions.We also extend our gap prediction up to g = 2.5, where our es-  1] up to O(J 13 ) and from our approach to the gap estimation (with the MPO method of evolution).The dashed line specifies the phase transition point Jc = 0.329.(b) The estimated spectral gap in the ferromagnetic phase compared with the series expansion [40, Table 1] up to O(g 20 ).The series diverge at g ≈ 1.5, and at larger g we show only the iPEPS results.2] and up to O(g 20 ) [3, Table 3], respectively.The vertical lines indicate the estimates of the critical value Jc = 0.188 [39,85] and Jc = 0.194 [3].
timates begin to diverge from the more accurate predictions provided by the more advanced (and computationally demanding) variational iPEPS approach detailed in Ref. [64].
Next, we present the results for the 3D quantum transverse Ising model.We emphasize that variational calculations with iPEPS are not currently available in 3D, while MPS approaches are ineffective due to the rapid growth of the entanglement entropy in 3D systems.The corresponding 3D tensor-network calculations rely primarily on the simple update optimization, accompanied by various methods for subsequent calculation of ground-state observables [31,36,86,87].The 3D Ising model has been previously studied using tensor network methods [39,85], making the simple update strategy a natural choice for the gap estimation in these systems.In Fig. 4(a), we compare the 3D iPEPS estimates of the spectral gap with the series expansion from Ref. [3] in the paramagnetic regime.The critical parameter J c = 0.194 is suggested by the series expansions [3], with an approximate J c = 0.188 from previous iPEPS calculations [39,85].Our results for the gap value agree with these predictions, as we observe a slowdown in the exponential decay in the regime J > 0.18, where the gap estimation becomes challenging.Deep in the paramagnetic phase, our results are consistent with the series expansion predictions.However, in the critical regime, assessing the accuracy of our results is difficult, since the series expansions overestimate the gap, predicting ∆ SE > 0 at J = 0.2 (see Fig. 4), which exceeds the critical point estimated by other series expansions.In Fig. 4(b) we also show the results in the ferromagnetic phase.As in the 2D case, the series expansions diverge long before the critical point.Hence, we compare our results with the series expansions only for low g, where we find agreement between our results and the series expansions.
So far, our discussion has focused on Ising systems in both symmetry-broken and paramagnetic phases.To verify the applicability of our proposed method to different types of systems, particularly those with symmetryprotected topological order, we have extended our analysis to calculate the gap for the 1D Haldane model, as introduced in Sec.II C. Our estimated numerical value of the gap, ∆ = 0.410, is in close agreement with previously published numerical results, which reported ∆ = 0.410479(1) [44].The convergence of the numerical derivative of the commutator expectation value is illustrated in Fig. 5.It is also worth noting that, for 1D systems, the simple update method for calculating observables is exact.

IV. CONCLUSIONS
We have presented a tensor network realization of spectral gap calculations via Eq.( 4).Its implementation is straightforward, and only required a few extra lines to be added to an existing code for the simple update method.The runtime of our method scales as O(D z+1 ), where z is the lattice connectivity (i.e., z = 4, 6 for our 2D and 3D illustrations, respectively).Concretely, it takes about 10 min with MPO on a typical desktop computer to estimate one spectral gap in 3D [pointmarker corresponding to the bond dimension D = 3 in Fig. 4(b)].Our method is complementary to the more powerful but computationally very demanding variational iPEPS approach [64], which scales as O(D 10 ) − O(D 12 ); moreover, the latter method is only applicable to infinite regular 2D lattices and, unlike our approach, cannot be currently extended to 3D.
Looking ahead, there is significant potential for further improvement of the methodology presented.For instance, it can be generalized to arbitrary graphs with low connectivity [37,82] and to arbitrary unit cells without requiring any symmetries in the system.The use of modified environments, such as cluster environments or nearest-neighbor updates, may lead to more accurate calculations, at the expense of increased computational costs.These improvements would result in more precise calculations, particularly in critical regions.This work is part of a broader program aimed at using the simple update method to study gapped systems, which, in particular, was recently applied to simulate the IBM kicked Ising experiment [81,82].
Finally, we add that although the present work has focused on the spectral gap as a fundamental property in quantum and material sciences, it also plays an important role in vibration analysis, graph theory and network analysis, and data science.Inspired by the interdisciplinary success of tensor networks [88][89][90][91], adaption of our method to these other domains should be a subject of subsequent investigations.

FIG. 1 .
FIG. 1. Tensor network computations in 2D: (a) the network diagram of the iPEPS wave function, which consists of the tensors T placed on the sites of the square lattice.Each tensor has one physical index and four auxiliary virtual indices, which are geometrically contracted with the bond matrices λ.Both T and λ are in the superorthogonal canonical form [37, 74, 75].(b) The simple update approximation to the computation of the expectation value of the two-site operator O [36].(c)The PEPO optimization scheme, where we apply the PEPO operator (consisting of elementary tensors WP ) to the iPEPS wave function.The PEPO tensors approximate the operator exp(−Hdτ ) and then we truncate iPEPS with a superorthogonalisation[74,76].(d) The iPEPS optimization with consequent applications of horizontal and vertical MPOs (with individual tensors WM ).We construct MPO with the W II method[76] and also truncate iPEPS according to the superorthogonal canonical form.

= 3 FIG. 2 .
FIG. 2. (a) The decay of the expectation value C(τ ) = ln |⟨[H, O]⟩(τ )| with the imaginary time τ for the 2D transverse Ising model in the paramagnetic phase with J = 0.2, dτ = 0.2, and the iPEPS bond dimension D = 8 for PEPO/MPO and D = 3 for the Trotter gates (TG).(b) The numerical derivative of the commutator decay for the same parameters and methods, including TG scheme with CTMRG at χ = 10 for comparison.

FIG. 3 .
FIG. 3. (a)The gap in the paramagnetic phase obtained from the series expansion [41, Table1] up to O(J13 ) and from our approach to the gap estimation (with the MPO method of evolution).The dashed line specifies the phase transition point Jc = 0.329.(b) The estimated spectral gap in the ferromagnetic phase compared with the series expansion[40,   Table1] up to O(g 20 ).The series diverge at g ≈ 1.5, and at larger g we show only the iPEPS results.

FIG. 5 .
FIG. 5.The numerical derivative of the expectation value C(τ ) = ln |⟨[H, O]⟩(τ )| with the imaginary time τ for the Haldane model.The operator O = i S y i S z i+1 .