A Cheap Alternative to the Lattice?

We show how to perform accurate, nonperturbative and controlled calculations in quantum field theory in d dimensions. We use the Truncated Conformal Space Approach (TCSA), a Hamiltonian method which exploits the conformal structure of the UV fixed point. The theory is regulated in the IR by putting it on a sphere of a large finite radius. The QFT Hamiltonian is expressed as a matrix in the Hilbert space of CFT states. After restricting ourselves to energies below a certain UV cutoff, an approximation to the spectrum is obtained by numerical diagonalization of the resulting finite-dimensional matrix. The cutoff dependence of the results can be computed and efficiently reduced via a renormalization procedure. We work out the details of the method for the phi^4 theory in d dimensions with d not necessarily integer. A numerical analysis is then performed for the specific case d = 2.5, a value chosen in the range where UV divergences are absent. By going from weak to intermediate to strong coupling, we are able to observe the symmetry-preserving, symmetry-breaking, and conformal phases of the theory, and perform rough measurements of masses and critical exponents. As a byproduct of our investigations we find that both the free and the interacting theories in non integral d are not unitary, which however does not seem to cause much effect at low energies.


Introduction
A renormalization group (RG) flow is fully specified by the UV fixed point where it starts, and by the perturbing relevant operator which gives the initial direction. Once this UV data is given, everything else about the flow must follow, including the IR properties, and should be computable. Is there a general algorithm suitable for performing such a computation? We have in mind in particular strongly coupled situations, when perturbation theory is not applicable.
For strongly coupled flows starting at a free theory, IR physics is accessible by Monte Carlo simulations on the lattice. What about flows which start at a strongly interacting conformal field theory (CFT)? One could try to put this CFT on the lattice first, as a critical point of some lattice model, or as a fixed point of another RG flow starting at a higher scale from a free theory. However, this is baroque, and may not always be possible. Moreover, this dodges the question. Nonperturbatively, the UV CFT is specified by the conformal data (the spectrum of its local operators and the structure constants of the operator algebra). The relevant perturbation is specified by the operators by which we are perturbing and their relative coefficients. As a matter of principle, all IR physics must be computable in terms of only this nonperturbative UV data.
The purpose of this paper is to thrust into the limelight one method designed to solve this problem-the Truncated Conformal Space Approach (TCSA). While the method is familiar in the statistical mechanics and condensed matter community, so far it has not had much attention from high energy physicists. The field of potential applications of TCSA is huge and poorly explored. In particular, all prior work has been done in d = 2 spacetime dimensions, and one of the goals of this paper will be to provide a generalization to d > 2.
Although the original purpose of the TCSA was to study flows starting at an interacting CFT, the method is perfectly applicable when the UV CFT is free and therefore represents an alternative to the lattice for studying such flows. In TCSA, statistical errors are absent and systematic errors are very different from the lattice Monte Carlo. Including fermions (both Dirac and chiral) is straightforward; including gauge fields is more complicated but seems doable.
Time will tell if TCSA is a viable alternative to the lattice for flows in more than two dimensions. To get an initial feeling, in this paper we study the Landau-Ginzburg theory-the free scalar theory in d > 2 dimensions perturbed by two relevant operators φ 2 and φ 4 . Depending on the relative value of the quartic coupling and the mass, this theory flows in the IR to a massive theory with the Z 2 -symmetry preserved or spontaneously broken. For a critical value of the coupling the theory flows to an interacting IR CFT in the universality class of the d-dimensional Ising model. Using TCSA, we are able to observe this phase structure, and to obtain reasonably precise predictions about the finite-volume spectrum of the theory, at weak, intermediate, or strong coupling. The required computational resources turn out to be quite minor-the total cost of all the computations in this paper is O(10) single-core days on a desktop. A bonus is that the theory for any d (including fractional d) can be studied within the same unified framework. All in all, our initial experience with the TCSA has been quite positive, and justifies further explorations of the method.
The paper is organized as follows. We start in section 2 with a general discussion of TCSA and of its connection to the Rayleigh-Ritz method in quantum mechanics. Section 3 is devoted to the free massless scalar in d dimensions quantized on the sphere. We focus on the CFT description of the spectrum of theory, via radial quantization, but we also discuss the relation to more conventional canonical quantization. As a byproduct of this discussion, we show that the free massless scalar in a fractional number of dimensions is actually a non-unitary theory-its Hilbert space contains negative norm states. This curious fact does not seem to have been noticed before.
We continue in section 4 with more details about TCSA for perturbations of the free massless scalar. In particular, we describe an OPE-based method which allows to efficiently compute the matrix entering the TCSA eigenvalue problem. This method should be useful beyond the Landau-Ginzburg flows studied in this paper.
Section 5 is central to the paper-it explains how TCSA results depend on the cutoff energy of the Hilbert space, and how this dependence can be reduced by applying a renormalization procedure. The results of this section are crucial for improving the accuracy of TCSA.
In sections 6 and 7 we apply TCSA to study the Landau-Ginzburg flow. We first discuss in section 6 the case when only the mass parameter is nonzero, i.e. we study the free massive theory as a perturbation of the free massless scalar by the mass term. This trivial theory allows us to test TCSA and the renormalization procedure in a controlled situation where the exact answer is known. After this test, we proceed in section 7 to study the general case when both the mass and the quartic are turned on. For reasons which will be explained below, we perform this study in d = 2.5 dimensions. In spite of this exotic spacetime dimensionality, we find all the traits expected from the Landau-Ginzburg flow for d = 3. In particular, depending on the values of the couplings, we observe both phases of the theory, separated by a continuous phase transition. We perform rough measurements of the mass spectrum of the theory in both phases, and of the critical exponents at the phase transition. We also comment about how the non-unitarity of the theory at fractional d manifests itself via some high-energy eigenvalues acquiring imaginary parts.
We conclude in section 8 with a discussion of future research directions. Appendices A and B contain a brief review of prior work, on TCSA in d = 2, and on other Hamiltonian truncation techniques in quantum field theory. Two more appendices contain derivations of auxiliary technical results.

Truncated Conformal Space Approach
We would like to study an RG flow obtained by perturbing a d-dimensional CFT by a relevant scalar operator V. 1 To set up the calculation, we will first of all need an IR regulator. The most natural IR regulator for a CFT is to put the theory on the "cylinder" R × S d−1 R , where the radius of the sphere R serves as an IR scale. The map to the cylinder amounts to Weyl-transforming the metric of the flat Euclidean space. The CFT dilatation generator then maps to the Hamiltonian on the cylinder, and the CFT local operators O i of dimensions ∆ i map to states |i on the cylinder whose energies are given by In the theories we will be considering here, there will be a unique ground state corresponding to the unit operator, with energy zero. 2 The Hamiltonian of the perturbed theory on the cylinder is The key idea is to think about this Hamiltonian as an infinite matrix in the Hilbert space of unperturbed CFT states |i . The CFT Hamiltonian in this basis is diagonal and simply related to the CFT operator dimensions: 3 i|H CFT |j = R −1 ∆ j δ ij . (2. 3) The matrix elements of the perturbation on the other hand are related to the CFT three point functions and will provide an off-diagonal piece of the Hamiltonian. Schematically, we will have is the coefficient of the CFT three point function of the unit-normalized local operators in flat space (it is thus R independent). Precise prefactors and normalizations will be discussed below. Notice that the dependence on R follows by dimensional analysis. The mass scale of the flow is given by (2.5) For R Λ −1 IR we are close to the UV regime, where V is a small correction to H CFT and perturbation theory is reliable. To probe the IR physics, we must take instead Here V cannot be treated as a small perturbation, and the right thing to do would be to diagonalize the whole Hamiltonian H CFT + V . But how can we do this given that this matrix is infinite?
The trick is to introduce a UV cutoff Λ UV and to truncate the Hilbert space keeping only the states below this maximal energy: If the cutoff is chosen so that Λ UV Λ IR , (2.8) we can hope that the IR physics is not much affected. Let us furthermore assume that the UV CFT has a discrete spectrum, which will be true for most CFTs of interest. In this case, the truncated Hilbert space is finite-dimensional, H is a finite matrix and can be numerically diagonalized. One then learns a wealth of information about the theory in the IR. For example: • the ground state dependence on R gives the vacuum energy density; • by looking at the number of exponentially degenerate ground states we can infer the symmetry breaking pattern; • the excited states give the massive spectrum of the theory, including one-particle, manyparticle, and bound states; • studying the dependence of two-particle states energies on R we can extract the S-matrix; • for flows ending in conformal fixed points we can extract the spectrum of IR operator dimensions.
This, then, is the essence of the TCSA, first proposed in 1991 by Yurov and Al. Zamolodchikov [1]. In practice, one tries to take Λ UV as high as possible, but this will be limited by the rapid growth of the number of states with energy. The success of the method depends on whether we can get reasonable results with numerically tractable Hilbert space sizes.
Both the original paper [1] and all the subsequent TCSA literature known to us consider d = 2, but here we presented directly the general d case because the basic logic is very similar. Our focus in this paper will be on d > 2. The current status of the TCSA research in d = 2 is summarized in Appendix A, to which the reader will be referred from time to time.

A Digression: Truncation in Quantum Mechanics
When one first hears about the TCSA, the usual reaction is incredulity. How can such a naive method solve such a hard problem? To demystify it a bit, recall that very similar techniques are routinely used in quantum mechanics under the name of the Rayleigh-Ritz method. For example, consider the problem of finding the spectrum of the anharmonic oscillator: Let us express this Hamiltonian as an infinite matrix in the Hilbert space spanned by the eigenstates |n of the harmonic part H 0 . Then truncate the basis by keeping only the states |n below some cutoff, n N . The claim is that in the limit N → ∞ the eigenvalues of H diagonalized within the truncated Hilbert space tend to the exact anharmonic oscillator energy levels, as obtained e.g. by solving the Schrödinger equation. This exercise is extremely easy to carry out (we recommend it to the reader!), since the matrix elements of q 4 ∝ (a + a † ) 4 are known in closed form. One also finds that in this example the convergence is exponentially fast. 4 The method works equally well for both small and large λ. 5 As is well known, perturbation theory would be divergent for the anharmonic oscillator problem for any value of the quartic, and requires Borel resummation. However, the Rayleigh-Ritz method is completely immune to this. In fact, Rayleigh-Ritz is probably the most practical method to find energy levels of the anharmonic oscillator. Double well potentials can also be considered, by including a negative frequency term to the perturbation: (2.10) Rayleigh-Ritz is used everywhere in quantum mechanics where perturbation theory is insufficient. As a less trivial example, consider the helium atom Hamiltonian: where H 1 and H 2 describe the electrons in the Coulomb field of the nucleus, and V 12 is the electron-electron interaction. Here one could e.g. work in the Hilbert space spanned by the tensor products of the H 1 and H 2 eigenstates, compute the matrix elements of V 12 , truncate and diagonalize. In practice, it helps to choose a different basis, which does not diagonalize H 1 + H 2 , but takes into account that the two-electron wavefunction will be singular at the coincident points (this improves the convergence rate). This is how the helium energy levels were obtained with a precision adequate to test QED quantum corrections [2].
In fact, Rayleigh-Ritz is nothing but a version of the time-honored variational method. The variational method is usually used for the ground state, and one is often content to get a reasonable estimate for its energy. With Rayleigh-Ritz, one usually aims at a precision determination, and gets the excited states as well as the ground state (convergence being fastest for the low-lying states). The variational method approximates the ground state energy from above. By the minimax principle, the same is true for all the truncated Rayleigh-Ritz eigenvalues: (2.12) In quantum mechanics, truncation methods can often be put on solid ground by establishing rigorous results about its convergence. Simplest results of this type can be found e.g. in [3], section XIII.2, and many more are scattered through the literature. As in the helium example above, the convergence rate is influenced by how well singularities of the exact wavefunction in the coordinate representation are reproduced by the functions of the trial basis. For quantum field theory we are not aware of general mathematical results about the convergence of truncation methods, but there is a large body of evidence that TCSA does converge. The evidence in d = 2 is summarized in Appendix A, and examples in d > 2 will be studied in this paper.

A Case Study for TCSA in d Dimensions
In this paper we will apply TCSA to study the Landau-Ginzburg theory, i.e. the free massless scalar theory perturbed by a linear combination of : φ 2 : and : φ 4 : operators. This is perhaps the simplest d-dimensional flow. A priori, we are interested in 2 d < 4. However, in this first work we will have to stay away from the extremes of this range, since the TCSA analysis becomes complicated near these extremes.
The reason why d close to 2 is hard for TCSA is that the scalar dimension ∆ φ = (d − 2)/2 approaches zero in this limit, and the free scalar spectrum becomes dense and eventually continuous in d = 2. 6 To have a sufficiently sparse spectrum, we will keep d not too close to 2.
The reason why d close to 4 is hard is that φ 4 becomes marginal in this limit. On the contrary, as we will see, TCSA works best for strongly relevant perturbing operators V. The more relevant the operator is, the better-behaved perturbation problem is in the UV. The best situation is realized when which for the perturbations considered here means (2.14) When (2.13) is satisfied, the perturbation is simply UV-finite. At ∆ V = d/2 the vacuum energy becomes divergent, as can be seen at second order in perturbation theory. Other UV divergences appear if we further increase ∆ V , and these also affect the couplings of nontrivial local operators (including V itself). These short-distance divergences have to be handled in the usual QFT wayby adding counterterms. In this first work we would like to avoid dealing with UV divergences, so we will stay within the bounds (2.14). This does not mean however that we will altogether ignore cutoff dependence. Even in the range (2.14) when there are no UV divergences, the accuracy of the method for a finite cutoff will be influenced by power-suppressed corrections. This important issue will be discussed below.
As the reader must have noticed, we are considering the case of fractional d on equal footing with the physically interesting integer d. We will see that the TCSA problem allows a natural continuation to general d.

Free Scalar in d Dimensions
In this section we will discuss the UV CFT at which our RG flows will be starting-the free massless scalar CFT in d dimensions. These results presented here lay the groundwork for the numerical investigations and for the renormalization, studied in the subsequent sections.

Scalar Operators
The local operators of the free boson theory are built by taking products of the fundamental field φ and of its derivatives, e.g.
: φ 2 ∂φ ∂∂φ ∂∂∂φ : , where some or all of the vector indices on the derivatives may be contracted. The operators are all inserted at the same point, and the normal-ordering sign means as usual that we don't consider Wick contractions within the operator when computing its correlation functions with other operators.
We can classify the operators according to their spin, i.e. their representation under SO(d). When we put the theory on the cylinder, the spin of an operator becomes the spin of the state into which it maps under the state-operator correspondence. Eventually we will perturb the theory by adding to the Hamiltonian an integral of a scalar operator over the sphere, as in Eq. (2.2). Since this perturbation preserves rotation symmetry of the sphere, the Hamiltonian matrix will split into blocks corresponding to the states of the same spin. The scalar sector contains most of the states we are interested in: the ground state, one-particle states at rest, and two-particle states in the center-of-mass frame. In the large R limit, many of the states of higher spin will correspond to spin 0 states slightly boosted along the sphere. 7 In principle, there could also exist additional states with intrinsic spin, which could be thought of as bound states of fundamental scalars at strong coupling. This would be analogous to vector mesons in gauge theories with matter. In this paper we however focus exclusively on the scalar sector.
Operators in the scalar sector will look like (3.1) with all the indices contracted. In (3.1) there is only one way to contract the indices to get a scalar operator, but in general there may be several inequivalent ways. It is useful to encode a scalar operator by a graph where φ with n derivatives corresponds to a vertex with n edges sticking out of it, and contracting two derivatives means joining two vertices with an edge. Notice that we can ignore operators containing contractions of derivatives acting on the same φ, since ∂ 2 φ = 0 by the equation of motion. On the other hand, two vertices may be connected by more than one edge, i.e. our graphs can have parallel edges. In graph theory, graphs obeying these conditions are called multigraphs without loops. 8 Depending on how derivatives are contracted, the graphs may have one or several connected components. In particular, each φ without derivatives will give an isolated vertex, see figure 1. It's also clear that isomorphic graphs correspond to identical operators, and should not be counted separately.

Parity-Odd and Null States
The above discussion was incomplete in two ways. First, we could also consider contractions involving the -tensor, which correspond to the parity-odd operators, as opposed to the parity-even operators discussed above. These two classes of operators won't mix under the parity-preserving φ 2 and φ 4 perturbations we will consider. In this paper we focus on the P = +1 sector. Notice that the P = −1 scalar operators have pretty high dimensions, e.g. the lowest dimension one in Consequently, the P = −1 states in the IR are also likely to be heavier than for P = +1.
Second, it's important to realize that for integer d different graphs sometimes give rise to the same operator. The simplest example in d = 3 is This can also be also written as where we consider M µν = φ ,µν as a symmetric 3 × 3 matrix, traceless by the equations of motion. Then (3.4) is just a statement about symmetric polynomials built out of its eigenvalues.
As one goes higher in energy, one encounters infinitely many relations of this type. For example, for d = 3 all graphs tr M n , n 4, are expressible via linear combinations of the products of tr M 2 and tr M 3 . There are also many other similar relations. One then has two alternative ways to proceed. Either one may decide to find all such relations below the UV cutoff one is working at, and explicitly eliminate all graphs which are not independent. Alternatively, one may decide not to eliminate anything, and work in an extended Hilbert space containing one state per graph. In this Hilbert space, Eq. (3.4) would be interpreted as a null state condition, i.e. it means that a certain linear combination of states has zero overlap with any other state (and in particular zero norm).
In this paper we will follow the second approach, because we would like to treat integer and fractional d within the same formalism. However, null state relations only exist for integer and finite d. At fractional d, all non-isomorphic graphs give rise to inequivalent states. In future work focussing on integer d, it will probably make sense to eliminate the null states, in order to reduce the dimension of the Hilbert space and speed up the subsequent matrix diagonalization.

Non-Unitarity at Fractional d
The above discussion of null states raises an interesting question-what is the precise fate of the null states when one passes from integer d to a nearby fractional d? As we mentioned, these states are then no longer null, but are they positive-or negative-norm? We claim that some of these states acquire a negative norm.
To see a concrete example, let us take the operator O in the LHS of Eq. (3.4). By an explicit computation, its two point function for a general d is given by Notice that the zero at d = 3 is first order. For 2 < d < 3 the two point function of O is negative, and so O must have an overlap with a negative-norm state. This example can be easily generalized to show that there are negative norm states for any fractional d, in fact infinitely many of them. 9 The presence of negative-norm states means that the free scalar theory in fractional d is not unitary. To our knowledge, this observation has not been made before, although theories in fractional dimensions have been extensively studied, especially in relation to critical phenomena, where they form the basis of the -expansion. As we will see below in section 7.4, the lack of unitarity will lead to the presence of complex energy eigenvalues, once the free theory is perturbed by the quartic coupling. However, the mass term alone leaves all energy levels real-see section 3.5.
It has to be said that the first negative norm state has a pretty high dimension: This must be the reason why they have not been noticed until now. A few negative norm states at high dimensions, hidden among lots of positive-norm states of comparable dimensions, probably do not have a strong effect on the low-energy physics. In a recent conformal bootstrap study of the Wilson-Fisher fixed point in fractional dimensions [7] it was assumed that these theories were unitary, and very reasonable results were obtained. 10

Primaries and Descendants
CFT local operators can be divided into primaries and descendants. For d = 2 descendants are obtained by acting on the primaries with the raising part of the Virasoro algebra. For general d considered here, descendants are simply derivatives of the primaries.
The basis for the Hilbert space on the cylinder includes of course all states, those corresponding to primaries and to descendants. If the UV CFT is strongly coupled, the matrix elements (2.4) between primary states are part of the non-perturbative conformal data, while matrix elements involving descendant states can be computed using the conformal algebra. This situation often occurs in 2d TCSA applications (see Appendix A).
On the other hand, in our case the UV CFT is free and we will follow an alternative approach. Instead of relating the descendant matrix elements to those of the primaries, we will simply evaluate all the necessary matrix elements using the fact that our operators are built out of the fundamental scalar field. This procedure is much faster, because the additional work required to classify the states into primaries and descendants is significant. In particular, many scalar states will be descendants of primaries with spin, whereas our approach completely avoids introducing states with spin.
Let us finally remark that the current mathematical understanding of the underlying algebraic structure for CFTs in fractional d appears to be rather incomplete. This is particularly relevant for the primary/descendant classification of local operators, which would require defining so(d) algebras and their representations in fractional d. 11 Fortunately we need not be concerned with this in practice. All computations involving scalar operators, which have all indices contracted, can be done directly for fractional d by setting ∆ φ = (d − 2)/2 and δ µ µ = d.

Relation to Canonical Quantization
So far, our discussion of the free massless scalar on the cylinder has been based entirely on the radial quantization and the state-operator correspondence. Here we would like to comment on a more low-brow approach-canonical quantization. The free scalar in curved space is described by the action As is well known, the theory is Weyl-invariant for m 2 = 0 and the coupling to the Ricci scalar ξ = d−2 4(d−1) . We now quantize canonically on the cylinder metric The field φ is expanded in the eigenfunctions of the Laplacian on the sphere of radius R (spherical functions): φ l,n , l = 0, 1, 2 . . . , (3.9) where l is the angular momentum quantum number, and n numbers the states in the multiplet. The Laplacian eigenvalue is l(l + d − 2)/R 2 . To each of these modes we will associate a harmonic oscillator of frequency where we used the fact that Ric = (d − 1)(d − 2)/R 2 for a round d − 1 dimensional sphere of radius R. The Hilbert space is the Fock space of these oscillators.
The free massless scalar is recovered setting m → 0. In this limit, ω l reduces to (l + ∆ φ )/R, which is the right energy for the state corresponding in the radial quantization to the operator obtained by acting on φ by l derivatives.
What are the relative disadvantages and merits of the canonical vs radial quantization description of the free scalar Hilbert space? One definite merit of canonical quantization is that it allows to take the mass into account nonperturbatively from the start. On the contrary, in the radial quantization approach we have to perturb around m 2 = 0. In section 6 below we will do the exercise of reproducing the free massive scalar spectrum within this framework.
On the other hand, an advantage of perturbing around the conformal point m 2 = 0 is that the matrix elements of the perturbation have a simple overall power-law dependence on the radius R, see Eq. (2.4). So we don't have to recompute the matrix when we change the radius. On the contrary, in canonical quantization with a general mass, matrix elements will depend on R nontrivially via the 1/ √ ω l normalization factors accompanying the oscillators. Thus the matrix will have to be computed separately for every R.
In this paper, we will perturb around the conformal point and use exclusively TCSA (i.e. radial quantization). In the future, it would be interesting to see if canonical quantization with m 2 = 0 gives better results. In upcoming papers [10,11], we will use canonical quantization to study the Landau-Ginzburg flows in d = 2. As explained in section 2.2, TCSA is not directly applicable to these flows (see however footnote 6).

The Gram Matrix
As made clear in the above discussion, we will be working in the Hilbert space of scalar states on the cylinder, in the basis which in the radial quantization can be identified with scalar operators acting on the vacuum: As already mentioned in footnote 3, this basis in general will not be orthonormal and not even orthogonal. Rather, we will have a nontrivial Gram matrix. This Gram matrix is an essential ingredient in the existing implementations of the TCSA in d = 2 (see Appendix A).
In our case, the Gram matrix will not play a crucial role. In fact, as we will see below, the perturbed spectrum computation can be organized without using the Gram matrix at all. Nevertheless, the Gram matrix is a conceptually important object, so we would like to discuss in some detail its definition and evaluation.
First we need a map from the states to their conjugates. As usual in radially quantized CFT, this map is defined with the help of the inversion transformation R : where the conjugate operator [O i (x)] † is inserted at the point Rx. The rules for construction of the conjugate operators are as follows (φ is the fundamental scalar; A, B any two fields in the theory): Starting from Rule 1 and using Rule 2 repeatedly we can conjugate all derivatives of φ. Then by applying Rule 3 we can conjugate all normal-ordered products of derivatives, and in particular all scalar operators forming our basis.
Computation of the Gram matrix is thus reduced to evaluating two-point functions of operators made of several φ's acted upon by various derivatives. In principle, this is straightforward to do using Wick's theorem. The number of Wick contractions to perform can be dramatically reduced by using selection rules. To begin with, the only nonzero entries are those for which This rule can be easily understood using the relation with canonical quantization described in the previous section. In this description, N l maps to the total occupation number of oscillators with angular momentum l. To get nonzero overlap, all angular momentum modes should have the same occupation number.
Putting it all together, the Gram matrix is thus evaluated as follows. First one computes the overlaps ∂ l {µ} φ|∂ l {ν} φ , (3.14) by using the above prescription, or by using the conformal algebra, as explained e.g. in [12]. These are particular invariant tensors, symmetric and traceless in both groups of indices {µ}, {ν}.
Overlaps between general scalar states are then computed by contracting the basic overlaps (3.14) between their constituents.
Using this direct algorithm, we could compute the Gram matrix up to a rather high cutoff in operator dimension. However, we found it expensive to compute in this way the Gram matrix all the way up to the cutoffs we will be using in the numerical analyses below. An alternative, indirect, method for computing the Gram matrix will be described in section 4.3. That method is much faster and easily yields the Gram matrix up to the required values of the cutoff. In any case, as we will see below, the spectrum computations can be organized avoiding the use of the Gram matrix. In this paper, the Gram matrix has been used only in one instance-to count the number of negative norm states plotted in figure 20.

Simple versus Generalized Eigenvalue Problem
Let us formalize a bit more our problem. Energy levels on the cylinder are solutions of the eigenvalue problem H|ψ = E|ψ . (4.1) We will be looking for scalar eigenstates, expanding them in a basis of states |j : The states |j will be in one-to one correspondence with the scalar local operators of the UV CFT (in this paper, the free massless scalar theory). The Hamiltonian in this basis will be represented by a matrix: In terms of this matrix, Eq. (4.1) becomes a simple eigenvalue problem Notice that the matrix H i j is not hermitean. To transform the problem to a hermitean form, we consider the matrix elements We of course have where G ik = i|k is the Gram matrix discussed above. The matrix G ij is hermitean, and H ij is hermitean if the Hamiltonian is hermitean, which is the case for the Landau-Ginzburg flows with real couplings considered here. Actually, for the operator bases considered in this work, these matrices will be real symmetric. We then have an equivalent symmetric generalized eigenvalue problem: In the existing d = 2 TCSA implementations (see Appendix A), one starts by computing the matrices G ij and H ij , which naturally leads to the generalized eigenvalue problem (4.7). One then usually multiplies both sides by G −1 and transform to (4.4). 13 Here, we will choose an alternative path. Namely, we will directly compute the matrix H i j and find eigenvalues from (4.4). The method for computing H i j is described below in section 4.3.

Working in Presence of Null States
As mentioned in section 3.2, we will be working in a basis which, for integer d, will contain null states. In presence of null states the above discussion needs to be reconsidered. In particular, the Hamiltonian matrix is then ambiguous since we can add an arbitrary null state to the RHS in (4.3): Also, the eigenvalue problem (4.1) has to be considered modulo addition of an arbitrary null state in the RHS. In practice, however, we won't have to deal with these subtleties. We will compute the Hamiltonian matrix as if there were no null states, 14 and solve the original eigenvalue problem (4.1). Our final spectrum for integer d will thus contain both physical and null state eigenvalues. It's easy to see that the physical state eigenvalues are the same as in the more rigorous treatment. 15 The null eigenvalues are unphysical-they have to be separated and thrown out. There are many ways to do this in practice: one can follow a null eigenvalue from the UV where its value is known; one can detect it by the presence of crossings with physical states (physical eigenvalues don't cross in RG flows which are not integrable); one can check the nullness of the corresponding eigenvector.
For the low-lying spectrum this issue does not even arise, since the first null state has a pretty high dimension.

Matrix Element Evaluation: OPE Method
The CFT piece of the Hamiltonian matrix is diagonal: The nontrivial part is to compute the matrix of the perturbation. We compute this by using the following OPE method. Namely, in the radial quantization we are supposed to compute where V is the perturbing operator (in our examples it will be : φ 2 : or : φ 4 :), and O j is the operator corresponding to the state |j . Consider the OPE where A k (0) are local operators inserted at the origin, and C k (x) are c-number coefficient functions.
Since we are in a free theory, this OPE can be worked out explicitly. Notice that while V and O j will be scalars, many of the operators A k will be tensors, and {µ} stands collectively for their indices, contracted with those of C k . Now to evaluate (4.10) we just integrate the OPE term by term, which amounts to integrating the coefficients: (4.12) 14 We perform this computation in Mathematica keeping d as a free parameter, and set d to the desired value before the diagonalization. 15 A key to this argument is that null states can only be mapped into null states by the Hamiltonian. In principle, the fact that we don't solve (4.3) modulo the appearance of a null state could lead to some physical eigenvectors disappearing, due to the Jordan block phenomenon. However, this is very non-generic and would be easily detectable as we are varying parameters such as couplings and the radius of the cylinder. We have never observed it happen. By rotation invariance, the integrals will produce invariant tensors, i.e. a number of Kronecker deltas connecting the indices in {µ}. Contracting these with the indices of A {µ} k will give scalar operators. Expressing the RHS of (4.12) in the original basis, we read off the matrix V i j . In the above discussion we were effectively setting R and g to unity. To restore the dependence in these parameters, we need to multiply the resulting matrix by This, then, is how we compute the matrix entering the eigenvalue problem. Notice that this approach is more economical (involves fewer Wick contractions) than the direct computation of the three-point functions i|H|j .
The matrices V i j computed by the OPE method can be subjected to a check. We know that if we multiply them by the Gram matrix as in (4.6), the resulting matrix V ij must be symmetric. As mentioned at the end of section 3.6, we can compute the Gram matrix directly up to a rather high cutoff. Up to this cutoff, we can then check the symmetry of V ij for the φ 2 and φ 4 perturbationsthis check works.
For still higher cutoff, we found it expensive to compute the Gram matrix directly. However, we can ask the following question: given our matrices V i j , is there a symmetric matrix G ij , subject to the selection rule (3.13), and with the property that V ij = G ik V k j is symmetric, for both the φ 2 and the φ 4 perturbations? It turns out that up to the highest cutoffs explored in this work, such a matrix always exists and, moreover, is unique! This provides an alternative, indirect, method to compute the Gram matrix. It is by this method that we computed the Gram matrix used to count the negative norm states in figure 20.

General Remarks
In the following sections we will see in concrete examples that TCSA converges as Λ UV is increased. We will also see that the rate of convergence is power-like. We would like to examine here why the method converges, and how its convergence rate can be improved. We thus have to understand the effect of removing the high energy states from the Hilbert space on the low energy spectrum. Effective field theory intuition tells us that this effect should be small, and can be corrected for.
The basic equations are as follows. We work in the Hilbert space of the unperturbed CFT on the cylinder, which is divided into the low (l) and high (h) energy parts: where H l includes all states of energy up to Λ UV . The full Hamiltonian is a block matrix: where H ab maps H b into H a . The TCSA truncated Hamiltonian is the upper left corner: The full eigenvalue problem is Let us now eliminate c h by using the second equation. We get: This exact equation should be compared to the truncated equation used in TCSA: Here, we writeĒ,c rather than E, c to indicate that these are solutions to the truncated equation rather than Eq. (5.3).
Conclusion: TCSA will converge if the matrix correction in (5.5) can be neglected in the limit Λ UV → ∞. Naively, this seems likely since it is suppressed by H hh − E, and we are assuming that E belongs to the low-energy spectrum, while the eigenvalues of H hh will be presumably large.
However, the precise statement will depend also on the size of the matrix elements mixing H h into H l . This mixing being due to the perturbation, we can expect that the importance of corrections will depend on ∆ V .
Let us view the problem from a practical angle. Suppose we know an eigenvalueĒ and the corresponding eigenvectorc of the truncated problem (5.6). How can we correctĒ to get closer to the solution of the exact eigenvalue equation? Let us write the full Hamiltonian as We took into account that the off-diagonal elements H hl and H lh are associated only with the perturbation V . The eigenvalues of H 0 are known-these are the TCSA eigenvalues and the unperturbed eigenvalues of the diagonal H CFT,h . We will now view H 1 as a perturbation and compute corrections to the TCSA eigenvalues. By the usual Rayleigh-Schrödinger perturbation theory we get: Further corrections terms are simple to write down. For example, the next one is given by We will only use the term shown in (5.9) in this paper, but in the future increasing the accuracy of the renormalization procedure will likely require mastering (5.10) and perhaps even further terms.
Our job is not yet finished, since evaluating the correction term (5.9) requires an infinite summation over the states in H h . It would be desirable to find a simplified approximate form for this correction: where V c act simply on H l . For example, V c might be of the same form as V itself, i.e. an integral of a local operator V c over the sphere. Then adding ∆H to the TCSA Hamiltonian can be thought of as renormalizing the couplings. In the language of usual perturbative quantum field theory, the V c might be called counterterms. The difference is that in perturbation theory, we usually worry only about the counterterms which diverge when the UV cutoff is taken to infinity. Here we care also about the correction terms which are power-suppressed-we will want to add them in order to improve the accuracy of the method. The Hamiltonian with added correction terms can be called an improved TCSA Hamiltonian. This is analogous to "improved actions" in Lattice QCD.

Computation of ∆H
To find the correction terms, we examine the matrix element of ∆H between two states i, j ∈ H l : where E n = ∆ n /R stands for the unperturbed CFT energy. We will estimate the large energy asymptotics of M n . The key idea is to consider the correlation function: Inserting the resolution of unity, this correlation function can be represented through the same M n as: 14) The large energy behavior of M n can then be extracted from the part of C(τ ) which is non-analytic as τ → 0, since the low energy states give rise to an analytic contribution.
A moment's thought shows that nonanalyticity for τ → 0 can appear only from the region where the nonintegrated correlator has a singularity, i.e. from n close to n . In this region we can use the OPE To the accuracy needed below, it will be sufficient to use only the shown leading term in the OPE. Moreover, we will be considering only scalars in the RHS of the OPE. With a Poincaré-invariant cutoff, non-scalar operators are not induced in the renormalization group flow. However, the TCSA regulator is more subtle. We break the Poincaré group to SO(d) times dilatations. Furthermore, since we are working in a Hamiltonian formalism, we may find integrals of tensorial operators induced by the RG flow. As an example, the appearance of the stress tensor T µν on the RHS gives (after integrating over the sphere) a contribution of the form so it leads to a renormalization of the coefficient of H CFT in the TCSA Hamiltonian. 17 However, since the stress tensor and other operators with spin have high dimension, their effects will be suppressed compared to the effects of the scalars by a higher power of Λ UV .
Each term in the OPE will give rise to a term in the τ → 0 asymptotics of the correlator. The prefactor will be given by the matrix element of V c integrated over the sphere, while the dependence on τ will come from the integral of the OPE kernel. Up to O(τ 2 ) accuracy we have (see Appendix C): .
This non-analytic behavior can be reproduced provided that the large-dimension distribution of the coefficients M n contains a component with a power law: It should be kept in mind that M n is a discrete sequence, and so the given continuous distribution is supposed to approximate it only on average. Below we will discuss the accuracy of this approximation in more detail. Also, for the renormalization of the φ 2 flow we will work out the asymptotics of the sequence M n via an alternative method.
For the moment, to get an expression for ∆H, we introduce the shown asymptotics into (5.12) and perform the sum approximating it by an integral. Gathering all the prefactors, reinstating the dependence on the coupling constant and on R, we obtain the following formula for the correction term: For very large Λ UV we are allowed to drop the corrections due toĒ and ∆ i +∆ j in the denominator. However, below we will find it useful to keep track of these subleading corrections, at least approximately.

Renormalization Group Improvement
In the above discussion we were assuming that ∆H is very small, and correcting eigenvalues by the leading-order perturbation formula (5.8) is adequate. For this, Λ UV has to be taken sufficiently large so that the renormalizations of all the couplings implied by (5.20) are small compared to their values in the bare TCSA Hamiltonian. This condition is rather restrictive and in fact in our main example below-the Landau-Ginzburg flow-we will not be able to satisfy it, as the mass renormalization due to the quartic will sometimes be comparable to the bare mass.
The way out in such a situation is to perform an RG improvement of the correction procedure. This is inspired by the usual RG in perturbative quantum field theory, which resums large logarithms. Here we don't have logarithms but power-suppressed terms, but the logic is very similar.
As usual in RG, we imagine performing a sequence of Hilbert space reductions with cutoffs Λ 1 > Λ 2 > . . . To first order in ∆H, Eq. (5.8) is equivalent to diagonalizing H + ∆H. So instead of using (5.8) we will just keep adding ∆H to the original Hamiltonian, and iterate. We can imagine changing the cutoff in each step from Λ to Λ − δΛ for δΛ Λ. Concretely, using the form of ∆H given in (5.20) and assuming that the corrections due to ∆ i + ∆ j andĒ can be ignored, this procedure results in RG equations of the form: In this way we obtain a flow in the space of Hamiltonians, which we can integrate all the way down to the desired cutoff Λ UV . It may be expected that, under certain circumstances, the final 'resummed' Hamiltonian obtained by such a procedure will have a larger range of applicability (i.e. work for smaller Λ UV ) than the first-order correction formula. This will be the case if the subleading on the right-hand side of (5.9), such as (5.10), are less important than the terms we are proposing to resum. Again, we can draw analogy from the usual perturbative RG, when the beta function is the sum of one-loop, two-loop etc terms. Our RG improvement is like integrating the one-loop beta function while dropping all higher-loop terms.
As already mentioned, in the examples considered below we will want to keep track of the corrections due to (∆ i + ∆ j )/R andĒ in (5.20). These corrections are state-dependent, and taking them into account completely would require a separate RG flow for every value of these parameters-a complication that we wish to avoid. Instead, we would like to find a practical way to represent them by operators. The easiest way to do so is to expand in powers of the inverse cutoff and keep only the first-order terms. In that case we can replace ∆/R by H CFT andĒ by H. For example, in the case where V c = 1, the first of these two subleading corrections can be thought of as 'wave function' renormalization of the coefficient of H CFT in the TCSA Hamitonian, 18 while the second correction becomes a uniform overall rescaling of all couplings. Both of these can be taken into account easily by a slight modifications of (5.21). The situation is more complicated if V c = 1. In this case, the expansion generates terms of the form 19 Not only are these terms not present in the original Hamiltonian, they are also of a qualitatively different type-they are not given as an integral of a local operator over the sphere. In other words, these terms are nonlocal. While this may seem confusing, a moment's thought shows that this was to be expected. The reason is that the TCSA regulator-throwing out all states above a certain energy-is not a fully local UV regulator. 20 So we have to learn to live with nonlocal correction 18 This shows once again that corrections due to integrals of non-scalar operators, T τ τ in this case, can be induced by the flow with the TCSA cutoff. In the previous subsection, we pointed out that the correction due to the direct appearance of T µν in the OPE would be suppressed, since the corresponding coefficient h is quite large. However, here we are discovering another way for the appearance of this correction-as a subleading term accompanying the unit operator in the OPE. 19 Notice that these corrections, as written, preserve the hermiticity of the Hamiltonian. 20 The TCSA regulator reproduces exact correlators as long as the insertion points are separated in the time direction by Λ UV −1 . In particular, correlation functions on a constant time slice are not faithfully reproduced no matter how far the points are separated in the space direction. By a fully local UV regulator we mean a regulator which reproduces exact correlation functions as long as points are separated in some direction, time or space, by Λ UV −1 . E.g. the point splitting procedure, used in conformal perturbation theory, is a fully local regulator.
terms. Fortunately, from the practical point of view the terms (5.22) pose no problem. First of all, they are easily computable, since they are given by products of matrices which we anyway have to compute in the earlier stages of the TCSA procedure. Secondly, although in principle the non-local terms would appear also on the right-hand side of the RG flow equations, which would substantially complicate the flow, in practice we found that they remain rather small compared to the local terms. This happens because their running is suppressed by one extra power of the cutoff. Therefore, in this work we will ignore backreaction of the non-local terms on the other running couplings.
Although the above procedure correctly takes into account the leadingĒ/Λ dependence, we realized that expanding inĒ/Λ is actually not a reasonable thing to do at large R. The point is that the ground state energy E 0 grows at large R like R d−1 , and even for moderately large R becomes non-negligible compared to Λ UV . Whether this is a problem depends on the sign of E 0 . If E 0 were to become large and positive, there would be no magic way out-the correction procedure would break down as soon as E 0 ∼ Λ UV , as seen e.g. by the blow up of the integral in (5.20). Fortunately, the ground state energy density at large R is usually negative. 21 In this case, although E 0 becomes large in absolute value, nothing bad occurs with the correction in (5.20); it even decreases with respect to the E 0 = 0 case. However, were one to expand inĒ/Λ, one would unnecessarily introduce large corrections even in this benign case.
We will therefore adopt the following prescription. We will replace the estimateĒ in (5.20) by E r + (Ē − E r ) where E r is a convenient reference energy that we estimate to be close to the expected value ofĒ. For example, we may choose E r to be around the ground state energy as obtained by extrapolation from lower values of the radius, or around the energy of the first excited state. In fact, the end results for the spectrum should not depend much on the chosen value of E r , which provides a consistency check for the method. We then expand not inĒ/Λ but instead in the difference (Ē − E r )/Λ, which is not expected to become large in the large volume limit. The RG evolution is then performed keeping track of the exact dependence on E r (no expansion) through the simple substitution Λ d−h+1 → Λ d−h (Λ − E r ) in the denominator of (5.21). Since we will expand in (Ē − E r )/Λ, the leading correction in (5.22) should be modified by replacing H → (H − E r ). Below we will see a concrete example of how this works, when discussing the Landau-Ginzburg flow.

Comments on Earlier Treatments of Renormalization
Cutoff dependence and renormalization have been discussed in the context of the d = 2 TCSA studies, most importantly in [13] (following [14,15]) but see also Appendix A for other references. In particular, Section 3 of [13] discusses in detail how the cutoff dependence can be analyzed using the OPE, and gives renormalization group equations similar to our (5.21) for the couplings of the local operators. At leading order, then, their results are basically equivalent to ours. 22 Ref. [13] also initiated a discussion of subleading terms. For example, the first of the two subleading terms in (5.22) may be discerned in their equations, for the special case where V c is the identity operator. However, significant differences do exist between us and them at how these subleading effects are implemented.
According to the prescription in Section 4.2 of [13], on top of leading RG improvement, each IR state should get a subleading correction factor computed from the conformal perturbation theory applied to a UV state from which the IR state in question originates. This prescription, as well as a more recent detailed discussion in Section 3 of [16], are designed to fix up, order by order in the coupling, the discrepancies between TCSA and conformal perturbation theory. On the other hand, our discussion uses from the very beginning the fact that the true expansion parameter is the inverse cutoff rather than the couplings, which become large in the IR.
Let's illustrate the differences by looking at the correction in Eq. (5.8). Our derivation demonstrates clearly that one should compute ∆H with the nonperturbative energyĒ and take the matrix element between the nonperturbative statesc. A similar correction in Eq. (3.11) of [16] uses the UV energy and the UV state in place ofĒ andc. At small R the two methods would give very similar results, but at large R the difference will be significant. Indeed, the IR states at large R will have a complicated composition, in which the original UV state carries little weight (see e.g. figure 7). Also the energy in the IR will get a very large correction, implying a large change in the denominator in (5.9). As a result the whole correction may be modified at O(1). Out of curiosity, we compared our method to that of [16] for the φ 2 flow discussed in the next section. We found that at large R our method is more effective in reducing the discrepancy from the exact results.
For completeness it should be noted that [16], using their renormalization prescription, achieved an excellent agreement of TCSA data with the results obtained by exact integrability methods applicable for the model they studied. This success is puzzling to us, since as we explained we believe that their prescription is problematic at large R. This question deserves further analysis. 6 The φ 2 Flow

Theoretical Expectations
We are now ready to do our first TCSA calculation in d dimensions-the flow starting at the free massless scalar and perturbing by the pure mass term 1 2 m 2 : φ 2 :. Needless to say, this RG flow is considered for illustrative purposes only. We expect to find the free massive scalar theory, whose spectrum in canonical quantization was discussed in section 3.5. We will restrict ourselves to the spin 0 states of that spectrum, corresponding to one or several particles at rest or in relative motion along the sphere in such a way that the total angular momentum adds up to zero.
In addition to the spectrum, another observable is the ground state energy. For the canonically 22 A factor 1 2 seems to be missing in their Eq. (3.7). Even having corrected this misprint, we did not manage to reproduce their figure 1(c). quantized massive scalar it's given by the zero point energy of all oscillators where D d (l) is the size of the spin l symmetric traceless representation of SO(d) (see e.g. [17]): Here we will be computing the ground state energy by perturbing a CFT, and so our expected answer is not E 0,can but a simple modification thereof. First of all, in our treatment the unperturbed CFT ground state energy on the sphere is set to zero (see footnote 2). Furthermore, the O(m 2 ) term in the ground state energy must vanish, being proportional to a CFT one-point function which is zero. We expect however that all terms higher-order in m 2 will agree between CFT and canonical quantization. Thus, the perturbed CFT ground state energy should be given by (6.1) with terms of the zeroth and first order in m 2 dropped. 23 This can be written as follows: where µ = mR and Because of subtractions, the general term in the series behaves at large l as l d−5 . So the ground state energy is finite for d < 4, in agreement with the criterion (2.13), (2.14). 24 For general R, we can compute E 0 by numerically summing the series (6.3). In the large volume limit µ 1, the sum will be dominated by large l terms and can be approximated by an integral. The leading behavior in this limit scales as the volume of the sphere, with a constant density set by the mass: One can show that the first correction in this formula arises at order 1/R 2 . The same is true for the rate of approach of masses of particle states to their infinite volume limit, as can be seen from the explicit formulas in section 3.5. The presence of these power-like in R corrections is due to the curvature of the general d-dimensional sphere. Here we are observing them in a free theory, and we expect them to be present in an interacting situation as well. This can be contrasted with what happens when a QFT is put on a torus T d−1 (which for d = 2 is of course the same as S d−1 ). In this case it's been observed long ago [18] that masses in an interacting theory are affected by terms which are exponentially small in the size of the torus.

TCSA Setup
Which value of d shall we choose in our numerical study of the φ 2 flow? As already mentioned in section 2.2, the case d → 2 is expected to be difficult, as the CFT spectrum is becoming dense in the limit. For d > 4 the vacuum energy will be divergent. Here we will show results for the physical value d = 3. We have also performed checks for other nearby values of d and they work equally well.
We will construct the truncated Hilbert space by including all scalar operators below a certain dimension For practical reasons, this maximal dimension will be held fixed when varying R. This means that we will be working with a sliding UV cutoff TCSA can be expected to reproduce the IR spectrum roughly below this cutoff. As we increase R, the sliding cutoff decreases and eventually becomes comparable with m. At this point TCSA results can no longer be trusted. 25 The number N 0 (∆) of scalar, parity even states in the Hilbert space as a function of ∆ is shown in figure 2. The dotteed line gives the number of physical states, counted using group theory. The total number of states in a d-dimensional CFT grows with ∆ exponentially [17] 26 : where C is a theory dependent constant related to the prefactor in the free energy density dependence on the temperature. 27 It's not hard to see that the number of scalar states will also grow exponentially with the same exponent, although with a smaller prefactor. In figure 2 we can clearly see this exponential growth. The fast growth of the number of states implies that it will be hard to increase ∆ max . The success of any TCSA calculation will depend on whether reasonable results can be obtained with a manageable ∆ max . As we will see, achieving numerical accuracy for such ∆ max will require the use of renormalization corrections discussed in section 5.
As discussed in section 3.2, in this paper we will be working in an extended Hilbert space which for integer d is somewhat larger than the physical Hilbert space, since it includes some null states. The extended Hilbert space states are in one-to-one correspondence with non-isomorphic multigraphs; their number is shown by blue squares, while red dots show null states. We see that the first null state occurs at ∆ = 10; this is the state (3.4).

Numerical Results
We will now show our numerical TCSA results and compare them with theoretical expectations. The TCSA computation starts by constructing the truncated Hilbert space. We consider cutoffs up 25 The range of validity of TCSA will be somewhat extended due to the fact that the induced ground state energy is negative, which reduces the size of correction terms, as we discussed in section 5  to ∆ max = 18, which corresponds to 4573 scalar P -even states. We then construct the Hamiltonian matrix H i j . The CFT part is given by the diagonal matrix (4.9). The perturbing part is computed using the OPE method from section 4.3, for the operator V = : φ 2 :. We then diagonalize the Hamiltonian matrix to find the spectrum of the perturbed theory.
Since the perturbation preserves the Z 2 symmetry which maps φ → −φ, the Hamiltonian matrix does not mix Z 2 -even and Z 2 -odd states. The two sectors have roughly equal number of states, and it makes sense to do the computation separately in each of them, reducing the size of the matrices to be diagonalized by factor ∼ 2. The ground state belongs to the Z 2 -even sector.
In figure 3, we plotted E 0 as a function of R. In this and other plots in this section, we set m = 1, which means that we measure R in units of m −1 and energies in units of m. The black solid curve shows the theoretical prediction for E 0 (R) obtained by summing the series in Eq. (6.3). Let us focus first on the 'raw' TCSA results, i.e. the results obtained without any renormalization corrections (blue curves marked 'raw'). We see that the agreement is good up to R ∼ 1, while for larger R there are noticeable deviations. As the cutoff is increased (we show ∆ max = 12 (18) in dashed(solid) blue), the numerical results are moving towards the theoretical prediction, but the convergence is not very fast.
Cutoff dependence of TCSA predictions was discussed in section 5. As we have seen, the errors induced by omitting the states with E > Λ UV are expected to go down as a series of power-laws with known exponents. These errors can be understood analytically and then subtracted away, greatly improving the accuracy of the TCSA calculations. The red curves marked 'ren.' in figure 3 have been produced using such a renormalization procedure (see section 6.4 for details). The agreement with the exact results is greatly improved; it now extends up to R ∼ 2.5. Notice that the corrected results also exhibit a smaller dependence on the cutoff. This is because we are subtracting the leading correction, and the remaining ones are suppressed by extra powers of Λ UV .
Is R ∼ 2.5 a large or a small radius? For two reasons, it should be considered as large. First of all, the corresponding sphere circumference, L = 2πR, is much larger than the inverse mass, so that there is plenty of room for a massive particle wavefunction to fit into the sphere. Second, it takes us much beyond the radius of convergence of conformal perturbation theory R c = (d−2)/2. We now turn to the excitations above the vacuum. In figure 4 we plot the energies of these excitations, subtracting the vacuum energy. We have two plots, one for the Z 2 -even and one for the Z 2 -odd sectors. To keep the plots from cluttering, we show the lowest five eigenvalues in each sector. Notice that in both cases we subtract the same quantity E 0 , which is the lowest energy in the Z 2 -even sector. Blue dots are computed using TCSA for ∆ max = 18, while lines joining them are added to guide the eye.
In the same plot thin magenta lines show the exact free massive scalar spectrum, computed by combining oscillator energy levels from section 3.5. In the Z 2 -even sector the lowest state corresponds to two particles at rest, while the states above it correspond to two particles with some angular momentum on the sphere combined in a state of total spin zero. Then there comes the state with four particles at rest etc. In the Z 2 -odd sector we recognize one particle at rest, three particles at rest, then three-particle states in relative motion, etc.
In figure 5 we show the same but for the spectra computed using the renormalized TCSA eigenvalues. The details of the renormalization procedure will be discussed in the next section. We see from these plots that renormalization extends the range of R where TCSA is in agreement with the exact results from R 2 to R 3.

Renormalization Details
The general method of renormalization was presented in section 5. Here we will describe particular issues which arise when the procedure is applied to the φ 2 flow. The leading contributions to the   correction term ∆H are expected to come from the low-dimension operators in the φ 2 × φ 2 OPE: 2N d |x| d−2 : φ 2 : + : φ 4 : + . . . , (6.9) Here N d is the normalization factor in the two point function of the canonically normalized massless scalar: where S d = 2π d/2 /Γ(d/2) is the area of the unit sphere in d dimensions.
Now, curiously, although the operators φ 2 and φ 4 appear in the OPE (6.9), their contributions to the renormalization corrections vanish. Indeed, the coefficient B(h) given by (5.17) is zero for the corresponding h's. The reasons this happens are not difficult to understand; they are ultimately related to the fact that the UV CFT we are perturbing is free. Starting with φ 4 , notice that since the dimensions factorize, ∆(φ 4 ) = 2∆(φ 2 ), and the OPE kernel is just a constant. Clearly, the τ → 0 limit discussed in section 5.2 is perfectly analytic in this case, and so B(h) must vanish. For φ 2 , although the OPE kernel is singular, it is a harmonic function of x − y. By the mean value property, the integral of a harmonic function over a sphere is equal to its value at the center of the sphere. This implies that also in this case the τ → 0 limit is analytic, and B(h) = 0. 29 Thus the only leading non-vanishing correction is for V c = 1. 30 We will have to include this correction taking into account the subleading dependence on ∆ i + ∆ j andĒ. Indeed, were we to drop these subleading parts, we would get a constant counterterm which would shift all eigenvalues in the same way. This would have a chance to improve the agreement for the ground state energy, but would have no effect on the spectrum of massive excitations. However, the raw TCSA massive spectra in figure 4 do show noticeable deviations, which we would also like to improve.
In fact, we will be able to do even better. Not only will we include the above-mentioned subleading effects, but we will also take into account the discreteness of the sequence M n . Recall that the general formula (5.18) gives this sequence only on average. However, it turns out that for the φ 2 flow the tail of the M n sequence can be worked out explicitly, independently of the argument in section 5.2. As we show in appendix D, M n at ∆ n ∆ j is nonzero only if ∆ n −∆ j −∆(φ 2 ) = 2p is an even integer, in which case it's given by It's not difficult to see that on average this sequence does agree with the continuous distribution (5.18), which also provides a check for the general argument.
We next evaluate ∆H via Eq. (5.12). When doing the sum, we use the expression (6.11) for all terms. This is not quite true, since (6.11) was derived under the assumption ∆ n ∆ j , which does not hold for the external states i, j just below and n just above the cutoff. However, the caused error cannot be large, since the states i, j close to the cutoff will anyway contribute little to the renormalization, having small weight in the eigenvectorc; see figure 7 below. So we will tolerate this little imprecision. The infinite sum over p becomes a 4 F 3 hypergeometric sum, and specializing to d = 3 we obtain the digamma function ψ(z). Reinstating the coupling and radius dependence, we get: where K j is defined as the smallest odd integer such that ∆ j + K j > ∆ max .
The leading term in ∆H for large ∆ max is a state-independent correction ∝ ∆ −1 max . As mentioned above, keeping only this correction would not be adequate. Instead, we will use the full expression (6.12) to compute corrected ('renormalized') eigenvalues E ren from the raw TCSA eigenvaluesĒ via the formula (5.8): It is these 'renormalized' results which were used to produce figures 3,5. Herec is the eigenvector corresponding to the raw TCSA eigenvalueĒ. In this approach, each energy level is corrected separately.
Note that to apply formula (6.13) we need to compute both right eigenvectorc j , as in (4.4), and the left eigenvectorc i : 14) The eigenvectors are assumed normalized viac ic i = 1. Of course these two eigenvectors are related, up to normalization, via the Gram matrix:c i ∝ G ijc j . (6.15) As mentioned in section 3.6, computing the full Gram matrix may be expensive, although we did find an indirect way to do it, described in section 4.3. If one has access to the Gram matrix, one can use it to compute the left eigenvectors via (6.15). Without the Gram matrix, one simply finds c i from the second eigenvalue problem in (6.14). 31 Figures 3 and 5 do demonstrate that our renormalization procedure works-upon applying the renormalization corrections, discrepancy from the exact results is reduced compared to the raw TCSA data. Figure 6 demonstrates the same as a function of the UV cutoff: we show how the TCSA ground state energy and the massive spectrum converge to their exact values with the gradual increase of ∆ max . We do this plot for one value R = 2, but the picture is qualitatively the same for all R. This figure shows that not only the accuracy is greatly improved after the renormalization, but the convergence rate is also improved. This is because the error terms remaining after the leading renormalization subtractions are suppressed by higher powers of 1/Λ UV .
One last aspect we would like to discuss here is an assumption implicit in the entire procedure of renormalization, namely that the contribution of high energy states to low energy observables 31 It should be noted that the nonsymmetric eigenvalue problems are somewhat more difficult to solve numerically than the symmetric ones, and more prone to numerical instabilities. In our work we overcome the instabilities by applying the transformation H → (H −σ) −1 to the matrix H before diagonalization. This transformation focuses on the eigenvalues nearest to σ. We also checked some of our results by working at a higher number of digits. In future work, it would be interesting to keep looking for other, more numerically efficient diagonalization procedures.  In figure 7, we plot this distribution for the lowest Z 2 -even massive excitation (the one which is interpreted as a state of two particles at rest) and for several values of R. As expected, for small R the distribution is strongly peaked at ∆ = ∆(φ 2 ) = 1. As R is increased, the distribution becomes wider and wider, but its high-energy tail does remain suppressed. The same qualitative behavior is true for the other states. One can wonder what it would mean if for very large R the distribution becomes flat or even peaked at high ∆. Does this ever happen for CFTs perturbed by a relevant operator? Presumably the method would completely break down for such R, but for the values of R explored in this work this does not happen. 7 The Landau-Ginzburg Flow

Theoretical Expectations
In the free massive flow considered in the previous section, we could compare TCSA results with the exact theoretical answers for all observables. In a significant fraction of the d = 2 TCSA literature (see Appendix A), the method is applied to integrable flows, where the exact results are also available (of course, through much harder work than for the free massive scalar). We would however like to encourage the use of TCSA in situations when no other technique is readily available. This will be the case for most physically interesting strongly coupled theories, since exact integrability is possible only in d = 2, and even then it's an exception, not a rule.
In this spirit, we will now use TCSA to study the Landau-Ginzburg flow-the flow obtained by perturbing the free massless scalar CFT by The quartic λ should be positive to have a stable vacuum, while m 2 can be positive or negative. The IR physics of this theory depends on the dimensionless ratio The theory is not integrable even in d = 2, and here we will study it in d > 2.
The case of small quartic coupling corresponds to |t| 1. In this regime the theory in the IR describes weakly interacting massive particles, and predictions can be obtained from perturbation theory. For positive (and still large) t, the perturbative vacuum is at φ = 0, and the Z 2 symmetry φ → −φ is preserved. On the other hand, for negative t, perturbation theory is developed around one of two degenerate vacua of the double-well potential. Thus, the Z 2 symmetry is spontaneously broken.
It is then interesting to know what happens for t = O(1), when the IR theory is strongly coupled, and perturbation theory is not useful. 32 One generally expects that the Z 2 broken and preserving phases extend into the strongly coupled region, where they are separated by a secondorder phase transition at t = t c , see figure 8. At t = t c the theory is expected to flow in the IR to a CFT, belonging to the Wilson-Fisher family of fixed points in the Ising model universality class.  In this paper we will only study the m 2 > 0 (i.e. t > 0) part of the phase diagram. Instead of varying t as in figure 8, we will find it convenient to work in the units m = 1, and vary λ. Using TCSA, we will compute how the finite volume spectrum of the theory depends on λ. As we will see, for small λ the spectrum will be consistent with preserved Z 2 symmetry, while for λ > λ c our calculations will indicate spontaneous Z 2 symmetry breaking. Thus we will obtain qualitative confirmation of the phase diagram in figure 8, and quantitative information about the massive spectrum in the strongly coupled region. We will determine the critical value λ c with some precision. For λ = λ c we will observe the mass gap going to zero, indicating that the IR theory is conformal. We will be able to get a rough estimate of the leading critical exponents at the phase transition point and compare them with the known values in the Ising universality class.
Notice that since our calculations indicate a positive value of t c = 1/λ (4−d)/2 c , the whole region t < 0 is expected to be in the Z 2 -broken phase. However, we have not explored this region numerically.

Numerical Results
We will now perform TCSA analysis of the Landau-Ginzburg flow in d = 2.5 dimensions, working with cutoff up to ∆ max = 17, which corresponds to 5494 (4907) Z 2 -even (odd) states. As already mentioned above, we will set m 2 = 1. The spectrum will depend on λ and the TCSA radius R. We will explore the region R 3 and λ 1.15 (see figure 9). Raw TCSA without renormalization corrections would give converged predictions only in the lower left corner of this region, corresponding to weak coupling and small physical volume. To extend the range of applicability of the method, we will apply a renormalization procedure as described in section 5, with the theory specific details described in section 7.3 below. To reduce the number of plots, we will only show results with all renormalization corrections taken into account.

Spectrum for a fixed R and varying λ
To visualize the spectrum dependence, we will plot it along a number of vertical and horizontal sections in the two-dimensional range in figure 9. Let us start with plots at a fixed R and varying λ. In figure 10 we show the results for R = 2.5. The ground state energy E 0 is defined as the lowest energy in the Z 2 even sector. The excitation spectra are given by E i − E 0 , in the Z 2 -odd and the Z 2 -even sectors, respectively.
We see that as λ is increased, the excitation energies first decrease, and then, for λ > λ c ≈ 0.5 -0.6, start increasing again. An interesting feature of the spectrum at λ > λ c is an approximate double degeneracy of states in the Z 2 -even and odd sectors, well visible for the vacuum and the first couple of excited levels. This behavior is the telltale sign that the theory for λ > λ c is in the phase of spontaneously broken Z 2 symmetry. We then expect that the theory at λ = λ c is conformal. This expectation will be further tested below.
It may be somewhat counterintuitive that the Z 2 symmetry breaks for a positive value of m 2 (remember that we fixed m = 1). In fact, there is no paradox. The m 2 is a UV parameter defining the initial direction of the flow, while the breaking is an IR phenomenon. As we flow from UV to IR, m 2 is renormalized and the effective squared mass may become negative. 33 In other words, we may imagine that a double-well potential is generated by quantum effects. In this case there will be two degenerate vacua and all excitations above the vacua should be degenerate as well. The degeneracy would be exact in infinite volume. In finite volume we expect some mixing due to the potential barrier tunneling, 34 so that the exact eigenstates are Z 2 -even or Z 2 -odd and split by a small amount (exponentially small for large volume). The mixing and the splitting are expected to become more important for higher energy states, for which the tunneling is not suppressed. All these intuitive expectations are confirmed by figure 10.  Another interesting feature of the spectra in figure 10 is the absence of level crossingeigenstates belonging to the same Z 2 sector don't cross. There are several values of λ when a pair of same Z 2 -parity eigenstates come close to each other, but then repel. This should be 33 See section 7.3 for the RG equations for the Landau-Ginzburg flow. 34 Such tunneling effects were for example studied in TCSA in Ref. [21]. contrasted with the free massive flow spectra, which do show level crossings, reproduced by TCSA calculations. The difference stems from the fact that the φ 2 flow is integrable, while the Landau-Ginzburg flow is not. This way of distinguishing integrable and non-integrable flows has long been noticed in the d = 2 TCSA literature (see e.g. [22,23]), and here we are observing it in d > 2.

Mass gap as a function of λ and determination of λ c
We will now further test the expectation that the theory at λ = λ c is conformal. In figure 11 we plot the low-lying spectrum of excitations (just the first three states) for λ varying from 0 to 1.15 and for three values of R = 2, 2.5, 3. We see that the excitation energies are decreasing with R for λ near λ c . This is especially noticeable for the second and third excited level. Away from λ c the spectrum is relatively stable with R.  The decrease of the spectrum with R at λ = λ c is what one should expect if the critical theory is conformal. Indeed, for a flow ending in a conformal fixed point the excitation energies should behave at large R as ∆ IR i /R where ∆ IR i are the IR CFT operator dimensions. We will test this expectation in the next section.
Let us now determine the critical value of the coupling λ c with some precision. According to the standard renormalization theory, the mass gap for λ near λ c should depend on λ as where ν is a critical exponent, 36 which in the case at hand is related to the dimension of -the first Z 2 -even scalar operator at the Wilson-Fisher fixed point: 35 Or even slightly increasing. We observed that this slight increase of the spectrum with R is reduced when raising the cutoff, so it must be attributed to truncation effects. 36 This critical exponent ν should not be confused with the shorthand ν ≡ (d − 2)/2 that was introduced in Eq. (3.10).
In our case the mass gap is E 1 − E 0 for λ < λ c and E 2 − E 0 for λ > λ c . We will fit E 1 − E 0 in the region λ < λ c to determine λ c and ν. We exclude the region λ > 0.5 from the fit since it is clearly affected by finite R effects which smear out the expected power-law behavior. We also exclude the region λ < 0.3 since Eq. (7.3) is expected to be valid only in the λ → λ c limit. We thus perform the fit in an interval [λ 1 , λ 2 ], and to estimate the systematic uncertainty we vary λ 1 between 0.3 and 0.4 and λ 2 between 0.45 and 0.5. This gives the following rough estimates for the critical coupling and the exponent ν (see figure 12): with a positive correlation between λ c and ν. We will now compare our determination of ν with the results by other approaches. The dimension ∆ for d = 2.5 dimensions can be extracted from the Borel-resummed epsilon-expansion series [24]. It can also be determined from the conformal bootstrap under the assumption that the Wilson-Fisher fixed point lives at a kink in the region of the (∆ σ , ∆ ) plane [7]. The latter analysis was done under the assumption that the Wilson-Fisher fixed point in fractional dimensions is unitary, which as we now know is not quite true. However, as noticed in section 3.3, a small fraction of high-dimension negative-norm states should not have strong influence on the conformal bootstrap predictions. This probably explains why [7] found no disagreement with the results of [24]. Both analyses predict: which gives a value ν ≈ 0.755, close to the upper end of the confidence interval for ν determined by our fitting procedure above. Assuming this precise value of ν and repeating the fits leads to a somewhat more accurate determination of the critical coupling:

Spectrum for a fixed λ and varying R
We next present how the spectrum depends on R for a fixed value of λ ( figure 13). We pick three representative values of the quartic: λ = 0.3 for the Z 2 -preserving phase, λ = 0.55 near the presumed critical point, and λ = 0.9 in the Z 2 -breaking phase. We will now comment upon what we see in those plots, first for the ground state energy, and then for the massive spectra.

Ground state energy
The ground state energy is expected to grow for large R as a constant times R d−1 , corresponding to a finite energy density (cosmological constant) induced by the RG flow. This behavior is clearly visible in the data for λ = 0.3, 0.55, while for λ = 0.9 the fit is not so good and there are significant deviations for R 2. These deviations decrease with ∆ max and are thus due to truncation effects. Jumping a bit ahead, notice that there are no comparably flagrant deviations in the massive excitation spectrum for λ = 0.9 and large R. This is because the largest truncation effects are expected in the coefficient of the unit operator, which has the smallest possible dimension (0), and the unit operator affects the ground state energy but not the spectrum.
Excitations for λ = 0.3 Since the energies are observed to tend to finite nonzero limits for R → ∞, we conclude that the IR theory is massive. The lightest Z 2 -odd state E 1 is a scalar particle of mass The next two excitations, belonging to the Z 2 -even sector, are readily interpreted as two-particle states. The former, of mass ≈ 2M , must have both particles at rest, while in the latter the particles must be in relative motion with respect to each other, with total angular momentum zero. Higher up, we observe a state of three particles at rest, of mass ≈ 3M , and orbital excitations thereof.
The appearance of this hierarchy of states, quantized in units of the lowest excitation, is a nontrivial consistency test on the results. It is also a prediction for the absence of bound states. At weak coupling, the two-particle interaction is repulsive in the Z 2 -symmetric phase of the φ 4 theory, so we wouldn't expect bound states. Our results show that this conclusion remains valid at strong(er) coupling. Notice that the physical mass M is significantly less than the bare mass m, so that the theory we are examining is presumably moderately to strongly coupled.
It is interesting to study the rate with which the excitation energies approach their infinite volume limits. Focussing first on the lowest excitation, the leading correction is expected to arise from coupling to curvature and scale as 1/R 2 : where A is a (theory-dependent) constant. Indeed, when the theory is put in a weakly curved background, it should be possible to describe corrections to the lightest state energy by an effective Lagrangian of the same form as the free massive scalar Lagrangian (3.7) with m → M and an effective ξ which will, in general, be different from ξ Weyl in the UV. Then we get (7.9) with A = ξ/ξ Weyl . 37 In figure 14 we test Eq. (7.9) for the lowest excitation of the λ = 0.3 spectrum. We see that it describes the large R approach reasonably well up to R ∼ 1.5, where the truncation effects apparently kick in and make the error to increase rather than decrease with R. Fitting the correction in the range R = 0.4 -1.5 we determine M ≈ 0.57, A ≈ 1.05.
One could wonder why A is so close to one in the case at hand. As already mentioned, we don't expect that A should be universal. We will encounter A < 0 below in the Z 2 -broken phase. Moreover, the coupling to curvature will be suppressed if the state in question is a pseudo-Goldstone boson, as may happen for more complicated flows with a continuous global symmetry. So, for the pion in QCD we expect A ∼ (m π /Λ QCD ) 2 1.
We next discuss the rate of approach for the two-particle states, starting with the two-particle state at rest. The energy of this state can be approximated as E 2 − E 0 = 2(M + ∆M curv ) + ∆M scat , (7.10) 37 Notice that for d = 2 the curvature vanishes, and the modification of the mass spectrum is entirely due to boundary conditions. The leading correction in this case is exponentially small [18]: where the last correction is due to the interaction (scattering) between two particles put into a finite volume. Since the interaction is short-range, we expect the leading correction of this type to scale as the inverse volume of the box [25]. In figure 15 we plot E 2 − E 0 − 2M for the λ = 0.3 spectrum. We see that the difference is not well described by the finite-volume correction 2 ∆M curv alone. A much better agreement can be obtained including a correction with the scaling ∝ 1/R d−1 , as would be expected from a scattering correction. Notice that the sign of the so determined scattering correction is positive, corresponding to a repulsion between the constituent particles. Indeed, as we already mentioned above, at weak coupling in the unbroken phase, λφ 4 interaction is repulsive; here we see the same effect persisting at strong coupling. In principle, it should be possible to relate the size of the scattering correction to the scattering length, as was done for a flat torus by Lüscher [25]. In would be interesting to work out the corresponding theory for the sphere. Finally, in figure 16 we plot the difference E i − E 0 − 2M for the lowest two orbital excitations in the two-particle sector, which should consist of two particles moving in the l = 1 and l = 2 angular momentum modes, combined to have the total angular momentum zero. Thus their finite volume mass correction should have an extra orbital term (see Eq. (3.10)) As is clear from figure 16, E i − E 0 − 2M decrease way too slowly with R to be described in the asymptotic region by just the sum of the curvature and orbital corrections. It must be that the difference is due to the scattering correction, although it looks hard to make a quantitative conclusion using the data in the R < 3 region.
1.00 0. 50 2.00 0. 30 3.00 1.50 0.70 Figure 16: Blue curves: E i − E 0 − 2M for the lowest two orbital excitations in the twoparticle sector of the λ = 0.3 spectrum corresponding (log-log scale). We excised by hand a state which asymptotes to 4M and thus looks like a four-particle state at rest. Dashed black lines: 2(∆M curv + ∆M l ) for l = 1, 2.
Excitations for λ = 0.55 A very different behavior presents itself in the spectrum dependence on R for λ = 0.55, which is close to the critical coupling. Instead of energies tending to finite limits, we see them all gradually decrease with R.
As already mentioned in section 7.2.2, at λ = λ c we expect the excitation energies to scale at large R as ∆ IR i /R, where ∆ IR i are the IR CFT operator dimensions. To test this expectation, we plot in figure 17 the excitation energies times R. We vary λ in the range 0.53 . . . 0.56, roughly the range determined in section 7.2.2 to contain the critical coupling. Apart from the state in the Z 2 -even sector, of dimension ≈ 1.175 (see Eq. (7.6)), we expect to see a Z 2 -odd operator σ of dimension ≈ 0.305 (as extracted again from [24,7]). We also expect to see the states corresponding to operators ∂ 2 and ∂ 2 σ, of dimension two units higher. Finally, we may hope to see the next primary Z 2 -even operator , whose dimension for d = 2.5 is not precisely known but may be expected to lie between 3.5 and 4. 38 Interestingly, for R ≈ 3 we can observe all of the above mentioned states in the spectrum in figure 17, at the dimensions where they are supposed to be and with the right Z 2 quantum number. The agreement of theory and our numerical results remains imperfect in that the curves 38 It's 4 in d = 2 and in d = 4 − , and ≈ 3.83 in d = 3 [26][27][28][29]. . The curves must tend to the indicated finite limits. They do reach these limits for R ≈ 3, but would overshoot them (except for σ) for larger R.
don't really approach finite limits very well. Perhaps one could claim that for the lowest two states σ and , whose variation with R is not huge. However, the higher states definitely exhibit growth with R and would overshoot the theoretical prediction for their dimension, were we to extend this plot to higher values of R. We hope that this issue will get resolved in the future by improving the accuracy of the method (see section 7.3).
Excitations for λ = 0.9 Finally, we consider the spectrum for λ = 0.9. 39 The first eye-catching feature of the spectrum is the approximate degeneracy of Z 2 -even and odd states in the region of large R. This degeneracy is clearly visible for the first excitation, which becomes degenerate with the vacuum, and for two more pairs of states. The interpretation of this phenomenon was already discussed in section 7.2.1-it means that the Z 2 -symmetry is spontaneously broken. In the infinite R limit we expect a pair of degenerate vacua |0 L,R , each with its own tower of excitations |i L,R . The spectrum is thus exactly doubly degenerate. For a finite large R, the height of the potential separating the two vacua is no longer infinite-it is expected to scale as the volume of space R d−1 . This allows tunneling and the eigenstates become the Z 2 -even and odd linear combinations: of energies E i ± δ i where δ i is the tunneling matrix element. As usual in quantum mechanics, we expect that the Z 2 -even combination will have a smaller energy than the Z 2 -odd one. This is confirmed by the spectrum in figure 13-in each of the three approximately degenerate pairs, it's the Z 2 -even state which is the lower one.
The above tunneling argument seems to predict that the even/odd state pairs should be split roughly symmetrically with respect to the infinite volume limiting value. In fact, since the excitation energies are defined as E i − E 0 and E 0 belongs to the Z 2 -even sector, we expect the Z 2 -even excitations to shift down by (δ i − δ 0 ), while the Z 2 -odd ones to move up by (δ i + δ 0 ). Since the tunneling probability strongly depends on the energy, we expect δ 0 δ i , and the shifts should be roughly symmetric. However, that's not what we see in the λ = 0.9 plots in figure 13-it rather looks that the negative shift of the Z 2 -even excitations is much larger than the positive shift of the Z 2 -odd ones. For the first pair of massive excitations, it even looks like both the Z 2 -odd and the Z 2 -even state have a negative shift.
The most natural explanation of this phenomenon is that we are forgetting the modification of the mass spectrum via coupling to curvature, see Eq. (7.9). This effect goes as 1/R 2 and for low-lying states should be larger than the splitting, which is exponentially small in the volume of sphere. The fact that both states in the first Z 2 -odd/even pair have a negative shift can then be explained by taking A negative in Eq. (7.9).
In fact, it is not totally unexpected that A should be negative for the lowest massive excitation at λ = 0.9. The same occurs for the Landau-Ginzburg flow in the weakly coupled part of the Z 2 -broken phase, i.e. for negative m 2 and a small quartic coupling. We did not study this part of the phase diagram numerically, but it's easy to understand what happens analytically. The full mass parameter of the UV theory, including the coupling to curvature, is m 2 + ξ Weyl Ric. Since m 2 < 0, we have to reexpand the Lagrangian around the true vacuum, and when we do this, the mass parameter picks up the usual −2 factor. We thus conclude that ξ = −2ξ Weyl , giving A = −2 at weak coupling in the Z 2 broken phase.
This finishes the discussion of splittings at finite R. Next, we would like to say a few words about the overall structure of the massive spectrum at large R. We identify the two near-degenerate pairs of even/odd states with two massive excitations, of mass A very interesting feature of this spectrum is that M 2 < 2M 1 . This is unlike the spectrum at λ = 0.3, which was neatly quantized in the units of the lightest excitation. In the situation at hand, the state of mass M 2 should probably be interpreted as a bound state of two M 1 particles.
The appearance of such bound states was found long ago, and their masses measured, in the lattice simulations of the broken phase of the Ising model and of the φ 4 theory in d = 3 dimensions [30]. In the weakly coupled regime, the existence of these states follows from the fact that the λφ 4 interaction becomes attractive in the broken phase, the cubic exchange diagrams overwhelming the repulsive effect of the contact term interaction [31]. Their binding energy, exponentially small at weak coupling, is known in the leading and first subleading exponential approximation [31,32]. Apparently, here we are observing the same effect in d = 2.5 dimensions and at strong coupling.

Renormalization Details
The renormalization in the Landau-Ginzburg flow is determined through the leading OPEs of the deforming operators in (7.1), given by where In the RHS of (7.14) we omitted the operators whose associated function B(h) vanishes because of the remark in footnote 29. All the retained operators have nonzero B(h) and will be relevant for the renormalization. The effect of φ 2 and φ 4 in the RHS will be to make the couplings λ and m 2 non-trivial functions of the cutoff. We are in a position to use the RG-improved formalism of section 5.3.
A crucial ingredient in the renormalization is the relation between the above OPEs and the asymptotics of the matrix (M n ) i j in (5.12). For the φ 2 flow we were able to determine the asymptotics of this matrix exactly including the discrete structure, see equation (6.11). In the future, such exact asymptotics may be also worked out with the φ 4 coupling switched on, see the end of Appendix D for a discussion, although the task looks more challenging. In this work, we will use the continuum approximation (5.18). We checked the accuracy of this approximation for many choices of external states i, j against the exact expression for M n within the range ∆ n ∆ max where we know V and can compute M n numerically. These checks convinced us that the approximation is adequate.
One such check is shown in figure 18, where we plot both the exact and the approximate behavior of (M n ) i j for V = : φ 4 : and for O i = O j = : φ 2 :. The blue dots represent the individual values (notice that M n is nonzero only for half-integer ∆ n ). The red dashed line shows moving average of these values within the interval [∆ − 1, ∆ + 1]. The solid black line is our approximation as given in (5.18), including the contributions of all three leading operators in the φ 4 × φ 4 OPE shown in (7.14). We see that the agreement between the moving average and the approximation becomes very good at ∆ ∼ 17, which is also the cutoff we used in this study. The agreement for other choices of O i and O j is similarly good.
From the OPEs (7.14) we can directly generate the RG equations discussed in section 5.3 for the couplings of the local operators. They are given by: where we denoted by g 0 , g 2 and g 4 the coupling associated to 1, : φ 2 :, and : φ 4 :. We also introduced with the OPE coefficients f abc given in (7.14) and B(h) was defined in (5.17). The renormalized couplings are then found by integrating these equations numerically from Λ = ∞ to the desired value of the cutoff Λ = Λ UV = ∆ max /R. We impose boundary conditions at infinity such that g 0 = 0 and g 2 and g 4 are given by their bare UV values: As we explained in section 5.3, the above RG equations depend on a reference energy E r . In our study we found it convenient to choose E r to be around the energy of the first excited state in the Z 2 even sector. An estimate for this energy was obtained by extrapolating the earlier obtained results for lower values of the radius or the coupling, or by performing a quick computation with a smaller ∆ max .
We also discussed in section 5.3 the subleading dependence on (∆ i +∆ j )/R and onĒ −E r . This dependence is taken into account by adding to the correction Hamiltonian the additional non-local terms given in (5.22). Their coefficients, which we denote as g The equations for other i are determined by modifying the corresponding equation in (7.15) in a similar manner. We integrate these equations with boundary conditions zero at infinity. Notice that we ignore the backreaction of these terms in the sense that they do not appear on the righthand sides of the flow equations. 40 In figure 19 we present an example of the flow of several couplings. We see that g 2 receives substantial corrections, demonstrating the need for the RG improvement. On the other hand, we see that the relative change in g 4 is small, and that g (H CFT ) 0 remains small compared to 1 (coefficient of H CFT in the bare Hamiltonian) throughout the flow. In this example we set R = 3, and λ = 0.7 and m 2 = 1 in the UV. We usedĒ r = −6.
With the correction terms described up to now, the results we obtained already looked reasonable. We were however able to take into account one further correction, which turned out to give a noticeable improvement mainly for small values of λ. Namely, we constructed a new nonlocal counterterm which completely takes into account the dependence on (∆ i + ∆ j )/R and (Ē − E r ), beyond expanding to first order as above. The coefficient of this correction is computed by running the RG flow again (once more ignoring the backreaction of this nonlocal term) but separately for each value of κ = 1 2 (∆ i + ∆ j )/R andĒ. In equations, this means that we update this extra nonlocal correction in each RG step as follows, .
Inside curly brackets, we subtract the zeroth-and first-order terms in κ andĒ − E r , since these terms were already taken into account above.
Let us recap. We integrate all the above RG equations from Λ = ∞ to Λ = Λ UV and obtain the correction terms. We divide them into four groups: 40 For g (HCFT) 0 and g (H) 0 , which multiply local terms in the Hamiltonian, it would be possible to incorporate such a backreaction rather easily. At every step of RG one should factor out the modified coefficient of H CFT , which leads to an overall rescaling of the remaining couplings. The product of all rescalings Z should be stored separately to undo the rescaling at the end of the computation. This procedure resembles wavefunction renormalization in perturbative RG. We implemented it, but found the numerical effect of this improvement to be very small.
• ∆H loc which reflects the change in all local couplings; • ∆H 1 and ∆H 2 which include the nonlocal terms proportional to It is now time for numerical diagonalization. We add ∆H loc and ∆H 1 directly to the bare TCSA Hamiltonian, and diagonalize. Let's call the resulting eigenvalues and eigenvectors E n and c n . In principle, we would also have liked to add ∆H 2 before the diagonalization. Unfortunately, we found that doing this destabilizes the numerics. This instability must have its origin in the fact that the factor (H − E r ) is not small for states of high energy, and even for states of low energy it's not manifestly small, being a difference of two separately large quantities. We therefore chose to add the effect of ∆H 2 only after the numerical diagonalization. We found it necessary, and sufficient, to do this to the second order in ∆H 2 . The correction is computed by the usual Hamiltonian perturbation formula: The sum over m in the second-order term is rapidly convergent, and it's enough to sum over the first few eigenstates. Notice that one has to appropriately insert the right and left eigenvectors. This is not reflected in the notation but explained in detail in section 6.4.
Finally, we compute one last correction due to ∆H nl , which turns out to be very small, so doing it to first order is sufficient: When evaluating this correction, we are supposed to setĒ in the definition of ∆H nl to the energy of the state we are correcting. Also recall that ∆H nl depends on κ = 1 2 (∆ i + ∆ j )/R, and this dependence comes into play when evaluating the scalar product.
This completes the description of the renormalization procedure used to produce the plots in section 7.2. As the above discussion shows, an efficient implementation of the leading-order truncation effects given in (5.21) is subject to various subtleties, mostly due to the non-negligible presence of (∆ i + ∆ j )/R andĒ in the correction terms. In this exploratory paper we have not aimed to present a complete analysis of these effects. Instead, we discussed various recipes for dealing with them at a practical level. The details we provided should be sufficient to reproduce our results. In the future it would certainly be interesting to perform a more systematic study of all the subtleties.

Non-Unitarity and Complex Energy Levels
As we observed in section 3.3, the free massless scalar theory in fractional d is not unitary-its Hilbert space contains negative norm states. In d = 2.5 the lowest negative-norm state occurs for ∆ = 9. In figure 20 we show the total number of scalar states and the number of negative norm states as a function of ∆. What are the consequences of having these negative norm states? One expected consequence is that once we perturb the theory, we will get complex eigenvalues. The purely massive perturbation 1 2 m 2 φ 2 is an exception, since in this case we expect that the spectrum agrees with the canonical quantization spectrum from section 3.5, and thus is real. 41 What if we turn on λφ 4 ? As we saw in the previous sections, numerics indicate that the low-energy spectrum is still real. This may not be so surprising, since the negative norm states are secluded at high energies. So to see complex eigenvalues, we may expect to have to go to high energies. We will now present several computations which show that complex eigenvalues do occur.
Let us first of all examine the case of very small R. In this limit we can treat m and λ as perturbations, with dimensional couplings m 2 R 2 and λR d−4ν . The second coupling decreases less slowly as R → 0, and will dominate at very small R. The effects of the perturbation is to split the degenerate energy levels of the CFT Hamiltonian. The splittings are proportional to the eigenvalues of the perturbation diagonalized within each degenerate subspace. In high energy subspaces, which contain negative norm states, some of the eigenvalues may and do turn out to be complex. We find that this happens for the first time at ∆ = 11.5, which is an 88 dimensional subspace with 7 negative norm states. We find that the matrix of the φ 4 perturbation within this subspace has one pair of complex conjugate eigenvalues ≈ 1.85 ± 0.04i. This implies that for very small R the energy levels will be complex.
As a side remark, we note that it would be interesting to carry out a similar computation in 4 − dimensions. Free massless scalar theory in 4 − dimensions perturbed by λφ 4 flows to a weakly coupled Wilson-Fisher conformal fixed point. This is a short flow, and the energy levels at the IR fixed point, which are in one-to-one correspondence with the IR operator dimensions, are computable in perturbation theory. Once again, we expect that some of the IR operator dimensions will be complex. Likely this will happen already to first order in , and it would be interesting to identify the first operator dimension for which this happens. The first null state in d = 4 has dimension 15. To get a complex anomalous dimension to first order in we have to go to even higher ∆. This computation is a bit more difficult to perform than the above computation in d = 2.5, because in d = 4 − there are O( ) splittings already in the unperturbed spectrum (like between (∂φ) 2 and φ 4 ). So one should use the near-degenerate rather than degenerate perturbation theory. We will come back to this question in future work [33].
Going back to d = 2.5, the above argument is confirmed numerically in figure 21, where we show the spectrum around ∆ = 11.5 for m 2 = 1, λ = 0.55, and 0 < R < 0.15. We see precisely one pair of complex conjugate eigenvalues emerging out of the ∆ = 11.5 group for small R. For larger R, the spectrum shows intricate structure. We see many beautifully resolved level crossing avoidances in the real part of the spectrum. We also see a second pair of complex conjugate eigenvalues appearing at R ≈ 0.04 and then disappearing at R ≈ 0.07. Zooming in on this line of complex eigenvalues, one notices that it joins collision points for pairs of real eigenvalues. This last observation may seem to create a minor paradox. Didn't we say that the Landau-Ginzburg flow is not integrable, and that in non-integrable flows energy levels don't cross? The resolution is that this last statement requires a qualification in presence of negative norm states. If two energy levels which head for a collision are both positive-norm (or both negative-norm), they will generically repel. However, in a subspace with non-sign-definite Gram matrix, no-levelcrossing rule does not apply. To see this, consider a toy-model 2 × 2 symmetric generalized eigenvalue problem where σ = ±1 depending on whether we are dealing with a subspace of positive or non-sign-definite norm. We are assuming that the Hamiltonian matrix is symmetric and real. The distance between the two eigenvalues is controlled by the discriminant: For σ = 1 the discriminant is a sum of two squares, and level crossing cannot generically happen. On the other hand, for σ = −1 the discriminant is not positive definite, and can readily change sign if the off-diagonal matrix element increases beyond a critical value. When this happens, we go from having two real eigenvalues to a complex conjugate pair.
As another illustration, in figure 22 we show the spectrum with R E i ∈ [8,11] for the same couplings as above but in a wider range 0 < R < 2. We clearly see several points where real eigenvalues collide and form a complex conjugate pair, which sometimes reemerges as a pair of real eigenvalues for a slightly larger value of R. The most prominent such collision happens at R ≈ 0.7. Figure 22: The spectrum at R E i ∈ [8,11] for m 2 = 1, λ = 0.55, and 0 < R < 2 with a step of 10 −3 . We are plotting energy levels multiplied by R. Light blue: real eigenvalues. Black: real part of eigenvalues with nonzero imaginary part. These are raw TCSA data with ∆ max = 11. The jittering spread noticeable in some of the eigenvalue curves at R 1.5 is due to numerical instabilities in the Mathematica diagonalization routine.
To resolve the multitude of eigenvalue curves in figures 21, 22, we had to compute the spectrum with a very small R step. For reasons of speed and numerical stability, we have performed these bulky computations with a relatively small ∆ max . Since the complex eigenvalues observed in these plots lie relatively close to the cutoff, their energies are likely to shift considerably when the cutoff is increased. However, we don't expect the complex eigenvalues to disappear. In fact, we performed checks for a few selected values of R, computing the spectrum with a higher cutoff and, for higher numerical stability, with a higher number of digits rather than at machine precision. The complex eigenvalues were always present.
Notice that the eigenstates corresponding to the complex eigenvalues necessarily have zero norm (computed with respect to the Gram matrix). In a unitary theory a state of zero norm has zero overlap with any other state. Such a state is unphysical; it can be kept in the Hilbert state or thrown out without physical consequence. This was the situation with the scalar theory in d = 3 in section 6, whose extended Hilbert state included some null states, but only as a matter of convenience. In a non-unitary theory, as the one we are discussing now, states of zero norm do not in general have zero overlap with other states. They cannot be removed from the theory without modifying it.
To summarize, the Landau-Ginzburg theory in d = 2.5 dimensions is a non-unitary interacting quantum field theory. Its spectrum on a sphere of finite radius contains negative-norm states with real energies, as well as zero-norm (but physical) states with complex energies. The negative-norm and zero-norm states belong to the high-energy part of the spectrum, and so their effect on the low-energy physics may not be huge, but the mere presence of these states is a proof that the theory is not unitary. We expect complex eigenvalues to be present also in the limit R → ∞. In particular, the critical point of the theory should have operators with complex scaling dimensions. The same should be true for theories in any fractional d. As mentioned above, we expect that these facts are not difficult to check in 4 − dimensions.

Discussion
We have presented an initial study of the feasibility of the TCSA for theories in more than two spacetime dimensions. Both for the massive flow and for the Landau-Ginzburg flow we have obtained promising results. It appears that the TCSA does work, and that its region applicability extends far beyond the realm of perturbation theory.
The most important challenge to the TCSA is the exponential growth in the number of states, see figures 2 and 20 for examples. This exponential growth cannot realistically be overcome by simply employing more computational resources. Consequently, it appears difficult to raise the cutoff much beyond the values used in this paper. Improving the accuracy of the TCSA method will therefore hinge on our ability to invent techniques that circumvent this problem. This was realized already in investigations of the TCSA in d = 2, see appendix A for a summary. In this paper we have approached the issue by following standard renormalization principles, and found that adding the right, analytically computed, counterterms can drastically improve the numerical results. Our methods worked well in the given examples, but a more systematic investigation is certainly called for.
The number of possible problems which can be attacked with the TCSA is huge. Where to start? In this paper we considered the Landau-Ginzburg theory in a sufficiently low number of dimensions (d = 2.5) where it is UV finite. This was reflected in the fact that the counterterms that we had to add were proportional to negative powers of the UV cutoff. The next natural step would be to carry out the same study in d = 3 dimensions, where the perturbing operator φ 4 will have dimension 2. Since this is above d/2, the ground state energy will require infinite renormalization (it will be linearly divergent). In addition, the mass term will be logarithmically divergent. The presence of these divergent leading counterterms implies that also the first subleading counterterms, suppressed with respect to the leading ones by only one power of Λ UV , will be numerically more important than in d = 2.5. It would be interesting to see if good numerical accuracy can be achieved in spite of these complications. This will likely require further improvement of the renormalization procedure.
After φ 4 in d = 3, the next step might be to consider Yukawa interactions. In principle, dealing with fermions (including chiral fermions) in TCSA is completely straightforward. Because in d < 4 the Yukawa interaction is less relevant φ 4 (in d = 3 it has dimension 2.5), the cutoff dependence for Yukawa theories is expected to be more severe than for the Landau-Ginzburg flow, and achieving good accuracy will probably be more challenging.
After the Yukawa theory, one could try to use TCSA to study the 3d gauge theories. The easiest problem in this class might be three-dimensional QED-the U (1) gauge theory with fermionic matter. The interaction term in the Lagrangian, A µψ ψ, has dimension 2.5, and so the cutoff dependence will be as severe as for the Yukawa interactions. In addition, it may not be entirely straightforward to treat the gauge interaction term in TCSA. Such issues obviously need to be resolved before one can start attacking nonabelian gauge theories in d = 4.
Another problem for the future is as follows. In this paper, we were focusing on extracting the mass spectrum of strongly interacting QFTs. Once the accuracy of the method is improved, and the mass spectrum is under total control, it will make sense to turn to the problem of recovering the scattering matrix, from Lüscher-type corrections to the masses on a sphere of finite radius. As we mentioned in section 7.2.3, the theory of these corrections is not yet available, and it would be interesting to work it out.
Let us conclude on a philosophical note. Few methods are available to study non-supersymmetric strongly coupled quantum field theories. Some of these, like the gap equation, are analytical but merely qualitative, and hardly more reliable than dimensional analysis. On the other hand there are lattice measurements, which are performed from first principles but require vast computer resources. In this paper we pointed out that there exist alternative algorithms, which are computationally much cheaper but nevertheless defined with a level of mathematical rigor equal to that of the lattice. Here we highlighted the promise of the TCSA. However, the TCSA is just one representative of a family of Hamiltonian truncation methods in quantum field theory. In Appendix B we review the existing work on related methods, some of which also look promising and are currently under active development. It seems worthwhile to explore these methods, perhaps not necessarily in order to compete with the lattice but rather with the hope that they can significantly reduce the time and resources required to answer interesting nonperturbative QFT questions. We must break free from the view that some questions can only be answered by the lattice, roll up our sleeves and start constructing alternative tools.
work. We will also point out the differences between the way the computations are done in d = 2 and in d > 2.
TCSA was introduced by Yurov and Al. Zamolodchikov in [1], where they used the method to study the Lee-Yang CFT M 2,5 perturbed by the only relevant scalar of the theory, of scaling dimension 42 ∆ = −2/5. Even though this theory is non-unitary, it was a perfect example on which to demonstrate the method, because: 1) the UV CFT spectrum is very sparse; 2) the perturbing operator is very relevant, so that excellent convergence is achieved for very modest cutoffs; 3) the flow is integrable, which allowed comparing the predictions to an exact solution in the IR.
The next important paper was [22] which studied flows originating at the tricritical Ising model M 4,5 . This model has several interesting relevant perturbations. Using TCSA, they managed to chart out the intricate phase structure in the IR largely in agreement with expectations. They were the first to observe level repulsion for non-integrable flows, and spontaneous global symmetry breaking (via exponential degeneracy of ground states in finite volume). They also pointed out that for a given UV cutoff the method fails above certain R and gave a practical criterion to determine this R, by looking for a change of exponent in the eigenvalue dependence on R.
One case where Ref. [22] encountered a difficulty was the perturbation by the subleading energy operator of dimension ∆ = 6/5, which is known to flow in the IR to the Ising model CFT M 3,4 . They observed a large UV cutoff dependence of the results which was obscuring the IR behavior. The basic reason for this was soon explained in Ref. [34]: only for ∆ < d/2 = 1 is the perturbation UV finite, and we can expect naive TCSA to converge as the cutoff is taken to infinity. For ∆ d/2 UV divergences appear, and the TCSA Hamiltonian needs to be renormalized. For ∆ just above d/2, the only UV divergence is in the vacuum energy density. Since this affects every energy level in the same way, a quick fix is to consider energy level differences. As ∆ is raised further, more nontrivial divergences are expected to appear; it was not discussed at the time how to deal with them.
After a several year's pause, a series of interesting papers appeared where TCSA was applied to perturbations of the free scalar boson theory (up to that moment only minimal model perturbations had been studied). In [35,36] the sine-Gordon perturbation was studied, and an agreement with the exact spectrum known from integrability was observed. In [37], the two-frequency sine-Gordon model was studied. This model is non-integrable and has a phase transition in the 2d Ising universality class, which the authors were able to locate and study using TCSA. An integrable SUSY sine-Gordon model was studied using TCSA methods in [38].
We would also like to mention a nice paper [39] (see also [40]), which first discussed how one can use Lüscher corrections to extract scattering phases from TCSA data.
How does one do TCSA computations in d = 2? The space of local scalar operators of the UV CFT is obtained by acting with raising operators on the Virasoro primaries: In minimal models, some of the states constructed in this way will be null; they are usually separated away. In early works [1,22], one also separated the states into quasi-primaries and higher Virasoro descendants. This is not strictly necessary, and the required computation effort may easily outweigh the subsequent speedup in the evaluation of the Gram and the Hamiltonian matrices. When we built a code to check some of the d = 2 results, we avoided doing this step.
In the free scalar theory, one can construct states using the U (1) Kac-Moody algebra rather than the Virasoro algebra; this is simpler since the U (1) Kac-Moody primaries are just the exponentials exp(iαφ) while there are many more Virasoro primaries.
One computes the Gram matrix for these states G ij = i|j using the Virasoro (or Kac-Moody) algebra to commute all lowering operators to the right until they hit the primary. Finally, one computes the matrix elements of the perturbing operator, i|φ(1)|j . Here one uses the commutation relations with a Virasoro primary [L n , φ(z,z)] = z n [(n + 1)h + z∂ z ]φ(z,z) , (A. 2) or an analogous Kac-Moody relation. Clearly the conformal algebra plays a big role in the d = 2 TCSA computations, while it did not feature prominently in our discussion of TCSA in d > 2.
The earliest TCSA computations were done in Hilbert spaces of a very modest size: 17 in [1], about 200 in [22]. The sine-Gordon papers cited above used much larger Hilbert spaces of a few thousand states. Trying to increase the number of states further by brute force is a game of rapidly diminishing returns, since the number of states grows exponentially with the cutoff, while the error goes down only as a power law.
How can one tame the growth of the number of states? It is natural to employ the ideas of renormalization in this context, and the current literature already contains several proposals. One approach [41] (described in more detail in [42]) is Numerical RG, inspired by the namesake method used in Wilson's famous solution of the Kondo problem. The idea is to add states to the Hilbert space in manageable batches, and after each addition rediagonalize the Hamiltonian and throw out the least important states so that the total number of retained states never grows more than a few thousand, while the total number of explored states may be several orders of magnitude larger. This procedure is not exact: some systematic error is accumulated because the subsequent batches of states are not allowed to talk to each other directly, but only through the retained states. Empirically, this error seems to be small, since numerical RG procedure is in good agreement with the exact results when available. It would be nice to better understand theoretically why this happens, and to get a quantitative estimate on the error.
An alternative opinion, which is also ours, is that one should invest more effort into purely analytic approaches to renormalization. The point is that, if the theory is weakly coupled at the UV cutoff, one should be able to integrate the energy states above the cutoff analytically, and explicitly construct the correction terms needed to improve the accuracy of TCSA. This idea was first studied in [14,15] (see also [43]). This was done in the context of boundary RG flows, which can be studied via a variant of TCSA. For bulk flows, analytic renormalization was discussed in [13] and more recently in [16] (see also [44], which in particular includes a study of marginally relevant perturbations). Our discussion of renormalization in section 5 was inspired by [13], although, as we discussed in section 5.4, our method is different in several aspects from [13,16].

B Other Hamiltonian Truncation Techniques
We would like to list here works which applied other Hamiltonian truncation methods in quantum field theory. These methods are conceptually close to TCSA even though their implementations may be quite different, for example because the UV theory is massive instead of conformal. The earliest such work known to us is [45], which studied a 2d Yukawa model-a massive scalar and a massive fermion with a Yukawa interaction yφψψ-in the Hilbert space of the free massive theory on a circle of length L.
Then, in a series of interesting papers [46,47] (see also [48]) this approach was applied to the massive φ 4 theory in d = 2 and d = 3 dimensions. The reader should compare the spectra in figure  8 of [46] (d = 2) and figure 1 of [47] (d = 3) with our figure 10 in d = 2.5. We believe that it is worth revisiting the approach of [46,47] and possibly improve on some of the implementation details. 43 As we mentioned in section 3.5 we are planning to do this at least in d = 2 where TCSA is not directly applicable [? ].
Then there is a series of papers by Fonseca and A. Zamolodchikov [49][50][51] on the "Ising field theory", which is the 2d Ising model perturbed by both relevant operators m + hσ. In principle, this model can be studied using TCSA. However, these papers find it more efficient to take the perturbation into account exactly from the start, using the fact that it corresponds to turning on mass in the fermionic description of the 2d Ising model.
Finally, we would like to mention a large body of work which apply Hamiltonian truncation methods to field theories quantized on the light front. A time-honored technique is DLCQ, which compactifies x − to get a discrete spectrum for the Hamiltonian evolution in the x + direction [52,53]. While this idea has been around for a while (see [54] for a review), it has not yet become a viable alternative to the lattice above d = 2 dimensions. In 2d it did lead to several interesting analyses, for example of the φ 4 theory [55], or of the large N QCD with matter in the adjoint [56].
Recently, another approach to Hamiltonian truncation on the light front was proposed in [57,58] in the context of 2d QCD (large or finite N ) coupled to massless matter. Instead of compactifying a light-cone direction, they use a discrete basis of multiparton wavefunctions, which are in one-to-one correspondence with the left-moving quasiprimary operators of the matter CFT. This appearance of conformal operators creates at least a superficial similarity with TCSA, although it seems that the light-front physics is rather different from that in a finite volume. At any rate, it turns out that the conformal basis of [57,58] leads to a dramatic improvement in the convergence rate of the method. Empirically, convergence seems to be exponentially fast. In comparison, DLCQ is always affected by 1/(lightcone volume) corrections. As we saw in this work convergence rate of TCSA is also power-like in the cutoff. It would be extremely interesting to understand if the successes of the light-front conformal basis can be extended to 2d theories where the right-and left-moving sectors talk to each other, and to theories in d > 2.

C Asymptotics of C(τ )
In this appendix we will give some details concerning the derivation of Eq. (5.17). We start from the definition of C(τ ) in Eq. (5.13). We do the Weyl transformation which maps the correlation function on the cylinder into one in flat space. In flat space the operation insertions look like: is the product of factors picked up by the operators under the Weyl transformation. To the shown order in τ , which is the one we need, the effect of w will average to zero when summed over a ↔ b. We next use the flat space OPE: Instead of inserting the RHS operator at the middle-point as in (5.15), we put it at one of the endpoints, since this facilitates the subsequent integration. However, the needed accuracy then requires the inclusion of the shown first subleading term.
We now have to do the integral over y running over the sphere of radius 1/r = e −τ /2 . By rotation invariance we can take x = ( 0, r). The nonanalytic behavior of the integral as τ → 0 will come from the region of the y sphere closest to x, i.e. its northern cap, which we parameterize as y = ( ρ, e −τ − ρ 2 ) , |x − y| 2 ≈ ρ 2 (1 + τ ) + τ 2 , (C.4) where we kept the approximation needed to get C(τ ) to O(τ 2 ). We are led to evaluate the integral D M n sequence for φ 2 × φ 2 The purpose of this appendix is to derive Eq. (6.11), which gives the exact asymptotics for the M n sequence in the case of the φ 2 perturbation. We thus consider the matrix V i j defined as We will describe two methods to get the answer. The first one is direct: we will study the matrix elements V k j and V i k and identify the states of energy ∆ k ∆ i , ∆ j which contribute to the sum defining M n . To compute V k j , we consider the OPE φ 2 (x) × O j (0). By Wick's theorem, we can write: : φ 2 (x) : O j (0) = : φ 2 (x)O j (0) : Here in the second and third line we put terms where one or two φ's out of φ 2 (x) are contracted with the φ's making up O j , which can possibly carry several derivatives denoted collectively as (α), (β). The operatorsÔ j andÔ j are O j minus the contracted parts.
The operators in the third line have all dimension < ∆ j , so they are not relevant for the asymptotics. The operators coming from the second line will appear by expanding φ(x) under the normal ordering sign into the Taylor series and picking up terms which will not vanish upon integration over the unit sphere, with φ(x)∂ (α) φ(0) as a weight. A moment's thought shows that the only surviving operators will be : ∂ (α) φ(x)Ô j (0) :, i.e. where φ(x) is expanded at order α. These operators have dimension ∆ j and are also irrelevant for the asymptotics.
We conclude that the large dimension states appearing in the OPE are the states | p (φ 2 )O j , whose dimension is ∆ j + ∆(φ 2 ) + 2p. To complete the calculation, we need to compute V i k for k being one of these states and O i an operator of low scaling dimension. We have The second method of deriving Eq. (6.11) proceeds by analyzing the correlation function C(τ ), defined in Eq. (5.13). Let us focus on one particular contribution to this correlation function, when we contract the φ's in the φ 2 insertions among themselves, and the φ's in O i and O j among each other. Going to flat space, the integral over n and n in (5.13) can be performed exactly, giving a hypergeometric function 2 F 1 . It turns out that expanding this 2 F 1 in a series in r and reading off the coefficients gives exactly the sequence (6.11). This means that all other ways of contracting φ's give contributions to the correlation function whose expansion in r does not give rise to terms with arbitrarily high powers. These other options would involve only one or zero contractions between the two φ 2 operators. One can show that for these terms the 2 F 1 truncates to a polynomial, and therefore they do not correspond to an infinite sequence of scalar operators.
The second way of computing the asymptotic behavior of M n looks much more economical than the direct method given above. It must also be easier to generalize to other perturbing operators. This would allow one to construct a renormalization procedure for the Landau-Ginzburg flow which takes into account the discreteness of the sequence M n , and is thus more accurate than the one described in section 7.3.