High-Precision Calculations in Strongly Coupled Quantum Field Theory with Next-to-Leading-Order Renormalized Hamiltonian Truncation

Hamiltonian Truncation (a.k.a. Truncated Spectrum Approach) is an efficient numerical technique to solve strongly coupled QFTs in d=2 spacetime dimensions. Further theoretical developments are needed to increase its accuracy and the range of applicability. With this goal in mind, here we present a new variant of Hamiltonian Truncation which exhibits smaller dependence on the UV cutoff than other existing implementations, and yields more accurate spectra. The key idea for achieving this consists in integrating out exactly a certain class of high energy states, which corresponds to performing renormalization at the cubic order in the interaction strength. We test the new method on the strongly coupled two-dimensional quartic scalar theory. Our work will also be useful for the future goal of extending Hamiltonian Truncation to higher dimensions d>= 3.

Introduction. Many interesting strongly interacting Quantum Field Theories (QFTs) are not amenable to analytical treatment. Such theories are often studied via Lattice Monte Carlo (LMC) numerical simulations, starting from the discretized Euclidean action. However, LMC has some drawbacks, for example it cannot easily compute real-time observables, it is rather computationally expensive, and it cannot directly describe renormalization group (RG) flows starting from interacting fixed points. Therefore, it is worth exploring other numerical approaches to strongly interacting QFTs. One promising alternative is provided by the Hamiltonian methods, which look for the eigenstates of the quantum Hamiltonian. These methods use various finite-dimensional approximations to the full infinite-dimensional QFT Hilbert space. Notable examples are the methods using Matrix Product States [1, 2] and more general Tensor Networks [3] such as MERA [4] or PEPS [5]. In this paper we will be concerned with another representative of this group of methods-Hamiltonian Truncation (HT), also known as the Truncated Spectrum (or Space) Approach, which is a direct generalization of the variational Rayleigh-Ritz (RR) method from quantum mechanics. This method goes back to the seminal work of Yurov and Al. Zamolodchikov [6,7] and has since been applied in many contexts. See [8] for a recent extensive review and the bibliography.
The idea of HT is simple. The QFT Hamiltonian operator H is split as H 0 + V where H 0 is an exactly solvable Hamiltonian whose eigenstates form the basis of the Hilbert space. One quantizes at surfaces of constant time and works in finite volume so that the spectrum is discrete [9]. The Hilbert space is then truncated to the low-lying eigenvectors of H 0 . The matrix of H in this truncated Hilbert space is diagonalized exactly on a computer, to find the low-energy spectrum of interacting eigenstates.
As was understood early on [10], the numerical convergence of the HT depends crucially on the scaling dimension ∆ V of the interaction V . If the interaction is strongly relevant, in the RG sense, then HT converges fast, but convergence rate worsens as ∆ V increases. This is a limitation of the method. For interaction dimensions larger that d/2, naive HT actually diverges [10]. To ensure the convergence, we will assume here that (1) Another limitation of HT, as of many variational methods in general, is that the Hilbert space grows exponentially with the cutoff. Specifically, the dimension grows as exp(CE α T ), where C > 0 is a theory-dependent constant and E T is the energy cutoff on the H 0 spectrum. The exponent typically depends on the spacetime dimension as α = 1 − 1/d, so this problem becomes more severe in higher d. These two limitations are the main reason why the HT has been so far applied mainly in d = 2.
Motivated by the need to mitigate the above limitations, the recent works [11][12][13][14] (following notably [15]; see also [16][17][18]) started developing the theory of renormalized HT, in which high-energy modes are not simply truncated away, but integrated out to produce an effective low energy Hamiltonian. As a result the convergence is improved. Renormalized HT has been applied in several strongly coupled QFT studies in d = 2 [12][13][14][15]18] and in one study in d = 2.5 [11]. We hope that in the future Hamiltonian Truncation will develop into an accurate numerical method, applicable also in d 3. Here we will take another step towards this goal by proposing a novel and still more accurate approach to renormalization. A more detailed technical account of our work will appear elsewhere [19].
Setup. Consider the Hamiltonian H of a QFT in finite volume, which we assume can be split into a solvable part H 0 , whose eigenfunction and eigenvectors are known, plus an interaction V , whose matrix elements are computable in the basis of eigenstates of H 0 : The H 0 can represent a free or an integrable Hamiltonian, or the Hamiltonian of a conformal field theory (CFT) on the cylinder S d−1 × R. We assume that its finite volume spectrum is discrete, which is the case for most H 0 's of interest. Notice that although numerical HT calculations are performed in finite volume, infinite volume observables can then be extracted via controlled extrapolation.
The spectrum of the interacting theory is found by solving the eigenvalue equation: where E is the energy of a given state, and c the corresponding eigenvector living in the Hilbert space H spanned by the eigenstates of H 0 . Eq. (3) is infinitedimensional and cannot be solved on a computer. So we split H into a finite-dimensional "low-energy" part H l and a "high-energy" part H h . Motivated by effective field theory, a natural choice is to include into H l all the states with energy below a given cutoff E i E T , which plays the role of a UV cutoff. This should provide a good approximation for the interacting eigenstates with energy well below the cutoff. Different types of cutoff are possible but will not be considered here. We then project the eigenvalue equation onto those subspaces: where c = (c l , c h ) is the low/high energy split of c, i.e. c l = P l c, c h = P h c, where P l , P h are the projectors on H l , H h . Similarly, H ll = P l HP l , and so on. The raw HT consists in throwing out all the states in H h and solving the eigenvalue equation, By the min-max theorem, as the cutoff E T is increased, the eigenvalues E raw approach the exact eigenvalues E from above. As shown in [11], the raw HT numerical spectrum is expected to converge with polynomial rate 1/E ρ T , with ρ = d − 2∆ V > 0 by our assumption (1). This polynomial convergence must compete with the exponential growth of states in the Hilbert space.
It is possible to do better than in (6). Instead of simply truncating (4), we use (5) to express the high-energy part c h of the eigenvector in terms of the low-energy part c l : Plugging this back in (4) gives the equation where H eff is the effective Hamiltonian operator acting on H l . It is given by The solutions of (8) are equivalent to the solutions of the original eigenvalue problem (3). Eqs. (9,10) are the starting point of renormalized HT [20]. Integrating out the high-energy part of c h we correct or, as we say, renormalize H ll by ∆H. While in general ∆H cannot be computed exactly, the goal is to approximate it sufficiently well so that solutions of (8) become close to the exact eigenenergies. The hope is that this can be done keeping the cutoff E T , and therefore the dimension of H l , relatively low and manageable on a computer.
One natural way to approximate ∆H(E) would be via an expansion in powers of V : truncating it to a fixed order. This is what was done in the previous works [11][12][13], where (11) was truncated to the leading order (LO) n = 2, and ∆H 2 was computed in an analytic local approximation which will be briefly reviewed after Eq. (29) below. This was shown to improve significantly the numerical convergence of the spectrum in E T . However, in Ref. [14] it was shown that, first of all, this method is not easily generalizable to higher orders and second, increasing the accuracy of the approximation of ∆H 2 alone does not necessarily improve the convergence. Furthermore, the naive expansion (12) is not convergent and there will appear unbounded matrix elements as the power of V is increased [21]. We will now introduce the main novelty of the present paper-an approach to renormalize the truncated Hamiltonian that neatly avoids the problems pointed out by [14], and leads to a more accurate spectrum than any previous approach.
NLO-HT as integrating out tails. Let us rethink Eqs. (6-10). Eq. (6) can be viewed as an instance of the RR approach, where the full Hamiltonian has been projected on the finite-dimensional subspace H l ⊂ H l ⊕ H h . The high-energy Hilbert space H h is infinite-dimensional, but Eq. (7) implies that we don't need all of it. Indeed, this equation says that one could retrieve the exact result by truncating H h to a finite dimensional subspace H t spanned by the vectors Of course, these states are impossible to compute exactly, so let us approximate them by setting H hh ≈ H 0hh , i.e.
with E * a parameter that will be set close to E. We call these |Ψ i 's tail states, as their linear combination approximate the high-energy "tail" of the eigenvectors. We next consider the eigenvalue equation for the Hamiltonian (3) projected on the space spanned by {|i , |Ψ i }: where G ij = Ψ i |Ψ j is the Gram matrix of the tail states, which are not orthonormal, and Here ∆H 2 and ∆H 3 are the same as above with E → E * . Assuming that the operators (18,19) can be evaluated to high accuracy, one can diagonalize (16,17) numerically on a computer and obtain the Rayleigh-Ritz eigenvalues E RR . By construction, these eigenvalues have variational interpretation with an ansatz enlarged with respect to the raw HT, implying via the min-max theorem that E E RR E raw [22].
Let us transform equations (16,17) further by integrating out the tail states. Substituting c t from (17) into (16) we get an equivalent equation for the RR spectrum: where ∆ H is given by In our calculations we will have E * ≈ E RR [23]. So we will neglect the last term in the denominator and will use Now observe that the power expansion of this expression agrees, up to third order in V , with (11): This key observation reveals the connection of the discussed method with the renormalization idea from the previous section. Although this was not obvious from the start, Eq. (23) means that ∆ H implements a next-toleading (NLO) renormalization correction. The presence of ∆H 3 in the denominator of (23) is crucial to address the problems originating from the naive truncation of the expansion (12). [24] We will refer to the spectrum obtained via this method as NLO-HT.
Testing NLO-HT in the (φ 4 ) 2 theory. In the rest of the paper we will apply NLO-HT to one particular strongly coupled relativistic QFT-the φ 4 theory in 1+1 dimensions. We stress however that the basic ideas of NLO-HT and of its implementation described below are general and can be used for many other theories.
We introduce here the (φ 4 ) 2 theory very briefly; see [12,19] for details. The theory is defined by the normalordered Euclidean action We quantize it canonically with periodic boundary conditions, expanding the field into creation and annihilation operators: where Here x is the coordinate along the spacial circle of length L, while τ ∈ R is the Euclidean time. From now on, we will use the units m = 1.
In terms of normal-ordered operators, the Hamiltonian is a sum of the free piece and the quartic interaction: The ellipsis in H in (26) refer to the Casimir energy and other exponentially suppressed corrections needed to correctly put the theory in finite volume. They are discussed in detail in [12] and defined in Eqs. (2.10, 2.18) of that paper. The Hamiltonian H acts in the free theory Fock space. There are three conserved quantum numbers: total momentum P , spatial parity P (x → −x), and field parity Z 2 (φ → −φ). We will focus on the invariant subspaces H ± consisting of states with P = 0, P = +, Z 2 = ±. The states in H + (resp. H − ) contain even (resp. odd) number of free quanta. The basic problem is to find eigenstates of H belonging to H ± . The two subspaces do not mix, and the diagonalization can be done separately.
Let's describe briefly how the matrices entering the NLO-HT eigenvalue equation (20) are computed in practice. The matrix elements of H ll are known in closed form and are straightforward to evaluate, taking advantage of the sparsity for efficiency. The matrices ∆H 2,3 in (22)  involve infinite sums over states in H h . We approximate ∆H 2 to high accuracy by splitting it as [14] ∆H 2 = ∆H < 2 + ∆H > 2 .
Here the matrix ∆H < 2 involves a finite sum over the states in H h of energies E T < E i E L which is evaluated exactly. On the other hand, the matrix ∆H > 2 involves an infinite sum over the states with E i > E L , for which we use a local approximation [11,12]: Here the O i are a finite number of local Lorentz-invariant operators; for the (φ 4 ) 2 theory these are 1, : φ 2 :, : φ 4 :. The coefficients κ i (E L ) are known analytically. This approximation is most accurate for matrix elements (∆H 2 ) ij such that E i , E j ≪ E L . Its validity is justified by the operator product expansion. The original local approximation in Eq. (13) was given by the same formula (29) but with E L = E T . So it was not accurate for states close to the cutoff. Instead, the error in evaluating ∆H 2 via (28) can be made arbitrarily small throughout the low-energy Hilbert space H l by raising E L above E T . In our calculations we find that E L = 3E T provides a sufficient approximation. The error can also be further reduced by including subleading (higher derivative) operators in (29).
The strategy for computing ∆H 3 is analogous. We break down the matrix into various contributions. Some of those involve a finite sum over elements in H h close to the cutoff and are computed exactly. The remaining pieces contain the contributions of the states much above the cutoff. Those are approximated by a sum of local operators, with analytically known coefficients [19].
Numerical results. The basic features of the lowlying φ 4 spectrum are as follows. The lowest eigenstate E 0 belongs to H + and is the ground state in finite volume (the interacting vacuum). The second-lowest eigenstate belongs to H − and is interpreted as the one-particle excitation at zero momentum. The excitation energy over the ground state E 1 − E 0 measures the physical particle mass m ph . The above is true for moderate quartic coupling g < g c ≈ 2.8, when the vacuum preserves the Z 2 invariance. At g = g c the particle mass goes to zero and the theory undergoes a second order phase transition to the phase of spontaneously broken Z 2 symmetry, with critical exponents given by the 2d critical Ising model.
We will now use the NLO-HT method to provide accurate non-perturbative predictions for E 0 and m ph as functions of the coupling g. Notice that perturbation theory ceases to be accurate for g 0.2 ( [12], appendix B). We will only study here the Z 2 -invariant phase. The Z 2 -broken phase at g > g c was studied previously in [13,25,26].
While here we will focus on the vacuum and the first excited state, we stress that higher excited states and other observables are both possible and interesting to study using the HT. E.g. one can extract the S-matrix from the volume dependence of the two-particle state energies [7].
The first step is to compute the spectrum as a function of E T for fixed g and L and to extrapolate E T → ∞. Our NLO-HT calculations explored the couplings g 3 and the volumes L 10, while E T was fixed for each L to have about 10 4 states in H l . For comparison, we will also report raw and local LO renormalized HT calculations, which were pushed to much higher E T , corresponding to about 10 6 states. As an indication of the needed computer resources, our most expensive NLO-HT data points (L = 10, E T = 20) required 40 CPU hours and 80 Gb RAM per coupling value.
Empirically, the NLO-HT spectrum was observed to converge with cutoff as 1/E 3 T . A representative plot, for the vacuum energy at g = 2, is in Fig. 1(left). This is much faster than the raw and the local LO renormalized HT predictions for the same observable, which show ∼ 1/E 2 T convergence, although LO renormalization reduces the prefactor significantly, Fig. 1(right). The smooth behavior of the NLO-HT data with E T allows us to extrapolate to E T = ∞. For this we fit the NLO-HT data points with the function F (E T ) = α+β/E 3 T +γ/E 4 T , with α, β and γ free parameters, and use F (∞) = α [27].
Next we discuss how the spectrum depends on L. There are precise theoretical expectations for this dependence, which allows us to perform interesting consistency checks, and helps to extrapolate the mass gap and the vacuum energy density to their infinite volume limits (for g not too close to g c ). For the mass gap at Lm ph ≫ 1 we expect, in a 1+1 dimensional QFT with unbroken Z 2 symmetry [28,29]: where σ √ 3, and S(θ) is the S-matrix for 2 → 2 scattering, with θ the rapidity difference.
We neglect the third term in the r.h.s. of (30), while we approximate the second one as follows. In this work we will not measure the S-matrix, [30] but we will instead parametrize it by replacing S(θ + iπ/2) with a series expansion around θ = 0. This is reasonable because the integral in ∆m(L) is dominated by small θ. Eq. (31) then implies: The Bessel function comes from the constant term of the S(θ) expansion, while the second term comes from doing the integral via the steepest descent of the θ 2 term (the linear term vanishes in the integral). Further corrections are suppressed by additional powers of m ph L.
In Fig. 2 the above expectations are compared to the g = 2 NLO-HT data. We include the NLO-HT data points at the highest E T we could reach for the given L (blue), and the NLO-HT data extrapolated to E T = ∞ as discussed above (red error bars). We also include the fit of the extrapolated data using Eq. (33) (green curve). The fit has three parameters (m ph , b, c) and works well in the whole range of L. We extract the value of m ph at L → ∞ from the fit, with the uncertainty determined by fitting the upper and lower ends of the error bars. We have done analogous L → ∞ extrapolations for all couplings g 2.6 in steps of 0.2. These are shown in Fig. 3 (red error  bars), where the L = 10 results extrapolated to E T = ∞ are also shown for comparison (green error bars). A few L = ∞ results are also reported in Table I. For g > 2.6, close to the critical point, the described fitting procedure  cannot be used, as the physical mass approaches zero, and the condition Lm ph ≫ 1 is not satisfied.  [12].
Also in Table I, we report analogous measurements of the infinite volume vacuum energy density Λ (the cosmological constant). The NLO-HT data for E 0 (L)/L are extrapolated to E T = ∞ and then are fitted with the theoretical expectation at Lm ph ≫ 1: where a = O(1). This formula is valid in any massive quantum field theory in 1+1 dimensions in absence of bound states [12,31]. Coming back to Fig. 3, we see by eye that the mass gap vanishes somewhere close to g c ≈ 2.8, signaling a quantum critical point. This is in accord with previous theoretical [32] and numerical [12,26,[33][34][35][36] studies [37]. For a better estimate of g c , we fit the L = ∞ data points in the range g 2.6 with the rational function with fit parameters r, g 1 , g 2 , g 3 , g c , and ν. We have f (g c ) = 0 by construction. We impose g 1 , g 2 , g 3 > 0 so that f (g) has poles on the negative real axis. The critical coupling estimate from this fit is [38] g c = 2.76 (3) .
The ν parameter in the above fit is a critical exponent. Assuming the Ising model universality class for the phase transition, we expect ν = (2 − ∆ ǫ ) −1 = 1, using ∆ ǫ = 1, the dimension of the most relevant non-trivial Z 2 -even operator of the critical Ising model. In the fit leading to (35) we fixed ν = 1. Relaxing this assumption gives the same central value with slightly larger error bars. The rationale behind introducing the poles into the ansatz f (g) is that they are supposed to approximate the branch cut at g < 0 that the analytically continued function m ph (g) is expected to have. We checked that modifying our ansatz, and in particular increasing the number of poles, does not affect appreciably the confidence interval for g c . We also checked that the g 2 and g 3 coefficients of our best fit are roughly consistent with the perturbation theory prediction m ph (g) = 1−1.5g 2 +2.86460(20)g 3 +. . . [12]. With a more complicated ansatz, we found fits perfectly agreeing with perturbation theory. The resulting g c values are nearly identical to (35). This is not surprising, since most of fit power relevant for constraining g c comes from 1 g 2, not from the region of small g where perturbation theory is accurate.
Finally, we compare the NLO-HT results to the expectations for the finite volume spectra at the critical point. CFT predicts that the energy levels at g = g c should vary with L as where ∆ I are operator dimensions in the critical Ising model. This relation should hold at L ≫ 1, where corrections due to irrelevant couplings die out. In Fig. 4 we test it for the first three energy levels above the vacuum, which should correspond to the operators with dimensions ∆ σ = 1/8, ∆ ǫ = 1, ∆ ∂ 2 σ = 2 + 1/8. The error comes from extrapolating to E T = ∞ and (the largest contribution) from varying g in the range (35). We see reasonable agreement for σ and ǫ, while it looks like the agreement for ∂ 2 σ will be reached at higher values of L. This figure can be compared to Fig. 6 of [12] and Figs. 22, 23 of [26], which show similar behavior.
Conclusions. In this work we proposed a variant of renormalized Hamiltonian Truncation called NLO-HT. Its main idea is to integrate out exactly a certain class of high-energy states, which allows for variational interpretation, and furthermore implements the renormalization corrections up to cubic order in the interaction strength.
We tested NLO-HT by computing the low-lying spectra of the strongly coupled two-dimensional φ 4 theory. Numerical spectra in finite volume were found to converge rapidly with the Hilbert space cutoff E T , faster than for other existing versions of Hamiltonian Truncation, and allowing controlled extrapolation to the continuum limit E T = ∞. The finite volume corrections were then removed using the theoretical knowledge of these effects in QFT. In this way we extracted highly accurate predictions for the vacuum energy density and the physical mass in the infinite volume limit, for a range of non-perturbative coupling constants.
In the future NLO-HT will be used to perform accurate studies in other strongly coupled RG flows in d = 2. In particular, it can be applied to flows starting from an interacting CFTs. We also believe that our ideas will be useful to extend Hamiltonian Truncation to weakly relevant interactions, with scaling dimension in the range ∆ V > d/2 excluded in this paper, and in particular to flows in higher dimensions d 3, most of which fall into this category.