Krylov Spaces for Truncated Spectrum Methodologies

We propose herein an extension of truncated spectrum methodologies (TSMs), a non-perturbative numerical approach able to elucidate the low energy properties of quantum field theories. TSMs, in their various flavors, involve a division of a computational Hilbert space, $\mathcal{H}$, into two parts, one part, $\mathcal{H}_1$ that is `kept' for the numerical computations, and one part, $\mathcal{H}_2$, that is discarded or `truncated'. Even though $\mathcal{H}_2$ is discarded, TSMs will often try to incorporate the effects of $\mathcal{H}_2$ in some effective way. In these terms, we propose to keep the dimension of $\mathcal{H}_1$ small. We pair this choice of $\mathcal{H}_1$ with a Krylov subspace iterative approach able to take into account the effects of $\mathcal{H}_2$. This iterative approach can be taken to arbitrarily high order and so offers the ability to compute quantities to arbitrary precision. In many cases it also offers the advantage of not needing an explicit UV cutoff. To compute the matrix elements that arise in the Krylov iterations, we employ a Feynman diagrammatic representation that is then evaluated with Monte Carlo techniques. Each order of the Krylov iteration is variational and is guaranteed to improve upon the previous iteration. The first Krylov iteration is akin to the NLO approach of Elias-Mir\'o et al. \cite{Elias-Miro:2017tup}. To demonstrate this approach, we focus on the $1+1d$ dimensional $\phi^4$ model and compute the bulk energy and mass gaps in both the $\mathbb{Z}_2$-broken and unbroken sectors. We estimate the critical $\phi^4$ coupling in the broken phase to be $g_c=0.2645\pm0.002$.

Nonetheless a large number of questions remain open, and difficulties remain to solve.To give a few examples, it is generally difficult to compute scattering amplitudes and other matrix elements in a fully non-perturbative way, especially when the asymptotic states involve bound states or topological excitations.This holds even for a simple model like the 1 + 1d ϕ 4 theory.The control of such matrix elements is potentially important as baryons (and nuclei) appear as solitons in the low-energy chiral effective (e.g.Skyrme) models used in nuclear physics [30][31][32].Progress is possible with Hamiltonian methods, where the Lüscher formula provides direct access to elastic 2 → 2 scattering below the inelastic threshold, while spectral sums of form factors can be used to obtain scattering information above threshold [33].
It is also hard to analyze systems out of thermodynamic equilibrium.Circumstances in which this occurs include relativistic heavy ion collisions, the nucleosynthesis of the early universe [34], as well as experiments involving ultracold quantum gases [35,36].A simple example when the real time dynamics of a 1+1d dimensional QFT can be studied directly with Hamiltonian methods is the calculation of overlaps and expectation values after a sudden quench [37][38][39][40][41][42][43][44][45].
Finally, it is challenging to study systems at finite particle density and zero temperature.For example, the phase diagram of high-T c superconductor cuprates are often enriched by a quantum critical point at T = 0, around which the long wavelength excitations are described by the spin-fermion model: an effective 2+1D field theory involving scalar bosons interacting with non-relativistic fermions in the presence of a Fermi sea.The study of such models by path integral methods is hindered by the so-called sign problem, while conventional Hamiltonian approaches face difficulties due to the enormous dimensionality of the truncated Hilbert space needed to reach sensible accuracy.
It is anticipated [46,47] that these difficulties will be overcome with the eventual development of quantum computers.However, despite rapid progress in this direction, a nontrivial simulation of any quantum field theory on a physical quantum computer is still someways in the future.Development of new methods serves a double purpose: on the one hand, they may break new paths towards the solution of the above challenges.On the other hand, they are interesting in the pursuit of finding the optimal classical algorithms to benchmark the evolving quantum computers on simpler toy models.
Of the methods listed above, truncated spectrum methodologies (TSMs) are perhaps the least explored, particularly for QFTs in two or more spatial dimensions.TSMs are a Hamiltonian approach that is essen-tially the Rayleigh-Ritz variational principle applied to a QFT.Two distinct variants of TSMs exist, one where the standard equal time quantization is employed [6,7], and where lightcone quantization is employed [8,9].We focus here on the first.TSM is a numerical approach well suited to study the low-energy spectrum as well as vacuum expectation values and matrix elements of operators.As such, it aims to be an alternative to lattice Monte Carlo, especially in situations where the latter is plagued by the sign problem.
In particular, TSMs are well-suited to study the spectrum and matrix elements of bound states and topological excitations.For example, the kink-antikink elastic scattering phase can be obtained in the the strongly coupled regime of the broken phase ϕ 4 model [56].TSMs are also an ideal tool to study QFTs subject to sudden (or gradual) changes of the couplings in the Hamiltonian or the environment, called quenches [37,38,42].Real-time correlators in the presence of these quenches are available to be computed with TSMs [59].Periodic driving and Floquet dynamics was analyized in [72].False vacuum decay was studied in [45].
However, it is generally difficult to implement Hamiltonian truncations for systems featuring more than one space dimensions, multiple fields and/or gauge constraints.In this paper, we present an incarnation of Hamiltonian truncation that aims to improve on the majority of the aforementioned issues.
TSMs are typically applied to a problem where the Hamiltonian comes in two parts: Here Ĥ0 is some Hamiltonian over which one has complete control (i.e.Ĥ0 is that of a conformal field theory with known structure constants or of a free boson or fermion).The space of eigenstates of Ĥ0 often provides the computational basis under which TSMs play out.The second term of the Hamiltonian, V, is some interaction whose effects one is trying to understand.To begin to apply the TSM to understanding this problem, one first makes a division of the computational Hilbert space, H, into two parts, In the simplest (and original) form of TSMs, this division is predicated upon energy as determined by Ĥ0 .In this particular case, H 1 then consists of a finite basis of emphasized states and the effects of the (infinite dimensional) H 2 are ignored entirely.For example in the study of the relevant perturbation of the Lee-Yang conformal minimal model [6], the computational Hilbert space, H consists of the Lee-Yang CFT space of states, H 1 consists of conformal states with energy (as defined by the conformal Hamiltonian) less than some cutoff, and H 2 all states above this cutoff.In this way the TSM reduces the solution of a quantum field theory to an exact diagonalization.More sophisticated versions of the TSM build on this crude (but surprisingly often accurate) approximation by computing how the discarded part of the Hilbert space, H 2 , affects the physics.Taking into account the effects of H 2 can be done purely numerically in the spirit of a numerical renormalization group [73] or analytically where the effects of integrating out H 2 leads to a new effective Hamiltonian defined on H 1 [1, [74][75][76][77].It is at this stage, in how to deal with H 2 , that we offer a novel approach.
In most applications of TSMs, the dimension of H 1 is on the order of 10 3 to 10 7 .Here we instead explore keeping the dimension of H 1 small, less than 100.Keeping dimH 1 small mandates that we account for the effects of H 2 .Here we propose to do so by using a Krylov subspace iterative approach.In principle, this iterative approach can be taken to arbitrarily high order and so offers the ability to compute quantities to arbitrary precision.While dimH 2 = ∞, at each Krylov order we do not need to introduce an explicit UV cutoff on the Hilbert space.We are able to avoid doing so because we compute the matrix elements that arise at a given Krylov iteration by employing a Feynman diagrammatic representation that is then evaluated with Monte Carlo techniques.In general these diagrams can be regularized with conventional UV cutoffs from the perturbative literature.In particular, in the 1 + 1d ϕ 4 model they are UV finite after normal ordering and do not require the introduction of a cutoff to be evaluated.Importantly each Krylov iteration is variational and is guaranteed to improve upon the previous iteration.The first Krylov iteration is akin to the NLO approach of Elias-Miró et al. [1].
The paper is organized as follows.In Section II we begin with an overview of the truncated spectrum methodology and its bipartition of the computational Hilbert space.We show how the notion of a 'tail' state, first introduced in a TSM context by Ref. [1], arises from this partition.We then show how the tail states can be computed through an iterative, continued fraction approach.This approach involves writing the tail states in terms of a Krylov subspace formed by the action of products of a resolvant operator T. We give a variational estimate to each iteration of the continued fraction, which is in turn a variational estimate of the ground state energy of each symmetry sector, i.e., guaranteed to form an upper bound on the exact energy, and that each iteration of the continued fraction improves upon that before.The discussion in Section II is model agnostic.
In Section III we then turn to applying it to the ϕ 4 model in 1 + 1d.There we provide a brief review of the properties of the model in both its unbroken and broken phase.We then turn to defining the computational space of states and the division of the Hilbert space, H into two parts, H = H 1 ⊕ H 2 .We discuss in particular a division termed "Zero-mode separation", where H 1 consists of the zero-mode eigenstates that arise from a ϕ 4 Hamiltonian shorn of its oscillator modes.
Having established the computational basis that we will use, we turn to numerical details on how the TSM Hamiltonian is evaluated in the Krylov subspace.The first detail that we discuss is how matrix elements of the Hamilitonian are computed.Here we show that these matrix elements can be written in terms of standard ϕ 4 Feynman diagrams and that these diagrams in turn can be evaluated through Monte Carlo.Formally, the exact ground state energy E * = E appears explicitly in the Krylov Hamiltonian.Thus in principle finding the eigenenergies of this Hamiltonian must be done in a selfconsistent manner.However we will argue that E * is essentially a parameter of the Krylov basis and can be traded for another set of states with no explicit dependence on E * .This will allow us to choose a basis arising as a natural tensor product of the zero-mode states and the states with oscillators.
Having finished with how we will implement the method in the context of ϕ 4 , we then turn to results.In Section IV we discuss results for the broken phase of the model.There we present results for the kink mass and the value of the critical coupling.In Section V we present results on ϕ 4 's disordered phase, both the ground state energy and the mass of the first excited state.In Section VI we demonstrate that our method likely provides a better estimate of the low-lying eigenvectors in the theory (as opposed to eigenenergies) than traditional truncated spectrum methods.
In the final section, Section VII, we turn to the discussion of the method.We consider two points.In the first, we can in principle compute the method at arbitrary high Krylov iterative order.However in practice, we are limited by the ability to perform the numerical evaluation of n-point functions of the ϕ 4 interaction.In this paper the highest Krylov order that we go to is 3.While we have some understanding of how to extrapolate in Krylov order, this issue remains in general open.In the second, we discuss the prospects of applying the method to higher dimensional field theories.We believe that the method, because it employs Feynman diagrammatic representations, has a good chance of being generalized to at least linear sigma models in 2 + 1d.We consider this to be an exciting possibility.While applications of TSMs to 2 + 1d QFTs have only recently appeared, TSMs have been applied to 1+1d QFTs for over three decades.

II.1. Overview
Quantum field theory (QFT) is essentially quantum mechanics (QM) combined with Lorentz invariance.More precisely, a QFT can be imagined as a limit of of a QM system as the number of degrees of freedom is taken to infinity.An immediate consequence of this is that numerical methods developed for quantum mechanics are potentially applicable to QFT problems as well.
One such method is the Rayleigh-Ritz-Löwdin approach coming from quantum chemistry [78][79][80].In the context of quantum field theories with infinite dimensional Hilbert spaces, this approach often comes with the moniker 'Truncated Spectrum Methods' or in its more narrow application to perturbed conformal field theories in 1 + 1d, 'Truncated Conformal Spectrum Approaches' [6,48].Let us write our infinite dimensional Hilbert space in which our quantum field theory operates as where we suppose the |c i ⟩ are orthonormal.Now let Ĥ be a Hamiltonian defined on this space.Because H is infinite dimensional, it is not possible numerically to find the exact eigenstates {|ψ i ⟩} ∞ i=1 and their corresponding eigenenergies E i : To circumvent this problem, one can divide the Hilbert space into two parts: Here H 1 is a finite dimensional subspace and has been chosen to contain, in some sense, the most 'important' states in the basis {|c i ⟩} ∞ i=1 for the problem.We will return to this notion of 'importance' later in the paper.
With this division, the Schrödinger equation for our Hamiltonian becomes, in block form: Here the vector (c 1 c 2 ) should be understood, schematically, as a vector in the bipartite Hilbert space: The simplest approximation that we can make is to suppose that Ĥ12 and Ĥ21 are small and thus can be ignored.
We then can solve the restricted, numerically finite problem, of finding the eigenvalues of Ĥ11 : Often this approximation is excellent.This can be seen in problems where the Hamiltonian comes in two parts: Here Ĥ0 is, typically, diagonal in terms of the basis while V is an operator that mixes this basis.With the Hamiltonian in this form, it is natural to make the division between H 1 and H 2 on the basis of the states |c i ⟩ as determined by Ĥ0 : H 1 is composed of a finite set of low energy states: while H 2 consists of the complement: This form of the Hamiltonian was present in two of the first uses of truncated spectrum methodologies to study properties of perturbed conformal field theories, [6,48].
Here the dimension of H 1 was kept to be very small, a few tens, yet the answers generated were accurate to several significant figures.Here the reason why the methodology worked so well was because the operator V was highly relevant, i.e.V12 was in some sense small and only mixed the two parts of the Hilbert space weakly.But it will often be the case that one wants to understand theories where V12 is not small.In cases where it is not possible to ignore the mixing of H 1 and H 2 induced by Ĥ12 / Ĥ21 , we need a strategy to incorporate the effects of this mixing.Such strategies can be numerical or they can be analytical.Numerically, it is possible to adapt a Wilsonian numerical renormalization group to the problem [73].Analytically, a common way to improve precision is to perform perturbation theory to leading order in V: This is of course not standard perturbation theory in V.This operator's effects in mixing states in H 1 are being accounted for exactly.What is instead being treated perturbatively in V is the mixing between H 1 and H 2 .
While computing this correction can sometimes improve the results, often it does not.In response, one could imagine going further and computing higher corrections in the mixing between the two parts of the Hilbert space.However, it can be that the perturbative series in V is asymptotic.This is the case in studying the ϕ 4 perturbation of a scalar real field in 1 + 1d [1].It is also the case in applying truncated spectrum methods to the sinh-Gordon model [60].We thus want to develop a framework where these corrections can be analyzed systematically and sense made of any asymptotic series.
The first step is to recast the exact problem so that it a finite dimensional one.The full Schrödinger equation reads when broken into components as We can eliminate |c 2 ⟩ from this equation leaving us with a single Schrödinger equation for the first N 1 ≡ dim H 1 eigenvalues and eigenstates of the problem: While we have recast the problem into a finite dimensional one, it is still in general intractable because we cannot compute the matrix elements of ∆H.
To begin to attack this issue, we take the stance of Ref.
[1]: we recognize that if one is interested in the obtaining the exact values of the first N 1 eigenvalues of the problem, one does not need to deal with the infinite part of Hilbert space, H 2 , in its entirety.Rather the subspace of H 2 defined by is sufficient.Following Ref.
[1], we call the states in this span, 'tail' states.
Considering the problem in the expanded (but still finite dimensional) space H 1 H t , we can write the Hamiltonian in the form where Pt is a projector onto H t and Ĝtt is the inner product matrix of the tail states -while the basis in H 1 is orthonormal, the tail states are not so.While this recasting of the problem does not in itself solve our problem -not being able to compute the matrix elements of ∆H amounts to be unable to write down the exact form of the tail states -it does provide us with a basis for finding systematic, iterative approximations for the tail states.This we do in the next section.

II.2. Tail States as Continued Fractions
Here we construct in an iterative fashion expressions for the tail states that have a continued fraction form.The basis on which we develop the continued fraction will be to exploit the division of the Hamiltonian into a 'free' part Ĥ0 and an interacting part V as in Eqn.(9).Our presentation here is based on Ref. [81].
The exact tail states, can be rewritten implicitly as We can think of |t l,0 ⟩ as a zero-order in V22 approximation to the full tail state |T l ⟩ while T |T l ⟩ encapsulates the first and higher order terms of |T l ⟩ in V22 .Our goal here is to develop such terms systematically.
To do so, we want to define an operator, Tl,1 , that separates from |T l ⟩ the contribution of |t l,0 ⟩, i.e.
Tl,1 |t l,0 ⟩ = 0. ( The operator that does this is We can use this definition to write the exact tail state as |T l ⟩ as Eq. ( 21) is almost what we are looking for in terms of isolating the next correction beyond |t l,0 ⟩ to the tail state |T l ⟩ in the form of a continued fraction.The only problem with it is that |T l ⟩ appears on both sides of this equation.
To evade this, we multiply eq. ( 21) with ⟨t l,0 | T from the left and rearrange terms to yield where we have defined |T l,1 ⟩ as At the next order of approximation, we can write giving us an approximate for the exact tail state |T l ⟩: We stress that the correction to |T l ⟩ beyond |t l,0 ⟩ is not merely one higher order in V but is some resummation of terms in powers of V. Unlike a naive perturbation theory, this expansion is well-behaved, i.e. is not asymptotic.We can continue this expansion by repeating the process with |T l,1 ⟩.We define an operator Tl,2 that projects from |T l,1 ⟩ the contribution of |t l,1 ⟩: As before we can then write Here we have the option of approximating |T l,2 ⟩ as |t l,2 ⟩, valid if We can continue this process iteratively, defining a sequence of states |T l,n ⟩ , |t l,n ⟩ and operators, Tl,n defined as This sequence describing the exact tail states can be terminated at finite order by making the approximation

II.3. A Variational Improvement
The takeaway point of the above calculation is that we have recast the exact tail states |T l ⟩ as a linear combination of states with basis vectors In the above, we have truncated our approximation for the tail states at order N T .The iterative procedure described in the previous section gives the form of the nonperturbative coefficients τ lk .However rather than compute the coefficients out to some order, we are going to treat them as variational parameters that are chosen so as to minimize the energy.We do so by starting from a tail-extended basis given by the combination of the original low energy Hilbert space, H 1 , with | tlk ⟩: This computational basis has dimension N 1 (N T +1  (11) G (12) . . .G (1N T ) 0 G (21) G (22) . . .
Recalling that the operator T was defined in eq. ( 18) as the matrix elements of the block matrices H 11 , H 1t , H (nm) tt are defined as and the non-trivial elements of the Gram matrices are defined by II.4.The NLO Approximation of Ref.
[1] and Beyond As we have stated, Ref.
[1] introduced the idea of working with a basis formed from H 1 and the span of the tail states, H t .However Ref. [1] did not attempt to construct the exact tail states of eqs.( 22) and (28).Instead they examined the leading approximation to the tail states of the form and set up the truncated Hamiltonian as They termed this approximation 'next-leader order' (NLO).In this paper we will also explore variational approximations that go two orders beyond this -taking the eigenvalue problem in (33) with N = 2 -we term this 2NLO.

  
H 11 H (1) 1t and N = 3 (we term this 3NLO) (11) G (12) G (13)   0 G (21) G (22) G (23)   0 G (31) G (32) To this point, this section has presented a general for-malism independent of particular models and particular choices of H 1 and H 2 .For a model with the simple structure of eq. ( 9), matrix elements of H (nm) tt,ll′ involve n+m+1 powers of the interaction: at most 3 at NLO order, 5 at 2NLO and 7 at 3NLO.In the next section, we turn to these issues in the context of the ϕ 4 scalar field theory.

III. IMPLEMENTING THE METHODOLOGY
ON THE ϕ 4 THEORY The discussion in Section 2 was completely general and was presented independent of any model.It also did not specify certain crucial implementation details that would be needed in any actual application of the methodology.In particular, it did not indicate how we intended to divide the computational Hilbert space into H 1 and H 2 .And we did not discuss how we are going to compute the matrix elements of H

and H
(nm) tt .We take up these tasks in this section in the context of the ϕ 4 theory.

III.1. Review of ϕ 4 Theory
The two-dimensional ϕ 4 theory corresponds to the normal ordered Hamiltonian where g 4 > 0 and the normal ordering is understood with respect to the mass scale m 0 and at infinite volume.The Hamiltonian is invariant against the global Z 2 transformation φ(x) → − φ(x) in addition to Poincaré invariance and spatial parity symmetry.Classically, separating the quadratic term into two pieces is of course redundant.
In the quantum theory, different choices for Ĥ0 are connected by finite a renormalization of the couplings.
In finite volume, L, Lorentz symmetry is lost but translation invariance persists.The periodic finite volume mode expansion of the field φ takes the form The Hilbert space in which Ĥ acts can be built up from the vacuum |0⟩ using the creation operators â † n .A general Fock basis vector takes the form In practice it is more feasible to work with operators normal ordered at finite volume L. Strictly speaking this yields a different regularization scheme.The schemes can be connected by applying finite, volume-dependent corrections to the coefficients g 2 and g 4 .A derivation for these corrections is provided in Appendix B. The explicit form of the Hamiltonian, expressed with finite volume normal ordering, reads where we introduced The functions z(x) and e 0 (x) are defined in eqs.(B10) and (B14), respectively.The model possesses two phases.Classically, for g 2 > 0, the ground state is invariant with respect to the parity Z 2 symmetry, while for g 2 < 0 the symmetry is spontaneously broken.In the quantum model, the situation is more complicated due to the presence of the additional energy scale m introduced by normal ordering.In particular, a duality emerges between a theory with g and the relation holds [55].This is called the Chang duality and is a weakstrong duality in the quartic coupling.A consequence is that the Ising critical point and the broken phase can be reached starting from the symmetric phase and increasing the quartic coupling g 4 at fixed g 2 .We provide a collection of estimates for the location of the critical point in the literature in Table I.

III.1.1. Symmetry preserving phase
In the symmetry preserving phase, there is a single, self-interacting excitation that correspond to the elementary fluctuations of the field.The bulk energy density, E 0 , has been calculated to 8 loops in [19], yielding This is an asymptotic series, which can be resummed using Borel techniques.
The physical mass M is given by the energy difference of the lowest-energy one-particle state and the vacuum.It can be obtained from the position of the pole in the interacting two-point function, which, in turn, is accessible through a coupling expansion.The expansion is given as which is again a Borel resummable asymptotic series.

III.1.2. Symmetry broken phase
In infinite volume, where symmetry breaking is exact, there is a doubly degenerate ground state.In finite volume, instanton effects restore the symmetry, but the energy splitting of the two vacua remains on the order of e −M K L , where M K is the kink mass, given semiclassically (g 2 < 0, |g 2 | >> g 4 ) as [88,89] (50) The broken phase is enriched by a number of mesons, resonances and topological kinks: The exact structure of the spectrum also depends on the coupling.The bulk energy density in the broken sector can be calculated by expanding in the fluctuations around either vacua.Its perturbative expansion in the coupling starts as [20] Λ This series is also Borel resummable, but the naive series provides a robust estimate in the weakly coupled branch (g 2 < 0, |g 2 | >> g 4 ) of the broken phase.

III.2. The Computational Space of States
We now turn to how we will set up our computational space of states.The first step in doing so is specify how we are dividing the Hilbert space of the ϕ 4 theory into two, i.e., how we are choosing H 1 and H 2 in eq. 5.
H 1 consists of states in the zero mode sector of the theory, i.e., states involving no oscillator modes, ân , n ̸ = 0. To explain this further, we follow [55,56] and single out the zero mode part of the field φ and its conjugate π: In eq. ( 52) and in the following, we apply the following notation.For an operator A acting only on the zero mode subspace, we omit the hat.For an operator B acting only in the oscillator subspace, we apply a wide tilde.Whenever a hatless operator is added to a wide tilded operator, it is understood that they are tensored with the unity in the other subspace, so that the addition is meaningful.
We then divide the Hamiltonian Ĥ into a zero mode H ZM and non-zero mode ĤNZM part: with the free parts defined by There is some freedom in how one distributes V between V ZM and V′ which we exploit in this work.One choice for this division is given in Section III.2.1.We first describe our choice for H 1 and H 2 as described in Eq. ( 5).Our choice for H 1 is given by In We now turn to specifying our tail state basis.For each |m⟩ ∈ H 1 , we define a sequence of tails given by (see Eq. 31) where P ⊥ is a projector onto H 2 .

III.2.1. "Conventional" tails
An obvious choice for separation of V into V ZM and V′ is to take where we used the notations In the above g ′ 2 and E 0 (m 0 , L) are defined in Eq. ( 46).In order to complete the computation we need to specify the matrix elements of our basis.At first tail order (see Eq. 38), we can substitute the projected operators with the unprojected ones in the relevant matrix elements.This is because the extra terms induced by the projectors always involve one-point functions of operators : ϕ n : for some n ≥ 1, which are all zero.We obtain the following matrix elements (G (11) ) We provide similar expressions for the matrix elements involved in the eigensystem at second tail order in Appendix C (see Eq. 39).In turn, these expressions can be transformed into Euclidean multipoint functions of the interaction V′ by using the integral representation In particular, for the matrix element (H In order to implement our scheme we need to make a choice for N ZM .This choice is unproblematic.The lowest eigenvalues of the full problem converge quickly in the number, N ZM , of kept zero-mode states.Typically N ZM can be kept to a number less than 10.And at least in certain cases, the zero-mode Hamiltonian, H ZM , encodes much of the physics of the full problem.In such cases the effects of the oscillator modes on the physics of the problem can be considered as perturbative.We have found that this choice of separation of H 1 and H 2 works in both the unbroken and broken phases of the theory.In Fig. 1 we provide an example of our use of these tails.There we compute the ground state energy in the unbroken phase up to the third tail order for a non-trivial value g 4 as a function of N ZM .[90] We see that after N ZM exceeds 6, the result reaches an asymptotic value.

III.2.2. Using numerically efficient tail states
The vanilla method described in the first part of this section suffers from three technical drawbacks.Firstly, the computational basis is not a simple tensor product of states in H 1 and states built by acting with oscillator modes above the oscillator vacuum.The evaluation of matrix elements using Monte Carlo results in a linear algebra that scales cubically with the number, N ZM of zero mode states kept.Secondly, the tail states depend on the coupling g 4 in a complicated way and so a separate computation is needed for each value of g 4 to be investigated.Thirdly, the tails depend explicitly on the exact energy level E * , which has to be tuned self-consistently, and is in principle different for the ground state and for excited states.These three drawbacks make the method as described computationally expensive.
We can avoid these issues if we define our tails with the following choice for V (U ) ZM and V′ (U )  .The superscript (U ) appears here to emphasize a different separation of V as compared to eqs. ( 57)- (57).We choose We also alter our choice for our computational basis.We change H 1 to Here E ZM p = m 0 p.Now H 1 consists of a basis of free oscillator states |p⟩.
Neglecting for a moment the cutoff, N ZM , in the zeromode subspace, tails of order K obtained with V (U ) ZM of eq. ( 62) can be written as a linear combination (l j ≥ 1) where c pK p ′ Q,{k},{l} are numerical coefficients and | 0⟩ is the ground state of H (m) 0,osc .Rather than work with the tail states | tpK ⟩, we will use the larger set |p ′ ⟩ ⊗ | t Q,{k},{l} ⟩ directly as a basis.One can see that this basis choice is a tensor product of the zero-mode space and the space of oscillator states.This will lead to a dramatic simplification in our numerical evaluation of matrix elements.
Using states of the form |p ′ ⟩ ⊗ | t Q,{k},{l} ⟩ represents a change of basis but it is not expected to spoil the convergence of the method (see also Appendix A).In fact this extension of the basis improves on the accuracy through the increase of the dimension of the variational basis.However one issue seen with this choice of basis is that the Gram matrix rapidly becomes degenerate, particularly at higher Krylov orders.To alleviate this issue we will impose an extra restriction on the tails and use only l j = 1 states (i.e.states that involve only single inverse powers of H (m) 0,osc in our numerical basis).We will refer to the resulting subspace as the universal tail space.
The further reduction of this tail space to one spanned by tails with at most K denominators will be called the universal tail space of order K. Accordingly, when it does not lead to confusion, henceforth we simplify the notation of the l j = 1 tails, specifying them by the list of operator powers k l only: | t Q,{k},{l} ⟩ ≡ | t {k1...k K } ⟩.Note that the resulting basis still has significantly higher dimensionality than the number of conventional tails.At order K, The values of the ground state energy for the first three tail orders computed at g4 = 1.5, L = 10 using the approximate tails for different choices of NZM .We again compare our results with those found in [1].Blue, orange, green dots correspond to 1, 2 and 3 Krylov orders, respectively.Also shown is the power law extrapolation using statistically generated eigenvalues (red) and its empirical improvement (purple).We see the results rapidly converge as a function of N kept .We compare our results to data provided to us by the authors of [1] (here shown as a dashed black line with error bounds shown as a thick gray line around it).Here NMC,0 = 10 8 , Niter = 50.Overlaps of normalized higher denominator tails with the universal tail space of order 2 (see main text).Even parity sector, m0L = 10.Tails | t Q,{k},{l} ⟩ corresponding to numbered points are as follows: {44},{21} ⟩.Note that points 3 through 14 only appear in K = 3 tails.
the total number of conventional tails is KN ZM , while the dimension of the full universal tail space truncated at order K is We reproduced the quantity depicted in Fig. 1 with universal tails in Figure 2. The error bars are smaller in the latter plot, as more integrand evaluations were possible due to the simpler numerics.In turn, this made possible to explore various extrapolation approaches in Krylov orders, which we discuss further in Sections III.3.3 and V.1.
With higher denominator tails omitted, the method still remains variational, but in principle prone to a small systematic error due to the neglected basis states.We show in Figure 3 that, after normalizing to unity, the overlaps of all l j > 1 tails appearing in eq. ( 64) up to K = 3 on the universal tail space of order 2 is close to 1.While Fig. 3 includes tails from only the even parity sector, the situation is entirely analogous in the odd sector.Thus the neglected states lie almost entirely in the space of states that we do keep.Further, since the tails | t Q,{k},{l} ⟩ are independent of the couplings, so are their overlaps.Therefore, no significant errors are expected to arise from omitting the higher power denominators.
As another test of this modified computation basis, we test our ability to represent the "conventional" tails in (56) (with V ZM defined in eq. ( 57)) in our new basis.In Figure 4, we see that the overlaps of first order "conventional" tails are close to one.These overlaps are defined as, where the indices I, J run over all universal tails up to order K, and G U is the inner product matrix between universal tails, This remains true for a wide range of couplings, E * parameters and volumes.This completes a numerical demonstration that replacing our original tail states with |p⟩ ⊗ | t {k} ⟩ is a very good approximation.
A particularly suitable form of the Hamiltonian is gained by emphasizing its tensor product structure relative to this new basis.Introducing the shorthand V n = L 0 : ϕ(x) n : m0,L dx, the Hamiltonian Ĥ can be written as with V ZM as defined in eq. ( 57) and g n defined in eq.(58).In our new basis, the needed matrix elements can be reduced into matrix elements of V n and 1 between the universal tails: where the oscillator space matrix elements take the form In the last line we transformed the energy denominators into multi-point functions with the aid of eq. ( 60).Note that the oscillator space matrix representation of H (m) 0,osc and the inner product matrix are actually diagonal between first order universal tails.The use of the basis vectors (64) in calculating the matrix elements of the effective Hamiltonian reduces to calculating general n-point functions of : ϕ n :.Matrix elements involving higher tail orders contain disconnected pieces.As with the original tail states, we provide the explicit expressions for the matrix elements needed to carry out the computation with second order tail states in Appendix D.
We show a demonstration of the precision that we obtain with these approximated tail states in Fig. 2. One can see that both the energy estimates and the error bars have improved compared to Fig. 1.The improvement in computed energies is due to the larger variational basis.The error bars have decreased in part due to auspicious error propagation from matrix elements to the lowest eigenvalues, in part due to cheaper evaluation costs as we were able to increase the number of integral evaluations by an order of magnitude.
The matrix elements that we need to evaluate can be represented in terms of Feynman diagrams.To see this, we introduce the multi-index quantity  4. Overlaps of normalized conventional first order tails in the broken sector with g4 = 1.5, mL = 10, and E * = −0.878,with the universal tail space of orders K = 1 (blue dots) and K = 2 (orange dots).Tails built upon the 10 zero-mode eigenstates with lowest energy are shown sorted by increasing ZM energy.(71) for an (n+1)-point function (n ≥ 1) of the integrated field powers.In terms of V in,••• ,i0 , the matrix elements can be written as In turn, the multi-point functions can be evaluated with Wick's theorem, yielding a diagrammatic structure similar to a perturbative calculation.We list the simplest matrix elements below, while the general Feynman rules are discussed in Appendix 3. V244, (e): V334, (f): V444 (see also eq. ( 73)).Permutations of the indices lead to the same diagrams with permuted τ arguments.
In the next subsection, we examine the explicit form of the propagator ∆(2) L , given in eq. ( 78) below.The diagrammatic representations of certain representative V in,••• ,i0 are given in Fig. 5.

III.2.3. Finite Volume Propagator of ϕ 4
In our implementation of Krylov space methods for ϕ 4 , we simultaneously use both the Hamiltonian formulation of ϕ 4 together with a Feynman diagrammatic description.We need to set the conventions of the the finite volume propagator that appears in such diagrams.
For the moment we will work in general space-time dimension.Our action for the Lagrangian of a free massive scalar field is The corresponding two-point function in infinite volume then takes the form: where T τ indicates ordering with respect to an arbitrarily chosen Euclidean time direction τ , r = |x 1 − x 2 | is the Euclidean distance, and K ν (x) is the modified Bessel function of the second kind.
In the special case D = 1, eq. ( 74) describes a single harmonic oscillator with frequency m 0 and mass 1.The propagator D (D) from eq. ( 75) reduces to As our computations for D = 2 will be done in finite volume, we need the corresponding finite volume propagator.We will assume that we have either periodic (P=0) or anti-periodic (P=1) boundary conditions.As derived in Appendix E, in D = 2 the finite volume mode expansion leads to the propagator L,∆ (τ, x) As we can see, in finite volume, the propagator ∆ (2) acquires corrections that are suppressed as e −m0L .When the zero mode is separated and the multipoint functions are restricted to the oscillator sector, we are required to substitute the complete propagator in eq. ( 77) with a modified one, In practical computations, we truncate the sum in eq. ( 77) to n ≤ 3, which is a good approximation for L > 3.
Finally, we note that for small volumes, it is possible to work in a different scheme where the unperturbed oscillators are massless.In this case the oscillator mass term is considered part of the perturbation and the modified two-point function is given by the closed formula The benefit of the massless scheme is that the volume dependence of the matrix elements factors out, so the numerical integrations do not have to be repeated for each studied volume.However, we found the precision of massless setup strongly volume-dependent and suboptimal for the volume range of our interest.Thus we do not elaborate upon this direction further.

III.3. Numerical Implementation
Now that we have specified the computational bases that we will use in evaluating ϕ 4 , we turn to a discussion of the details of the numerical evaluation of the eigenvalues of the Hamiltonian.
The numerical method proceeds as follows.We first compute the matrix elements of Ṽn and 1 of eq. ( 68) between order one universal tails.We then assemble the full Hamiltonian with a zero mode cutoff N ZM and solve the generalized eigenproblem in this 4N ZM dimensional basis.In the next step, we calculate all matrix elements involving universal tails of second order, and repeat the diagonalization in this larger, 13N ZM dimensional basis.Finally we also include third order tails and obtain the eigenpairs of the resulting matrices.We thus get a series of eigenvalues as the function of the largest order of tails kept.This set of points, possibly supplemented by a "zero-order" approximation with all tails neglected, serves as basis to extrapolate to infinite Krylov order.It is useful to note that both the conventional and universal tails are Z 2 parity eigenstates, so the even-odd particle number separation can be done even in the tail Hamiltonian.
For all couplings examined, eigenvalues converge exponentially quickly in increasing N ZM .This is exemplified for g = 1.5 in Figure 2. We found that the regimes most sensitive to the choice of N ZM are when either g 2 is large and negative and g 4 is small such that ϕ acquires a large vacuum expectation value, or when g 4 is very large.On the basis of these observation, in the following we fix N ZM = 40.
The details of this method can be grouped into three areas: i) the evaluation of the matrix elements through Monte Carlo; ii) generating a statistics of eigenvalues to estimate their errors, iii) the procedure by which we extrapolate to infinite Krylov order K = ∞.We take each in turn.

III.3.1. Evaluation of Matrix Elements
Matrix elements involving two-dimensional integrals are evaluated using the numerical global adaptive method Cuhre, in the implementation of the Cuba package [91].These matrix elements consist of H (1) 1t and G 11 .Their estimated relative numerical errors are negligible, on the order of 10 −8 .In turn, matrix elements requiring the evaluation of at least four-dimensional integrals are computed with the importance-sampling Vegas method of the Cuba library.This is basically a Monte Carlo method, in which the sampling ensemble is updated selfconsistently in a number N iter of iterations.Inside each iteration, a Monte Carlo-type probing of the integrand is performed with N M C,0 evaluations.The total number of integrand evaluations is thus N eval = N M C,0 • N iter .It is generally efficient to choose the value of N iter to be a few dozen [92].We fix N iter = 50.As a general rule, larger dimensionality of the integrand comes with slower convergence with respect to Monte Carlo evaluations.
The integrations are take place over a finite domain.The restriction in space occurs naturally because of our working in finite volume.We however introduce a "temporal" cutoff, fixed to be τ max = 50m 0 for V i1...in with up to 3 indices, and tau max = 40m 0 for those V i1...in 's with more than 3 indices.These choices are made so that the results are insensitive to further increasing τ max .This is reasonable as the integrand vanishes exponentially for large values of each of its τ arguments, which in turn follows from the asymptotics of the propagators and the presence of disconnected terms.
It is easy to check the precision of numeric integrals by imposing an artificial momentum cutoff p max .In the presence of a small enough p max , the required matrix elements can be exactly computed in an alternative way: building up the matrix representations of the operators in eq. ( 35) directly in a truncated Fock space with appropriate particle number and one-particle momentum cutoffs.On the other hand, such a cutoff is easily implemented in the integrals by replacing the propagator via restricting the sum in the mode expansion eq.(E3).The dependence of numerical error on the number of evaluations for the matrix element ⟨ t 444 | V j | t 44 ⟩ for various values of p max is shown in Figure 6.Other matrix elements behave similarly (or better).

III.3.2. Statistics for third order tail states
With the introduction of higher order tails, the inner product matrix becomes approximately degenerate.This issue is especially pressing in the case of universal tails.At the third tail order, the numerical precision of matrix elements is not sufficient for a direct computation, which manifests in (small) negative eigenvalues of the Gram matrix.
To alleviate this problem, we generate an ensemble of Ĥ and Ĝ matrices using the numerical integration errors of matrix elements and repeat the eigenvalue calculation several times, projecting out the directions corresponding to numerically negative eigenvalues of the Gram matrix.The result of this analysis for a sample size N s = 500 is shown in Figure 7.For each peak, we then refine the 0 1 × 10 9 2 × 10 9 3 × 10 9 4 × 10 9 5 × 10 9 -0.0010 -0.0005 0.0000 0.0005 0.0010 Relative error of the matrix element ⟨ t444| Vj| t44⟩ as function of the numerical integrand evaluations, in the presence of small momentum cutoffs: pmax = 2πL −1 (blue), pmax = 4πL −1 (green), pmax = 8πL −1 (red), at volume m0L = 10.The differences are measured and normalized with respect to the exact matrix element obtained by linear algebra in the explicit Fock basis for these cutoffs.The importance sampling method uses Niter = 50 iterations of updating of the sampling ensemble with NMC,0 = 10 8 evaluations inside each iteration.
energy bins until the maximum of the peak becomes less than ∆N ∆E treshold = 60.Separate Gaussian fits are then applied to the peaks.The standard deviation of these subsequent fits provide the errors of the eigenvalues at third order.Distribution of the lowest eigenvalues for 3 tail orders, g = 1.4 m0L = 10, based on an ensemble of Ns = 500 Hamiltonians, sampled assuming normal distributions for the matrix elements with standard deviations from errors derived from the Monte Carlo integration routines.The histogram is based on calculating the lowest 10 eigenvalues per parity sector for each matrix in the ensemble.Blue dots correspond to the even Z2 parity sector, while orange dots correspond to the odd sector.Energy bins of size ∆E = 0.005m 2 0 were used.

III.3.3. Extrapolations in K
After the eigenvalues from the three Krylov orders are obtained, we perform a power-law of the form As the input has an error, we generate a new ensemble of eigenvalues with the previously obtained errors.The fit is repeated for each member i of the ensemble, providing a set of extrapolations a i .a i are then averaged and their computed standard deviation provides the final error estimate (red error bars in Figure 2).The result of this extrapolation fits nicely to previous TSM results [1].

IV. BROKEN PHASE: "WEAK" COUPLING
We start our detailed presentation of our results for the ϕ 4 theory by focusing on its broken phase.We fix g 2 /m 2 0 = −0.25 and vary g 4 so that an Ising type quantum critical point is present at a relatively small quartic coupling.The spectrum and various physical quantities are well described by the perturbative treatment of fluctuations around the degenerate vacua.We consider the ground state energy, the meson mass gap, and the vacuum expectation value of the field, ⟨ϕ⟩, before turning to a finite volume scaling analysis to determine the location of the critical point.

IV.1. Ground state energy
Figure 8 shows the coupling dependence of the bulk energy (ground state energy per unit length) as a function of the quartic coupling g 4 = gm 2 0 for different system sizes L compared to the 4-loop perturbative result in infinite volume.There are two kinds of physical finite volume effects that alter the finite volume eigenvalues from their L → ∞ asymptotics.On the one hand, the doubly degenerate vacuum is split due to finite volume instanton effects, which is semiclassically exponentially small in the parameter M K L, where M K is the kink mass (eq.( 50)).On the other hand, there are finite volume corrections corresponding to the Casimir effect, which are exponentially small in the parameter M L, M being the smallest mass in the spectrum.At small quartic coupling, M is the lightest meson mass M 1 and is about unity.In the regime g > 0.1, M is again the kink mass and is approximately order 1/g, before it vanishes at the critical point, giving rise to a stronger, inverse-power volume dependence.In the entire broken sector, the numerics significantly outperforms raw Fock space truncation if no extrapolations are taken into account.On the other hand, we obtain the bulk energy in agreement with perturbation theory, up to finite volume effects.
Since the method does not rely on building up explicit Fock states, it is expected to perform especially well compared to raw truncation at larger volumes.This is indeed true.However, the numerical precision of matrix elements are observed to decrease in increasing the volume, which eventually restricts the applicability of the method in the infinite volume limit.

IV.2. Meson mass
Turning to the lightest meson, M 1 , it is again natural to compare to available perturbation theory results.Here the finite volume effects are largely analogous to the ground state energy, apart from additional exponentially small effects due to virtual particle pairs.F -terms arise due to vacuum polarisation: creation-annihilation of a virtual pair with a nontrivially winding trajectory on the cylinder, thereby scattering on the physical particle.µ-terms correspond to the dissolution of the particle to virtual constituents, which reassemble into the original particle after traveling around the cylinder [67,93].Numerical results are shown in Figure 9.The meson disappears from the spectrum at around g = 0.15, as then it becomes kinematically allowed to decay into a kinkantikink pair.Thus for g > 0.15 in the broken phase, the first excited state is actually a kink-antikink scattering state.[56] IV.3.ϕ VEV Since the method provides approximate eigenvectors, matrix elements of operators can be evaluated directly.
Here we focus on the vacuum expectation value of the field ϕ in the broken sector.In finite volume, this is actually the matrix element of φ between the even and odd parity ground states.To single out the nontrivial quantum corrections, we normalize with the classical expectation value ϕ cl = (8g) −1/2 .As the critical point is in the Ising universality class, ⟨ϕ⟩ is expected to vanish as (g − g c ) 1/8 at the critical point, we plotted its eighth power in Figure 10, against one-, two-, and threeloop perturbation theory.The third order results and the corresponding error bars are again obtained by creating an ensemble of N S = 500 Hamiltonians and Gram matrices as described in Section III.3.2, projecting out any negative Ĝ eigenvalues, and computing averages and standard deviations of ⟨ϕ⟩ over the ensemble.We do not attempt to extrapolate the vacuum expectation value through Krylov orders.The outcome of the analysis is that universal tail sets of orders K = 1, 2, 3 provide results consistent with a perturbative expansion to 1, 2, 3 loops, respectively, provided the coupling is not in the vicinity of the critical point.Around the critical point, finite volume effects dominate.

IV.4. Critical point
The broken and unbroken phases of ϕ 4 model are separated by a second order phase transition.Classically, the phase transition happens when the parameter g 2 changes sign.In the quantum model, fluctuations drive the transition and it can thus be reached via tuning g 4 = gm 2 0 for a fixed g 2 .Despite working in finite volume, where there is no phase transition per se, truncated spectrum methods are effective in measuring the critical g coupling.This is because finite volume effects in the vicinity of a critical point are well understood.
It is well known that the critical point of ϕ 4 theory lies in the Ising universality class.Therefore, in the vicinity of the critical coupling g c , we can model the system as described by the Ising Hamiltonian on the cylinder, modified by relevant as well as irrelevant perturbations: When the g coupling is fine-tuned to the critical point g c , there are still irrelevant perturbations affecting the finite volume spectrum.Fortunately, the irrelevant perturbations can be filtered by their symmetry properties and classified by their relevance.
In the first order of conformal perturbation theory, it is simple to estimate the corrections by the irrelevant terms.Since the operators are multiplied by L −2h+1 , with h being the chiral dimension, the least irrelevant operators dominate the large-volume spectrum, while more irrelevant ones are suppressed by powers of L. In our analysis, we approximate the perturbed Ising model by keeping the relevant thermal perturbation ϵ as well as the descendant ∂ ∂ϵ.We restrict our attention to the energy difference between the even and the odd ground states and consider the volume range L = 6 to L = 10.The rescaled volume-dependent gaps as function of the coupling, for different volumes L, are depicted in Figure 11.The linear coupling dependence of the gap is consistent with the assumption that in the vicinity of the critical point, the effective Ising couplings can be approximated as linear functions of the quartic ϕ 4 coupling g .
In the absence of irrelevant terms, one would conclude that the transition occurs at the point where the rescaled lines cross.But this has to be corrected by the effect of the irrelevant term.To improve the accuracy, we fit linear functions f L (g 4 ) = a L + b L g to the datasets in Figure 11.Equating these empirical gap functions with their theoretical approximations from conformal perturbation theory, and separating coefficients of powers of g 4 , we obtain a (generally overdetermined) system of linear equations for the critical coupling g 4c as well as the parameters α i .There are multiple ways to reduce the resulting system into a uniquely determined one, which provides slightly different results for g c .The variation in these results provides an error estimate.We so obtain an estimate of the critical point to be g c = 0.2645 ± 0.002.The order of performing the extrapolation in tail orders versus taking the difference of eigenvalues does not affect this final estimate.π We present here data for the determination of the critical value of g4.Rescaled E1 −E0 energy differences ("vacuum splitting") for volumes L = 6 (blue), L = 8 (orange) and L = 10 (green).The corresponding value in the Ising model is shown with a horizontal black line.

V. UNBROKEN PHASE AND STRONG COUPLING REGIME
Now we turn to the unbroken phase.Here we choose g 2 = m 2 0 /2 and vary g 4 .As per the Chang duality, the critical point appears at g ≈ 2.7, but higher orders of perturbation theory breaks down noticeably for below the critical point at g ≈ 0.3.In this phase, there is a single bosonic excitation in the spectrum.

V.1. Empirical study of cutoff dependence and tentative improvement on errors
Before turning to the measurement of the ground state energy and the mass gap, let us study the tail order extrapolation in more detail.It is difficult to determine the extrapolating function eq. ( 80) from first principles.However, it is interesting to see how the fitted power law exponent c varies as function of the coupling.We find that performing the fit without any restrictions on the parameters a, b, c, the resulting function c(g) is well described by a simple function form We fitted the parameters c i , i = 1, 2, 3 to the distribution of exponents in Figure 12, obtained by performing the fit (80) on a set of N S = 500 samples generated on the basis of the eigenvalue uncertainties.The proposed functional form eq. ( 82) is empirical and likely model-dependent.At the same time, we find that using the exponent c(g) determined from eq. ( 82) in fitting our data with eq. ( 80) results in good agreement with the data of Ref. [1].In the following plots, we report the results and errors both by fitting directly every sample with eq. ( 82) and alternatively, performing the fit with c constrained to be its eq.( 82) value.The truncated spectrum method (TSM) in its equal time formulation is essentially a finite volume method.The form of leading terms in the large volume expansion of the ground state energy are known.In case of an unbroken vacuum and a single type of massive particle, It has the structure In eq. ( 83), E is the bulk energy (density), M is the mass of the excitation, while the integral term is the first Lüscher (leading exponential) correction.The value of the bulk energy and the mass is taken from the numerics, but knowledge of these parameters provides the volume dependence of the ground state energy E 0 for intermediate to large volumes.This is checked against TSM numerics in Figure 13.
We show the computed ground state energy at volume m 0 L = 10 in the unbroken sector, as the function of the coupling g, in Fig. 14.

V.3. Mass gap
In the top panel of Figure 15 we show the validity of the perturbative expansion eq. ( 48).Consistently with the asymptotic nature of the perturbative expansion, higherorder results come with a more restricted domain of validity.It is interesting to see how this compares to Borel  resumming the perturbative series.While conventional perturbation theory breaks down spectacularly at around g = 0.2, our Krylov-enabled TSM is able to follow the resummed result over an extended range of couplings, as shown on the bottom panel.

VI. RESIDUAL ERROR
In performing any numerical analysis, it is vital to have as much control on the error as possible.In truncated spectrum methods, a commonly adapted measure of precision is to compute the eigenvalues at special points where the results are known "exactly" (with much more precision) from independent approaches.While this can give an idea on the expected performance of the method at large, it has a number of significant drawbacks: Ground state energy in the unbroken sector, m0L = 10.Show are our computations of E0/m0 as function of coupling g compared to [1] (black curve).The numerical results corresponding to K = 1, 2, 3 universal tail sets are marked by blue, orange and green dots, respectively.The power-law extrapolation is shown with a dashed purple line, with errors given by the purple shaded area.Extrapolation instead by constraining the exponent c(g4) to its value given by eq. ( 82) is marked by the data with red dots.
1. Information comes from special points, which are similar at best, but usually not identical to the actual point of interest.
2. Moreover, it is generally hard to quantify the distance from the reference point in theory space, and make a connection with TSM precision.
3. Finally, the error in the eigenvalues is a rather forgiving measure of the error, especially when actual eigenvectors are used to measure a quantity of interest.
In the literature of iterative eigensolver methods, it is customary to introduce the norm R of the residual error in the following way: Note the square inside the expectation value.In eq. ( 84), ψ comp is the computed (truncated) eigenvector, but Ĥ is the entire Hamiltonian without truncation.Without the square, one would simply get zero, as Inserting a partial resolution of the identity between the Hamiltonians results in a strong cutoff dependence, as depicted in Fig. 16.As the action of the operator Ĥ changes particle number by at most 4, it is easy to compute the projected matrix element to high energy cutoffs by introducing an additional particle number cutoff on the truncated space.(Here, for the raw truncation points, we do not separate the zero mode.)It is interesting to note the would-be dimensionality of the full Fock space with the corresponding energy cutoff E cut .The required dimensionality would be D = 6.5 * 10 5 for first (reference data shown as black curve).Numerical results corresponding to K = 1, 2, 3 universal tail sets correspond to blue, orange and green dots, respectively.The power-law extrapolation is shown with a dashed purple line, with errors given by the shaded area.Extrapolation with the exponent c fixed by eq. ( 82) results in the red dots.order tails, about D = 2.3 * 10 7 for second order tails, and D ≈ 10 9 for third order tails, giving a rough estimate of the raw truncation dimensionality needs to reproduce the precision of the Krylov tail numerics.This is impressive in itself, but still suffers from a significant deviation from the exact result.Fortunately, our approach is naturally suited to calculate the matrix elements of Ĥ2 directly, without any inherent truncation.
Let us note that when δE = E comp − E exact is small, the correction obtained from substituting E comp with E exact in eq. ( 84) will be O(δE 2 ), which is negligible compared to the effect coming from the square of the Hamiltonian.Intuitively, the residual error vector estimates the error resulting from acting with Ĥ on the truncated eigenvector instead of its true eigenvector with eigenvalue E exact .This measure can be particularly relevant in measurements which involve acting multiple powers of Ĥ to an approximate eigenvector, something that happens in computations involving the real time dependence of observables.At small coupling, the difference of the eigenvalues between the interacting and the free theory are indeed minuscule.Therefore, relying only on the precision of the eigenvalues is deceptive when the eigenvectors are needed.The residual norm provides a much more sensitive measure of the error.We show the residual error corresponding to the ground state energy in Fig. 17.A reference computation was also performed in a conventional massive Fock-space truncated space with an energy cutoff (with no separated zero mode).In that basis, we evaluated the matrix elements of Ĥ2 by expressing the squared interaction terms with normal ordered ones and performing the nontrivial spatial integrals where necessary.In this case, the residual error decreases approximately linearly with the energy cutoff, at least for the small cutoffs examined.Extrapolating this linear tendency again gives an idea on the Fock space dimensionality required to reproduce the precision of the Krylov computation.

VII. CONCLUSION AND DISCUSSION
In this paper we have implemented a Krylov-subspace method tailored to computing the low lying spectrum of quantum field theories, based on the iterative procedure originally described in [81].It can also be seen as a systematic improvement to the renormalization group method of [1, 54,77].We have demonstrated the method on the unbroken and broken sector of the twodimensional ϕ 4 model.
The method involves the computation of multi-point correlators in the "unperturbed" theory.In the present case, this translates to the numerical calculation of Feynman diagrams.Therefore we see this method as a bridge between perturbative diagrammatics and truncated spectrum methods (TSMs).FIG. 17. Norm of the residual error for the ground state.We point out that the square of the full, infinite-dimensional QFT Hamiltonian is sandwiched between the numeric eigenstates, making the computation difficult for large energy cutoffs.The largest raw energy cutoff corresponds to a 104 dimensional truncated subspace.
Throughout this work we focused on the simplest spectral quantities: the ground state energy and the difference between the two lowest two energy levels.However, we remark that other physical quantities are accessible to the method by modest and often obvious modifications.We succeeded in reproducing the bulk energy and the mass gap of the ϕ 4 model in the unbroken sector in agreement with Borel resummation and other TSM results.Remarkably, in a wide range of couplings, a truncated Hamiltonian of only dimension four (per parity sector) was sufficient to provide a good estimate over a wide range of couplings.In the unbroken sector, this can be contrasted to Borel resummation, which utilizes very similar Feynman diagram computations in an entirely dif-ferent way.We believe it is interesting that even though tail matrix elements suffer from the same asymptotic behavior as the perturbative series, the diagonalization of the Hamiltonian involving these same matrix element efficiently extracts convergent physics from these quantities, in analogy to Borel resummation.In the unbroken section we have also reproduced the correct finite-volume behavior of the ground state energy and demonstrated a novel way to measure the error for TSM eigenvectors.
Using our technique, we also analyzed the ϕ 4 model in the broken sector.We compared the bulk energy to raw TSM truncation and a 4-loop perturbative expansion and found reassuring agreement.We also used our code to provide a precise estimate of the critical point in the broken sector, 0.2645 ± 0.002.
While we have not pursued it here, the method indicates promising features for possible future real-time dynamics applications.Since the same universal tail basis provides a very good representation of lowest energy states over an entire range of couplings, it appears to be a natural basis to study quenches with both suddenly and gradually changing couplings.
This combination of Hamiltonian truncation and Feynman diagrammatics is a new method and our presentation leaves several issues that will require further elaboration upon in future work.Perhaps the most pressing issue is the extrapolation of data through Krylov iteration orders.At the moment we lack a clear physical understanding to motivate our extrapolation functions.We instead relied on the empirical quality of the fits and the consistency of the results with existing reference values from the literature.We do however discuss convergence in Krylov-like methods in Appendix A.
Our method easily generalizes to linear sigma models in two dimensions.In the unbroken sector, this would entail trivial modifications of the Feynman integrals by symmetry factors.From the analysis of the broken symmetry phase at strong coupling, it is presumably possible to probe nonlinear sigma models as well.
This method is easily extendable to higher dimensions and will serve as the next test of the method.For the ϕ 4 theory in 2 + 1d all that would be required would be to replace the form of the propagator.It may however be advantageous to implement the Feynman integrals in momentum space, as the real space propagator suffers from stronger singularities in higher dimensions.
The transition to momentum space could cause certain complications in the present setup, where the partial cancellation of disconnected pieces is implemented at the level of the integrand, leading to better temporal convergence.Strictly speaking the explicit projections leading to the extra disconnected pieces are not necessary.We might have omitted them at the cost of overlaps between the low-energy basis states and the tails.We opted for implementing the projectors in coordinate space as they act as a temporal cutoff, improving the convergence of matrix element integrals and decreasing their magnitude.Such cancellations in coordinate space do not obviously translate to momentum space.Fortunately, the structure of the disconnected pieces is known and can be taken into account by other means even if one opts for keeping them.The original idea of Krylov subspace methods is to compress the essence of an N × N dimensional operator A into an m ≪ N dimensional operator A (K) acting on K.The Krylov subspace K is spanned by the set {A k |r⟩ , 0 ≤ k < m}.A block Krylov subspace of block size n is obtained by the span of the set {A k |r i ⟩ , 0 ≤ k < m, 1 ≤ i ≤ n, n > 1} for linearly independent vectors r i , 1 ≤ i ≤ n.The solution to a linear system A |x⟩ = |b⟩ can be approximated by vectors |x K ⟩ + |δx (m) ⟩ for which A |δx (m) ⟩ ⊥ K. Eigenpairs (eigenvalue-eigenvector pairs) of A are approximated by eigenpairs of A (K) .
The case relevant to us is when the operator A is diagonalizable, i.e. related to a symmetric matrix Ā by a similarity transformation A = Q −1 ĀQ.If Ā is a bounded operator with a spectrum including at most a finite number ν of discrete negative eigenvalues, {λ − i } ν i=1 , the residual norm of the approximate solution after iteration m is bounded by the formula [94,95] where λ + max and λ + min are respectively the largest and smallest positive eigenvalues of Ā.
In our setup, the approximate tails are constructed iteratively from the linear system in eq. ( 18), and they are subsequently used as a variational basis for the determination of eigenvalues.Translating to the notation of the previous paragraph, in eq. ( 18) we would replace A with T ≡ 1 − T = ( Ĥ0,hh − E) −1 ( Ĥhh − E) and take |b⟩ ≡ |t l,0 ⟩ and |x⟩ ≡ |T l ⟩.One can then see the basis vectors in eq. ( 31) as linear combinations of vectors T i |t l,0 ⟩, 0 ≤ i ≤ k − 1. However unlike A above, in our case the operator T is infinite dimensional and not bounded.
Due to the unbounded character of T , It is difficult to make general statements about the convergence properties of the Krylov iteration (although recently progress has been made [96,97]).In particular, the bound (A1) formally gets ill-defined as there is no λ + max .In any case it is reassuring that even in this case the Krylov subspace only becomes exactly degenerate in a finite number of steps if the exact solution has been found.Indeed, if one can write . We will assume that E is strictly smaller than any eigenvalue of Ĥ0,hh .It follows that Q is positive.Moreover, when the parameter E is chosen to be the ground state energy in any symmetry sector respected by both Ĥ and Ĥ0 , then the symmetrized operator T = QT Q −1 is actually positive definite (in the given symmetry sector), as follows from the variational principle.
For concreteness, let us first consider the restricted version of eq. ( 18) in the presence of both a momentum cutoff Λ and a projection of the oscillator Hilbert space to a subspace spanned by Fock states with at most N particles (the zero mode is not affected by this cutoff).Then we have a restricted problem with a bounded operator T N,Λ (the truncated Hilbert space is finite dimensional).Two comments are in order.One is that in the limit Λ → ∞, the operator T N,Λ→∞ ≡ T N is, according our numerics, bounded, with a numerical behavior λ N,Λ max → λ N max + O(Λ −2 ).The second is that the limit T N is also strictly positive definite.The latter can be seen as a consequence of the lowenergy eigenvectors of H 0,ll having a finite weight in the ground state of H, while this ground state is separated from the rest of the spectrum by a gap.In the following we assume the existence of the above limit and suppress the Λ index.(In cases when the limit does not converge, it is still possible to reintroduce the cutoff.Then the statements below remain valid in the presence of a finite momentum cutoff.) The exact solution to the equation T N |T l,N ⟩ = |t l,0 ⟩ can be written as |T l,N ⟩ = |T l ⟩ + |δT l,N ⟩, where |δT l,N ⟩ is the difference between the regularized and the full (N → ∞, cutoff Λ) solution.It is safe to assume that ∥|δT l,N ⟩∥ → 0 for large N, as we expect to be able to approximate the eigenvectors arbitrarily well by systematically enlargening the truncated Hilbert space.
At iteration m, the vectors T k |t l,0 ⟩, 0 ≤ k ≤ m are identical to the vectors T k N =4(m+1) |t l,0 ⟩.This is because we build our tails on top of the oscillator vacuum, so |t l,0 ⟩ contains at most 4-particle states; on the other hand, the application of T changes the particle number by at most 4. Therefore, at iteration m, the approximate solution |T (m) l ⟩ can be written as l,4(m+1) ⟩.In turn, we can bound the error of the full iteration by looking at the regularized problem.Analogously to eq. (A1), we obtain the bound

∥T |δT
where λ min,m and λ max,m are the extremal eigenvalues of the operator T N =4m+4 with T N = QT N Q −1 .From the form of f (m) in eq.(A4) it follows that f (m) ≤ m.Depending on the large-m behavior of f (m), we can distinguish three cases: 1. f (m) → ∞: In this case the vanishing of the residual norm (with respect to T ) is guaranteed, with rate given by eq.(A4).This is because ∥Q −1 ∥ is bounded and ∥Q |t l,0 ⟩ ∥ = H (1) 1t,ll as defined in eq. ( 35).In particular, if f (m) ∼ αlog(m), this yields a power-law bound m −α .
(A4) yields a finite asymptotic bound, lim m→∞ ∥T δT 3. f (m) → 0: in this case the bound is essentially empty, stating simply that ∥T δT We note that the above bound provides a sufficient, but not necessary condition for convergence of the residual norm.To get an idea on the actual behavior of f (m) in the ϕ 4 model, we computed the extremal eigenvalues of T in a Fock space truncation for various momentum and particle number cutoffs.The results shown on Figure 18 indicate that we are consistent with case 1, at least numerically.Note that we calculate the eigenvalues for various particle number cutoffs N and use m = (N − 4)/4 and we include fractional values of m to increase the density of measurement points.
Many standard incarnations of Krylov subspace methods involve a step-by-step orthogonalization of the Krylov subspace.In higher orders, this orthonormalization procedure ensures numerical stability.In our setup, this orthogonalization is omitted.This is because we only consider low orders of the iteration and our matrix elements are provided with finite errors.As different matrix elements involve integrals of different dimensionality, the only way to perform a Gram-Schmidt-like procedure is to integrate first and propagate the error into the linear combinations, which is thus expected to increase.Since we deal with the generalized eigenvalue problem directly, in our case the explicit orthogonalization appears to be an unnecessary extra step.
We note in passing that there are alternative Krylov subspace methods that could be examined.In particular, we could have started from the original eigenvalue problem, and set up the Krylov subspace using powers of the full Hamiltonian Ĥ acting on eigenstates of the unperturbed part Ĥ0 .The eigenvalues could be approximated by a Lanczos method [98][99][100], that could be implemented in an approach very similar to ours.Note that the matrix elements of powers of Ĥ numerically grow significantly faster without the energy denominators inserted.The operator Ĥ is evidently unbounded in the Λ → ∞ limit, even if there is a projection to finite particle numbers.Alternatively, we could include the energy denominator as a preconditioner to the aforementioned Lanczos method.We do not cover these possibilities in the present work, but it may be useful to compare the performance of these alternatives on different models.
while in infinite volume, we have which leads to the relation where The explicit form of the Hamiltonian, expressed with fi-nite volume normal ordering, reads At this point the zero mode can be separated resulting in the following form where (B18)

Massless oscillators
When the massless oscillator basis is used, another renormal-ordering is necessary for the powers of φ: The conditionally convergent sum is understood in the presence of some cutoff |n| ≤ n max , and n max is subsequently taken to infinity.This results is a finite limit.The extra term on the RHS can be written in the following integral representation V α1α1 (T 3 , T 2 )V α2α2 (T 1 , 0) ⟨m| T g α1 (T 3 )g α1 (T 2 )g α2 (T 1 )g α2 (0)|m ′ ⟩ (C6) where T denotes time ordering so that larger time arguments are arranged to the left.The variables n 2 , n 3 , n 4 refer to the number of quadratic, cubic and quartic vertices excluding the one with zero argument, while Θ(T ) is the Heaviside unit step function.The compressed correlator V comp α1,...αn is obtainted from V α1,...αn in the following way.First we remove any disconnected Feynman diagrams appearing in V α1,...αn .Second, we distribute diagrams which can be transformed into each other by a permutation of the first n − 1 vertices into equivalence classes.We choose a representative diagram from each class so that α 1 ≤ α 2 ≤ • • • ≤ α n−1 , and multiply it with the number of diagrams in the class.Thus extra symmetry factors corresponding to this partial symmetrization are obtained by an explicit counting of equivalent diagrams.
The subset of Feynman diagrams containing only quartic vertices, for the lowest few orders are pictured in Figures 19 and 20.While restricting the evaluation to topologically different diagrams greatly reduces the com-plexity of evaluations, one has to keep track of the time arguments in the zero mode correlator, which have to be arranged "manually" into the correct time order in their numerical evaluation.
From the viewpoint of Feynman diagrams, the n-point functions contain disconnected pieces.(see e.g., Fig. 20 (a)).The explicit disconnected terms of I k1...kn n only partially cancel these terms, as the time ordering of the latter is partially fixed.Nonetheless, these terms combine in a nice way.Cancellation occurs whenever the disconnected correlator factors can be arranged into at least two subsets, so that the narrowest time interval encompassing all time arguments inside a single subset is disjoint from the corresponding intervals of the other subsets (Figure 21).In other words, there is a cutoff in the separation of disconnected pieces.Together with the exponential decay of the propagator, this feature ensures that the dominant contribution to the τ -integral comes from the the region with all τ variables being small.When the τ integrands are symmetrized under τ ordering, the regions where the disconnected terms give a contribution needs to be symmetrized explicitly.
The above method is suitable for calculating the ground state and zero-momentum one-particle states of the interacting theory.
Things become more complicated if one is to calculate general excited states.One could in principle include nontrivial oscillator states.When states having nontrivial oscillator content are included among H 1 states, one obtains Feynman diagrams with external lines, which adds to the complexity of calculating the associated matrix elements.

Choice of E *
All 'conventional' tail states depend on the energy parameter E * .Technically E * appears as the (a priori unknown) exact energy of the eigenvalue that we seek.In Fig. 1, the emphasis was on the numerical convergence in N ZM and so we took E * from reference data originally obtained with the method described in [1].In general it is tempting to perform a self-consistent iteration to fix this energy parameter.However, we argue that such an iteration is actually not favorable in general.The stage of the derivation in which E * is identified with the exact energy eigenvalue involves the full, non-truncated Hamiltonian.There is no reason why an iteration in a restricted Hilbert space should be optimal in obtaining the eigenvalue.A better point of view may be to consider E * as a variational parameter, and look for the optimal set of tails by minimizing the ground state energy with respect to E * .

Feynman diagrams
Practically, we generate the Feynman diagrams as follows.In our setup the relevant diagrams have no external legs, so we will restrict our attention to these "vacuum" diagrams.
1.For a given number of vertices n, generate all distinct ordered sets of vertex ranks (no. of legs) from which at least one diagram can be assembled.For the ϕ 4 model, 2 ≤ rank ≤ 4 for each vertex and the sum of all ranks has to be even.For example, there are six relevant ordered sets for n = 3. Henceforth we introduce the notation (r 1 , r 2 , . . ., r n ) for the set of diagrams having n vertices with ranks r 1 ,. . ., r n .
2. To each collection of vertices, e.g.(2, 2, 4, 4), construct a list of all possible Feynman diagrams by connecting the vertices in all possible ways, such that no line starts and ends on the same vertex.Disconnected diagrams are also kept.All possibilities are easily generated in Mathematica as long as the number of vertices is not too large.
3. Distribute the appropriate symmetry factors to the graphs.Each graph is multiplied by where R i is the number of legs of vertex i, and P ij is the number of propagators between vertices i and j.
4. If V comp is calculated, we select the connected Feynman diagrams appearing in V α1,...αn .We distribute diagrams which can be transformed into each other by a permutation of the first n − 1 vertices into equivalence classes.We choose a representative diagram from each class so that α 1 ≤ α 2 ≤ • • • ≤ α n−1 , and multiply it with the number of diagrams in the class.
5. Finally, the integrand corresponding to a given set of graphs is generated in real (Euclidean) space.

Disjoint interval
Overlapping intervals (b) An example of higher orders (arising in I 7 ): disconnected diagrams with disjoint argument configurations (like the one depicted) cancel.Here time increases from left to right (in particular, T 2 < T 1 .)The correlator ⟨V (T 6 )V (T 5 )⟩ is completely separated in its arguments from the other correlators which results in the cancellation by V (T 6 , T 5 )V (T 4 , T 3 , T 2 , T 1 , 0).

FIG. 21. Partial cancellation of disconnected terms
of diagrams approaches the order of 10 4 , the corresponding source code is in the order of megabytes, and it becomes convenient to divide the integrand into a number of separate C ++ functions and compile a static library.This greatly reduces compilation time.
and f3(Ecut) = a + bE −3 cut + cE −4 cut .We accept the result of the extrapolation with f4 as the central value, while the difference betweem fitting with f3 and f4 estimates the error.We opted for this to emphasize the uncertainty arising from the unknown exact form of the fitting function.Note that following the original error analysis in [1] yields an even sharper estimate.

4
, the normal ordering scales fixed to m 1 and m 2 , respectively, provided that the dimensionful quartic coupling agrees g

FIG. 1 .
FIG.1.The values of the ground state energy for the first three tail orders computed at g4 = 1.5, L = 10 for different choices of NZM .Blue, orange, green dots correspond to the first, second, and third Krylov orders, respectively.We see the results rapidly converge as a function of NZM .We compare our results to data provided to us by the authors of [1] (here shown as a dashed black line with error bounds shown in gray).
FIG. 7.Distribution of the lowest eigenvalues for 3 tail orders, g = 1.4 m0L = 10, based on an ensemble of Ns = 500 Hamiltonians, sampled assuming normal distributions for the matrix elements with standard deviations from errors derived from the Monte Carlo integration routines.The histogram is based on calculating the lowest 10 eigenvalues per parity sector for each matrix in the ensemble.Blue dots correspond to the even Z2 parity sector, while orange dots correspond to the odd sector.Energy bins of size ∆E = 0.005m 2 0 were used.

2 )FIG. 8 .
FIG.8.Bulk energy measurement in the broken sector from different volumes and compared to 4 loop L → ∞ perturbation theory (black).

FIG. 9 . 8 FIG. 10 .
FIG.9.Lowest energy gap measurement in the broken sector from different volumes and compared to 4 loop L → ∞ perturbation theory (black).At couplings g < 0.15 it corresponds to the B1 one-particle state, while for

FIG. 12 .
FIG.12.Fitted power law exponent to statistical ensemble of ground state eigenvalues Coupling g 4 = 0.4m 2 0 (M = 0.937, E = −0.00682).The corresponding result for the free boson with mass m 0 is shown for reference with a gray dashed line.

-
FIG.13.Volume dependence of the ground state energy in the unbroken sector, at Theoretical behavior including the leading Lüscher correction is shown with a continuous black curve .The results for the K = 1 and K = 2 universal tail basis computations are shown with blue and orange dots, respectively.
perturbative series truncated at various orders, compared to the Borel resummation of the perturbative series ([19]).Krylov TSM numerical data at m 0 L = 10, compared to[1]

FIG. 15 .
FIG.15.Mass gap in the unbroken sector, as function of coupling g.

2 0FIG. 16 .
FIG.16.Limitations of raw truncation (unbroken sector): squared norm of the projection of H |0⟩ to a truncated Hilbert space, normalized by its squared norm.Here |0⟩ is the unperturbed vacuum as a function of the TSM energy cutoff Etrunc, L = 10, g4 = 1.PT is a partial resolution of the identity in the truncated Hilbert space.Also shown are the squared norms of the projections to the universal tail spaces of 1st (blue) 2nd (orange) and 3rd (green, with light green uncertainty) orders.
Norm of the residual error for the ground state as the function of the coupling g in the unbroken phase.The results for the K = 1, 2 universal tails are shown with continuous blue and orange lines, respectively.Raw truncation with different energy cutoffs yield the set of dots at the top of the figure (Ecut = 5: light gray, Ecut = 7: gray, Ecut = 9: black).Norm of the residual error as function of the Fock TSM energy cutoff, for various values of g.
Appendix A: Relation to conventional Krylov subspace methods

FIG. 18 .
FIG.18.Numerical evaluation of the function f (m) of eq.(A4) for various momentum cutoffs Λ, at parameters g = 1, L = 10, E = −0.3941.We used Fock space truncation in the even parity sector in the presence of NZM = 10 kept zero mode states.

FIG. 19 .
FIG. 19.Lowest order Feynman diagrams contributing to the ground state energy: O(g 2 4 ) (top), O(g 3 4 ) (bottom).The vertex with the zero time argument is distinguished by being colored red.

4 (
Partial cancellation for the lowest order disconnected diagrams.Here T 3 ≥ T 2 ≥ T 1 ≥ 0. The top diagram is cancelled by the explicit subtraction in I ijkl {τ }), but the other two are not.

TABLE I .
Results for the critical point both in broken and unbroken phases from the literature.
116125964(91)g 4 + 0.3949534(18))g 5 −1.629794(22)g 6 + 7.85404(21)g 7 −43.02(21)g 8+ O(g 9 ) (48) forming H 1 , we use a finite dimensional subspace spanned by the first N ZM lowest energy eigenvectors of H ZM .H ZM is an interacting Hamiltonian with a discrete spectrum.The determination of the eigenstates |m⟩ is easily done numerically.To do so, we need to choose a basis in which to represent |m⟩.Here we choose a massive oscillator basis.Having specified H 1 , H 2 then consists of all states with non-trival oscillator content, any states involving a non-zero number of a † n .