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

TL;DR

This paper introduces a refined Hamiltonian Truncation method for strongly coupled 2D quantum field theories, reducing UV cutoff dependence and improving spectral accuracy through cubic order renormalization, with potential applications in higher dimensions.

Contribution

A novel Hamiltonian Truncation variant that integrates out high energy states at cubic order, enhancing accuracy and applicability in strongly coupled QFTs.

Findings

Smaller UV cutoff dependence in spectra

More accurate spectral results in 2D scalar theory

Potential extension to higher-dimensional QFTs

Abstract

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.

Tables

Table 1. Table 1: m ph subscript 𝑚 ph m_{\text{ph}} and Λ Λ \Lambda extracted with NLO-HT.

$g$	$m_{\text{ph}}$	$\Lambda$
0.2	0.979733(5)	$-0.0018166\left(5\right)$
1	0.7494(2)	$-0.03941\left(2\right)$
2	0.345(2)	$-0.1581\left(1\right)$

Equations

$${\Delta }_V<d/2.$$

$$H=H_0+V,H_0|i⟩=E_i|i⟩,V_{ij}=⟨i|V|j⟩.$$

$$H.c=\mathcal{E}c,c\in \mathcal{H},$$

$$H_{ll}.c_l+V_{lh}.c_h$$

$$V_{hl}.c_l+H_{hh}.c_h$$

$$H_{ll}.c_l={\mathcal{E}}_{\mathrm{raw}}{c}_l.$$

$$c_h=(\mathcal{E}-H_{hh}{)}^{-1}.V_{hl}.c_l.$$

$$H_{\mathrm{eff}}.c_l=\mathcal{E}{c}_l,$$

$$H_{\mathrm{eff}}$$

$$\Delta H(\mathcal{E})$$

$$\Delta H(\mathcal{E})=\sum _{n=2}^{\infty }\Delta H_n(\mathcal{E}),$$

$$\Delta H_n=V_{lh}\frac{1}{\mathcal{E}-H_{0hh}}{\left(V_{hh}\frac{1}{\mathcal{E}-H_{0hh}}\right)}^{n-2}{V}_{hl},$$

$$\Delta H\approx \Delta H_2^{\text{local}},$$

$$(\mathcal{E}-H_{hh}{)}^{-1}.V_{hl}|i⟩,|i⟩\in {\mathcal{H}}_l.$$

$$|{\Psi }_i⟩\equiv ({\mathcal{E}}_{\ast }-H_0{)}^{-1}.V_{hl}|i⟩,$$

$$H_{ll}.c_l+H_{lt}.c_t$$

$$H_{tl}.c_l+H_{tt}.c_t$$

$$(H_{lt}{)}_{ij}$$

$$(H_{tt}{)}_{ij}$$

$$[H_{ll}+\Delta \tilde{H}].c_l={\mathcal{E}}_{\mathrm{RR}}{c}_l,$$

$$\Delta H_2\frac{1}{\Delta H_2-\Delta H_3+({\mathcal{E}}_{\mathrm{RR}}-{\mathcal{E}}_{\ast })G}\Delta H_2.$$

$$\Delta \tilde{H}=\Delta H_2\frac{1}{\Delta H_2-\Delta H_3}\Delta H_2.$$

$$\Delta \tilde{H}=\Delta H_2+\Delta H_3+\dots .$$

$$S=\frac{1}{2}\int d\tau dx[:(\partial \varphi {)}^2+m^2{\varphi }^2:+g:{\varphi }^4:].$$

$$\varphi (x,\tau =0)=\sum _k\frac{1}{\sqrt{2L{\omega }_k}}(a_k{e}^{ikx}+a_k^{†}{e}^{-ikx}),$$

$$H=H_0+g{V}_4+\dots ,H_0=\sum _k{\omega }_k{a}_k^{†}{a}_k,$$

$$\displaystyle V_{4}=L\sum_{\sum k_{i}=0}\frac{1}{\prod\sqrt{2L\omega_{i}}}\Big{[}a_{k_{1}}a_{k_{2}}a_{k_{3}}a_{k_{4}}+\ldots\Big{]}\,.$$

$$\Delta H_2=\Delta H_2^{<}+\Delta H_2^{>}.$$

$$\Delta H_2^{>}\approx \sum _i{\kappa }_i(E_L){\int }_0^{L}dx{\mathcal{O}}_i(x).$$

$${\mathcal{E}}_1(L)-{\mathcal{E}}_0(L)=m_{\mathrm{ph}}+\Delta m(L)+O(e^{-\sigma m_{\mathrm{ph}}L}),$$

$$\Delta m(L)=-\frac{1}{8\pi m_{\mathrm{ph}}}\int d\theta e^{-m_{\mathrm{ph}}L\cosh \theta }F(\theta +i\pi /2),$$

$$F(\theta )=-4i{m}_{\mathrm{ph}}^2\sinh (\theta )\left(S(\theta )-1\right),$$

$$\Delta m(L)/m_{\mathrm{ph}}\approx b{K}_1(m_{\mathrm{ph}}L)+\frac{c}{(m_{\mathrm{ph}}L{)}^{3/2}}{e}^{-L{m}_{\mathrm{ph}}}.$$

$$\frac{{\mathcal{E}}_0(L)}{L}=\Lambda -\frac{m_{\mathrm{ph}}}{\pi L}{K}_1(m_{\mathrm{ph}}L)+a\sqrt{\frac{m_{\mathrm{ph}}}{16\pi L^3}}{e}^{-2m_{\mathrm{ph}}L}+\dots ,$$

Full text

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 White (1993); Perez-Garcia et al. (2007) and more general Tensor Networks Shi et al. (2006) such as MERA Vidal (2008) or PEPS Verstraete et al. (2008). 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 Yurov and Zamolodchikov (1990, 1991) and has since been applied in many contexts. See James et al. 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 111In relativistic QFTs one can also quantize on surfaces of constant light-cone coordinate. This light front quantization Brodsky et al. (1998) is also used in numerical solutions of strongly coupled QFTs via a version of HT; some recent work is Katz et al. (2014a, b); Chabysheva (2016); Burkardt et al. (2016); Katz et al. (2016); Anand et al. (2017). The structure of the unperturbed Hilbert space is different from the equal time case, which leads to important differences in the numerical procedure. All technical claims in this work will refer exclusively to the equal time quantization.. 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 Klassen and Melzer (1992), the numerical convergence of the HT depends crucially on the scaling dimension $\Delta_{V}$ of the interaction $V$. If the interaction is strongly relevant, in the RG sense, then HT converges fast, but convergence rate worsens as $\Delta_{V}$ increases. This is a limitation of the method. For interaction dimensions larger that $d/2$, naive HT actually diverges Klassen and Melzer (1992). To ensure the convergence, we will assume here that

[TABLE]

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}^{\alpha}})$, 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 $\alpha=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 Hogervorst et al. (2015); Rychkov and Vitale (2015, 2016); Elias-Miró et al. (2016) (following notably Giokas and Watts (2011); see also Feverati et al. (2008); Watts (2012); Lencses and Takacs (2014)) 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$ Giokas and Watts (2011); Lencses and Takacs (2014); Rychkov and Vitale (2015, 2016); Elias-Miró et al. (2016) and in one study in $d=2.5$ Hogervorst et al. (2015). We hope that in the future Hamiltonian Truncation will develop into an accurate numerical method, applicable also in $d\geqslant 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 Elias-Miro et al. (2017).

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}$:

[TABLE]

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}\times\mathbb{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:

[TABLE]

where ${\cal E}$ is the energy of a given state, and $c$ the corresponding eigenvector living in the Hilbert space ${\cal H}$ spanned by the eigenstates of $H_{0}$. Eq. (3) is infinite-dimensional and cannot be solved on a computer. So we split ${\cal H}$ into a finite-dimensional “low-energy” part ${\cal H}_{l}$ and a “high-energy” part ${\cal H}_{h}$. Motivated by effective field theory, a natural choice is to include into ${\cal H}_{l}$ all the states with energy below a given cutoff $E_{i}\leqslant 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:

[TABLE]

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 ${\cal H}_{l}$, ${\cal H}_{h}$. Similarly, $H_{ll}=P_{l}HP_{l}$, and so on.

The raw HT consists in throwing out all the states in ${\cal H}_{h}$ and solving the eigenvalue equation,

[TABLE]

By the min-max theorem, as the cutoff $E_{T}$ is increased, the eigenvalues ${\cal E}_{\rm raw}$ approach the exact eigenvalues ${\cal E}$ from above. As shown in Hogervorst et al. (2015), the raw HT numerical spectrum is expected to converge with polynomial rate $1/E_{T}^{\rho}$, with $\rho=d-2\Delta_{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}$:

[TABLE]

Plugging this back in (4) gives the equation

[TABLE]

where $H_{\rm eff}$ is the effective Hamiltonian operator acting on ${\cal H}_{l}$. It is given by

[TABLE]

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 222This has to be distinguished from the Numerical Renormalization Group (NRG) improvement of the HT Konik and Adamov (2007) à la Wilson’s NRG Wilson (1975). This method raises the cutoff by adding new chunks of the Hilbert space and tossing away the states which have low overlaps with the interacting eigenstates. Other ideas to extend the reach of HT include sweeping and reordering (see James et al. ). We have not used any of these interesting tricks in our work.. Integrating out the high-energy part of $c_{h}$ we correct or, as we say, renormalize $H_{ll}$ by $\Delta H$. While in general $\Delta 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 ${\cal H}_{l}$, relatively low and manageable on a computer.

One natural way to approximate $\Delta H({\cal E})$ would be via an expansion in powers of $V$:

[TABLE]

truncating it to a fixed order. This is what was done in the previous works Hogervorst et al. (2015); Rychkov and Vitale (2015, 2016), where (11) was truncated to the leading order (LO) $n=2$, and $\Delta H_{2}$ was computed in an analytic local approximation

[TABLE]

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. Elias-Miró et al. (2016) 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 $\Delta 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 333That can be intuitively understood as follows. For each $n\geqslant 2$, there will be states below the cutoff for which the matrix elements of $\Delta H_{n}$ grow as $\sim(cN)^{n}E_{T}$ in absolute value, where $N$ is the occupation number of the state and $c$ is a constant. For $N$ big enough, the expansion is therefore not convergent, as the truncated matrix element will outgrow the leading order contribution of $H_{0}$ growing as $\sim E_{T}$. For a detailed discussion of this point see Elias-Miro et al. (2017), appendix B..

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 Elias-Miró et al. (2016), 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 ${\cal H}_{l}\subset{\cal H}_{l}\oplus{\cal H}_{h}$. The high-energy Hilbert space ${\cal 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 ${\cal H}_{h}$ to a finite dimensional subspace ${\cal H}_{t}$ spanned by the vectors

[TABLE]

Of course, these states are impossible to compute exactly, so let us approximate them by setting $H_{hh}\approx H_{0hh}$, i.e.

[TABLE]

with ${\cal E}_{*}$ a parameter that will be set close to ${\cal E}$. We call these $|\Psi_{i}\rangle$’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\rangle,\,|\Psi_{i}\rangle\}$:

[TABLE]

where $G_{ij}=\langle\Psi_{i}|\Psi_{j}\rangle$ is the Gram matrix of the tail states, which are not orthonormal, and

[TABLE]

Here $\Delta H_{2}$ and $\Delta H_{3}$ are the same as above with ${\cal E}\to{\cal 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 ${\cal E}_{\text{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 ${\cal E}\leqslant{\cal E}_{\text{RR}}\leqslant{\cal E}_{\text{raw}}$ 444In this work, we introduce a tail state for each $|i\rangle\in{\cal H}_{l}$. This limits the number of states we include in the basis, as the full matrix ($\Delta H_{2}-\Delta H_{3}$) needs to be inverted over the space of tail states. On the contrary, the low-energy diagonalization of $H_{ll}+\Delta\tilde{H}$ is performed efficiently via the Lanczos method. This suggests that more efficient numerics could be achieved by reducing the number of tails states $|\Psi_{i}\rangle$; see Elias-Miro et al. (2017) for a discussion..

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:

[TABLE]

where $\Delta\widetilde{H}$ is given by

[TABLE]

In our calculations we will have ${\cal E}_{*}\approx{\cal E}_{\rm RR}$ 555In practice we fix ${\cal E}_{*}$ to the value given by the local approximation mentioned below Eq. (12). Further iterative improvements are possible, but their effect is negligible.. So we will neglect the last term in the denominator and will use

[TABLE]

Now observe that the power expansion of this expression agrees, up to third order in $V$, with (11):

[TABLE]

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 $\Delta\widetilde{H}$ implements a next-to-leading (NLO) renormalization correction. The presence of $\Delta H_{3}$ in the denominator of (23) is crucial to address the problems originating from the naive truncation of the expansion (12). 666By applying the power counting arguments mentioned in footonote 21, one can estimate $\Delta\tilde{H}\sim NE_{T}$, as opposed to $\Delta H_{2}\sim N^{2}E_{T}$, therefore taming the growth of the matrix elements. We will refer to the spectrum obtained via this method as NLO-HT.

Testing NLO-HT in the $(\phi^{4})_{2}$ theory. In the rest of the paper we will apply NLO-HT to one particular strongly coupled relativistic QFT—the $\phi^{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 $(\phi^{4})_{2}$ theory very briefly; see Rychkov and Vitale (2015); Elias-Miro et al. (2017) for details. The theory is defined by the normal-ordered Euclidean action

[TABLE]

We quantize it canonically with periodic boundary conditions, expanding the field into creation and annihilation operators:

[TABLE]

where $k=2\pi n/L\ (n\in\mathbb{Z})$, $\omega_{k}=\sqrt{m^{2}+k^{2}}$, $[a_{k},a_{k^{\prime}}]=0$ and $[a_{k},a^{\dagger}_{k^{\prime}}]=\delta_{kk^{\prime}}$. Here $x$ is the coordinate along the spacial circle of length $L$, while $\tau\in\mathbb{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:

[TABLE]

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 Rychkov and Vitale (2015) 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 $\mathbb{P}$ ($x\to-x)$, and field parity $\mathbb{Z}_{2}$ ($\phi\to-\phi$). We will focus on the invariant subspaces ${\cal H}^{\pm}$ consisting of states with $P=0$, $\mathbb{P}=+$, $\mathbb{Z}_{2}=\pm$. The states in ${\cal H}^{+}$ (resp. ${\cal H}^{-}$) contain even (resp. odd) number of free quanta. The basic problem is to find eigenstates of $H$ belonging to ${\cal H}^{\pm}$. 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 $\Delta H_{2,3}$ in (22) involve infinite sums over states in ${\cal H}_{h}$. We approximate $\Delta H_{2}$ to high accuracy by splitting it as Elias-Miró et al. (2016)

[TABLE]

Here the matrix $\Delta H_{2}^{<}$ involves a finite sum over the states in ${\cal H}_{h}$ of energies $E_{T}<E_{i}\leqslant E_{L}$ which is evaluated exactly. On the other hand, the matrix $\Delta H_{2}^{>}$ involves an infinite sum over the states with $E_{i}>E_{L}$, for which we use a local approximation Hogervorst et al. (2015); Rychkov and Vitale (2015):

[TABLE]

Here the ${\cal O}_{i}$ are a finite number of local Lorentz-invariant operators; for the $(\phi^{4})_{2}$ theory these are $\mathds{1}$, ${:\!\phi^{2}\!:}$, ${:\!\phi^{4}\!:}$. The coefficients $\kappa_{i}(E_{L})$ are known analytically. This approximation is most accurate for matrix elements $(\Delta H_{2})_{ij}$ such that $E_{i},E_{j}\ll 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 $\Delta H_{2}$ via (28) can be made arbitrarily small throughout the low-energy Hilbert space ${\cal 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 $\Delta H_{3}$ is analogous. We break down the matrix into various contributions. Some of those involve a finite sum over elements in ${\cal 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 Elias-Miro et al. (2017).

Numerical results. The basic features of the low-lying $\phi^{4}$ spectrum are as follows. The lowest eigenstate ${\cal E}_{0}$ belongs to ${\cal H}^{+}$ and is the ground state in finite volume (the interacting vacuum). The second-lowest eigenstate belongs to ${\cal H}^{-}$ and is interpreted as the one-particle excitation at zero momentum. The excitation energy over the ground state ${\cal E}_{1}-{\cal E}_{0}$ measures the physical particle mass $m_{\rm ph}$. The above is true for moderate quartic coupling $g<g_{c}\approx 2.8$, when the vacuum preserves the $\mathbb{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 $\mathbb{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 ${\cal E}_{0}$ and $m_{{\rm ph}}$ as functions of the coupling $g$. Notice that perturbation theory ceases to be accurate for $g\gtrsim 0.2$ (Rychkov and Vitale (2015), appendix B). We will only study here the $\mathbb{Z}_{2}$-invariant phase. The $\mathbb{Z}_{2}$-broken phase at $g>g_{c}$ was studied previously in Coser et al. (2014); Rychkov and Vitale (2016); Bajnok and Lájer (2016).

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 Yurov and Zamolodchikov (1991).

The first step is to compute the spectrum as a function of $E_{T}$ for fixed $g$ and $L$ and to extrapolate $E_{T}\to\infty$. Our NLO-HT calculations explored the couplings $g\leqslant 3$ and the volumes $L\leqslant 10$, while $E_{T}$ was fixed for each $L$ to have about $10^{4}$ states in ${\cal 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_{T}^{3}$. 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 $\sim 1/E_{T}^{2}$ 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}=\infty$. For this we fit the NLO-HT data points with the function $F(E_{T})=\alpha+\beta/E_{T}^{3}+\gamma/E_{T}^{4}$, with $\alpha$, $\beta$ and $\gamma$ free parameters, and use $F(\infty)=\alpha$ 777To estimate extrapolation errors, we fitted subsamples of the full data set, obtained by removing points at low $E_{T}$ in different combinations..

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_{\rm ph}\gg 1$ we expect, in a 1+1 dimensional QFT with unbroken $\mathbb{Z}_{2}$ symmetry Lüscher (1986); Klassen and Melzer (1991a):

[TABLE]

where $\sigma\geqslant\sqrt{3}$, and $S(\theta)$ is the S-matrix for $2\to 2$ scattering, with $\theta$ 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, 888As a further check of the method, the S-matrix (extracted via the volume dependence of the spectrum) could in the future be compared to the perturbative prediction for $g\ll 1$. See Coser et al. (2014), appendix B. but we will instead parametrize it by replacing $S(\theta+i\pi/2)$ with a series expansion around $\theta=0$. This is reasonable because the integral in $\Delta m(L)$ is dominated by small $\theta$. Eq. (31) then implies:

[TABLE]

The Bessel function comes from the constant term of the $S(\theta)$ expansion, while the second term comes from doing the integral via the steepest descent of the $\theta^{2}$ term (the linear term vanishes in the integral). Further corrections are suppressed by additional powers of $m_{{\rm 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}=\infty$ 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_{\rm ph}$, $b$, $c$) and works well in the whole range of $L$. We extract the value of $m_{\rm ph}$ at $L\rightarrow\infty$ from the fit, with the uncertainty determined by fitting the upper and lower ends of the error bars. We have done analogous $L\rightarrow\infty$ extrapolations for all couplings $g\leqslant 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}=\infty$ are also shown for comparison (green error bars). A few $L=\infty$ results are also reported in Table 1. 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_{\rm ph}\gg 1$ is not satisfied.

Also in Table 1, we report analogous measurements of the infinite volume vacuum energy density $\Lambda$ (the cosmological constant). The NLO-HT data for ${\cal E}_{0}(L)/L$ are extrapolated to $E_{T}=\infty$ and then are fitted with the theoretical expectation at $Lm_{{\rm ph}}\gg 1$:

[TABLE]

where $a=O(1)$. This formula is valid in any massive quantum field theory in 1+1 dimensions in absence of bound states Klassen and Melzer (1991b); Rychkov and Vitale (2015).

Coming back to Fig. 3, we see by eye that the mass gap vanishes somewhere close to $g_{c}\approx 2.8$, signaling a quantum critical point. This is in accord with previous theoretical Chang (1976) and numerical Schaich and Loinaz (2009); Rychkov and Vitale (2015); Bajnok and Lájer (2016); Milsted et al. (2013); Bosetti et al. (2015); Pelissetto and Vicari (2015) studies 999To compare with the critical coupling extractions using the light front quantization Chabysheva (2016); Burkardt et al. (2016); Anand et al. (2017) one has to perform nonperturbative mass renormalization Burkardt et al. (2016).. For a better estimate of $g_{c}$, we fit the $L=\infty$ data points in the range $g\leqslant 2.6$ with the rational function

[TABLE]

with fit parameters $r$, $g_{1}$, $g_{2}$, $g_{3}$, $g_{c}$, and $\nu$. 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 101010The central value corresponds to the smallest $\chi^{2}(g_{c})=\sum_{i=1}^{N}{(y_{i}-f(x_{i}))^{2}}/{err_{i}}^{2}$. The uncertainty interval was conservatively determined from the condition $\sqrt{\chi^{2}(g_{c})}\leqslant 3\sqrt{\chi^{2}(2.76)}$. Our determination is the best HT measurement of $g_{c}$. It is compatible with and has accuracy comparable to other available determinations Elias-Miro et al. (2017).

[TABLE]

The $\nu$ parameter in the above fit is a critical exponent. Assuming the Ising model universality class for the phase transition, we expect $\nu=(2-\Delta_{\epsilon})^{-1}=1$, using $\Delta_{\epsilon}=1$, the dimension of the most relevant non-trivial $\mathbb{Z}_{2}$-even operator of the critical Ising model. In the fit leading to (35) we fixed $\nu=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_{{\rm 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_{{\rm ph}}(g)=1-1.5g^{2}+2.86460(20)g^{3}+\ldots$ Rychkov and Vitale (2015). 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\lesssim g\lesssim 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

[TABLE]

where $\Delta_{I}$ are operator dimensions in the critical Ising model. This relation should hold at $L\gg 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 $\Delta_{\sigma}=1/8$, $\Delta_{\epsilon}=1$, $\Delta_{\partial^{2}\sigma}=2+1/8$. The error comes from extrapolating to $E_{T}=\infty$ and (the largest contribution) from varying $g$ in the range (35). We see reasonable agreement for $\sigma$ and $\epsilon$, while it looks like the agreement for $\partial^{2}\sigma$ will be reached at higher values of $L$. This figure can be compared to Fig. 6 of Rychkov and Vitale (2015) and Figs. 22, 23 of Bajnok and Lájer (2016), 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 $\phi^{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}=\infty$. 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 $\Delta_{V}>d/2$ excluded in this paper, and in particular to flows in higher dimensions $d\geqslant 3$, most of which fall into this category.