Truncated conformal space approach in d dimensions: A cheap alternative to lattice field theory?

We show how to perform accurate, nonperturbative and controlled calculations in quantum field theory in d dimensions. We use the truncated conformal space approach, 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 quantum field theory Hamiltonian is expressed as a matrix in the Hilbert space of conformal field theory 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 ϕ 4 theory in d dimensions with d being 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 nonintegral d are not unitary, which however does not seem to cause much effect at low energies.


I. 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 are 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 these 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 highenergy 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 Sec. II with a general discussion of TCSA and of its connection to the Rayleigh-Ritz method in quantum mechanics. Section III 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 nonunitary theory-its Hilbert space contains negative norm states. This curious fact does not seem to have been noticed before.
We continue in Sec. IV with more details about TCSA for perturbations of the free massless scalar. In particular, we describe an operator product expansion (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 V 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 Secs. VI and VII we apply TCSA to study the Landau-Ginzburg flow. We first discuss in Sec. VI 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 Sec. VII 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 nonunitarity of the theory at fractional d manifests itself via some high-energy eigenvalues acquiring imaginary parts.
We conclude in Sec. VIII with a discussion of future research directions. Appendixes 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 appendixes contain derivations of auxiliary technical results.

II. TRUNCATED CONFORMAL SPACE APPROACH
We study a 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 jii 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

VðxÞ: ð2:2Þ
The key idea is to think about this Hamiltonian as an infinite matrix in the Hilbert space of unperturbed CFT states jii. The CFT Hamiltonian in this basis is diagonal and simply related to the CFT operator dimensions 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 where f O † i VO j 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 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, (i) the ground state dependence on R gives the vacuum energy density; (ii) by looking at the number of exponentially degenerate ground states we can infer the symmetrybreaking pattern; (iii) the excited states give the massive spectrum of the theory, including one-particle, many-particle, and bound states; (iv) studying the dependence of two-particle state energies on R we can extract the S-matrix; (v) 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 1989 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. 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 jni of the harmonic part H 0 . Then truncate the basis by keeping only the states jni 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: 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 4 This is not a generic feature of the method. It is related to the fact that q 4 raises/lowers energy by at most a finite amount (four units). In more complicated quantum mechanical examples, and in quantum field theory, perturbation matrix elements will have powerlike tails at high energy. In this case the convergence rate will be powerlike in the cutoff, as we discuss below. 5 For large λ, exponentially fast convergence sets in once we include all harmonic energy levels below the anharmonic one we are trying to reproduce. 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 wave function 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], Sec. 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 wave function 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.

B. 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 the perturbation problem is in the UV. The best situation is realized when which for the perturbations considered here means 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 way-by 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.

III. 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.

A. 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.
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 6 One way to work around this problem is to compactify the scalar on a circle of large radius, considering a periodic approximation to the Landau-Ginzburg potential [4]. Note added. The paper [5] discussing such an approach was submitted to the arXiv the same day as our work. 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 Fig. 1. It is also clear that isomorphic graphs correspond to identical operators, and should not be counted separately.

B. 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 will not 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 d ¼ 3 is ϵ 123 δ 45 δ 67 ϕ ;1 ϕ ;24 ϕ ;356 ϕ ;7 : ð3:2Þ Consequently, the P ¼ −1 states in the IR are also likely to be heavier than for P ¼ þ1.
Second, it is 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 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 trM n , n ≥ 4, are expressible via linear combinations of the products of trM 2 and trM 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 nonisomorphic 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.
C. Nonunitarity 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 FIG. 1 (color online). The graph corresponding to the only scalar operator which can be obtained by contracting indices in (3.1). 7 And it may be interesting to check that their energies relate to the scalar energies in agreement with the d-dimensional Lorentz invariance which should emerge in the large R limit. 8 We can also replace n parallel edges by a single edge with n as a label. Then our graphs become simple edge-colored graphs. 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 ð3:5Þ 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 Sec. VII D, 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 Sec. III E.
It has to be said that the first negative-norm state has a pretty high dimension: 10 þ 5Δ ϕ ð3 < d < 4Þ: ð3:6Þ 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

D. 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 nonperturbative conformal data, while matrix elements involving descendant states can be computed using the conformal algebra. This situation often occurs in twodimensional 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.

E. 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 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 Sec. VI 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 do not 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= ffiffiffiffiffi ω l p 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 an upcoming paper [9], we will use canonical quantization to study the Landau-Ginzburg flows in d ¼ 2. As explained in Sec. II B, TCSA is not directly applicable to these flows (see however footnote 6 ).

F. 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: ð3:11Þ 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 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∶ x μ → x μ =x 2 . The Gram matrix is then defined as 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 (a) O i and O j contain equal numbers of ϕs, and (b) Δ i ¼ Δ j . These two rules are subsumed by the following much more powerful rule. Let N l ðOÞ be the number of times the l th derivative ∂ l ϕ occurs in the operator O (irrespectively of how its indices are contracted). Then the Gram matrix entry hO i jO j i can be nonzero only if 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. 12 To show this, start with ½AðxÞBðyÞ † ¼ ½AðxÞ † ½BðyÞ † , where the operators are inserted at the same radial quantization "time," so that no ordering issue arises. From here by induction in the number of fundamental fields we get ½∶AðxÞBðyÞ∶ † ¼ ∶½AðxÞ † ½BðyÞ † ∶, and Rule 3 follows by taking the coincident point limit.
Putting it all together, the Gram matrix is thus evaluated as follows. First one computes the overlaps h∂ l fμg ϕj∂ l fνg ϕi; ð3:14Þ by using the above prescription, or by using the conformal algebra, as explained e.g. in [10]. These are particular invariant tensors, symmetric and traceless in both groups of indices fμg; fνg. 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 Sec. IV C. 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 Fig. 20.

A. Simple versus generalized eigenvalue problem
Let us formalize a bit more our problem. Energy levels on the cylinder are solutions of the eigenvalue problem Hjψi ¼ Ejψi: ð4:1Þ We will be looking for scalar eigenstates, expanding them in a basis of states jji: jψi ¼ c j jji: ð4:2Þ The states jji 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: where G ik ¼ hijki is the Gram matrix discussed above. The matrix G ij is Hermitian, and H ij is Hermitian if the Hamiltonian is Hermitian, 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 transforms 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 Sec. IV C.

B. Working in presence of null states
As mentioned in Sec. III B, 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): Hjii → Hjii þ jnulli: ð4:8Þ Also, the eigenvalue problem (4.1) has to be considered a modulo addition of an arbitrary null state in the rhs. In practice, however, we will not 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 is 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 13 Strictly speaking, this is not necessary, since numerical methods for solving generalized eigenvalue problems are readily available. 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 do not 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 nongeneric 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. null eigenvalue from the UV where its value is known; one can detect it by the presence of crossings with physical states (physical eigenvalues do not 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.

C. 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 Z jxj¼1 VðxÞ O j ð0Þ; ð4:10Þ 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 jji. 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 fμg 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: By rotation invariance, the integrals will produce invariant tensors, i.e. a number of Kronecker deltas connecting the indices in fμg. Contracting these with the indices of A fμg 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 hijHjji.
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 Sec. III F, 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 perturbations-this 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 Fig. 20.

V. CUTOFF DEPENDENCE AND RENORMALIZATION 16
A. 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 powerlike. We 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: Let us now eliminate c h by using the second equation.
We get 16 The reader may wish to skip this section and come back to it while studying Secs. VI D and VII C below. 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 E 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 E ¼Ē þ hcjΔHjci; ð5:8Þ 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.

B. Computation of ΔH
To find the correction terms, we examine the matrix element of ΔH between two states i; j ∈ H l : ð5:12Þ 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 The large energy behavior of M n can then be extracted from the part of CðτÞ which is nonanalytic 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 0 . 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, nonscalar 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 Z 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 nonanalytic 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.

C. 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 belowthe 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 a RG improvement of the correction procedure. This is inspired by the usual RG in perturbative quantum field theory, which resums large logarithms. Here we do not 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 17 This term is the analogue of wave function renormalization in ordinary perturbation theory. 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 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 energyis not a fully local UV regulator. 20 So we have to learn to live with nonlocal correction 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. Second, although in principle the nonlocal 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 nonlocal 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 outthe correction procedure would break down as soon as E 0 ∼ Λ UV , as seen e.g. by the blowup 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 18 This shows once again that corrections due to integrals of nonscalar 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. 21 The second-order correction to the ground state energy is negative. Assuming that higher-order corrections do not change the situation, we may expect negative energy density at large R.
Studying many examples of RG flows known in d ¼ 2, this seems to be invariably true. The only exceptions happen when the dimension of the perturbing operator exceeds d=2. In this case the renormalized ground state energy density may be positive, although the nonrenormalized, divergent energy density is still negative. In both concrete examples of d > 2 flows studied in this work, the ground state energy density is negative at large R. 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.

D. 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 [11] (following [12,13]) but see also Appendix A for other references. In particular, Sec. III of [11] 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 Reference [11] 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 Sec. IV B of [11], 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 Sec. III of [14], 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 us 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 [14] uses the UV energy and the UV state in place ofĒ and c. 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. Fig. 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 [14] 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 [14], 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.

A. 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 Sec. III E. 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 quantized massive scalar it is 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. [15]): ð6:2Þ 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 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 Fig. 1(c). 23 Notice that this vacuum energy is not the same as the one considered in the studies of the Casimir effect, where one renormalizes by subtracting the vacuum energy density of the same theory in flat space. The latter procedure is relevant if one is interested in the dependence of the vacuum energy on the geometry keeping the mass fixed. Here we are interested in the dependence on the mass itself. 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 Sec. III E. The presence of these powerlike 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 has been observed long ago [16] that masses in an interacting theory are affected by terms which are exponentially small in the size of the torus.

B. TCSA setup
Which value of d shall we choose in our numerical study of the ϕ 2 flow? As already mentioned in Sec. II B, 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 ð6:6Þ 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 Fig. 2. The dotted line gives the number of physical states, counted using group theory. The total number of states in a d-dimensional CFT grows with Δ exponentially [15] 26 : where C is a theory dependent constant related to the prefactor in the free energy density dependence on the temperature. 27 It is not hard to see that the number of scalar states will also grow exponentially with the same exponent, although with a smaller prefactor. In Fig. 2 we can clearly In future studies one should perhaps separate the null states to speed up the numerics. 24 The two subtractions in (6.3) remove the divergences that originate from the normal ordering of the operators in the bare CFT Lagrangian and the bare ϕ 2 operator, respectively. These divergences are intrinsic to the CFT and not associated with the RG flow. 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 Sec. VC. 26 See [10] for a review. 27 C ¼ d½ζðdÞ 1=d =ðd − 1Þ 1þ1=d for a free massless scalar [15].
HOGERVORST, RYCHKOV, AND VAN REES PHYSICAL REVIEW D 91, 025005 (2015) 025005-14 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 Sec. V. As discussed in Sec. III B, 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 nonisomorphic 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).

C. 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 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 Sec. IV C, 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 2even and Z 2 -odd states. The two sectors have a 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 Fig. 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 Sec. V. As we have seen, the errors induced by omitting the states with E > Λ UV are expected to go down as a series of power law 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 Fig. 3 have been produced using such a renormalization procedure (see Sec. VI D 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 wave function to fit into the sphere. Second, it takes us much beyond the radius of convergence of conformal perturbation theory R c ¼ ðd − 2Þ=2. 28 We now turn to the excitations above the vacuum. In Fig. 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 Sec. III E. 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 Blue curves marked raw: raw TCSA results, i.e. before applying any correction. Red curves marked ren.: renormalized TCSA results; see Sec. VI D. Dashed and solid TCSA curves correspond to cutoff Δ max ¼ 12ð18Þ. 28 As determined by the leading singularity in the exact expressions for the vacuum energy density and the massive spectrum, located at m 2 R 2 ¼ −ν 2 .

TRUNCATED CONFORMAL SPACE APPROACH IN D …
PHYSICAL REVIEW D 91, 025005 (2015) 025005-15 one particle at rest, three particles at rest, then three-particle states in relative motion, etc.
In Fig. 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.

D. Renormalization details
The general method of renormalization was presented in Sec. V. 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: 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 hs. 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 Sec. V B 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.   4 (color online). A few lowest massive excitations from the raw TCSA spectra at Δ max ¼ 18 (blue dots connected with a line to guide the eye) vs exact spectrum (magenta lines). Left (right): Z 2 -even (odd) sector. The gray region indicates the sliding UV cutoff (6.7). 29 This argument shows that, more generally, corrections will vanish for the ϕ nþm and ϕ nþm−2 operators in the ϕ n × ϕ m OPE. This observation will be useful for the general Landau-Ginzburg flow in Sec. VII C.
Thus the only leading nonvanishing 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 Fig. 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 Sec. V B. 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 is given by ð6:11Þ It is 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 Fig. 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 stateindependent 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) E ren ¼Ē þc i ðΔHÞ i jc j : ð6:13Þ It is these renormalized results which were used to produce Figs. 3 and 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 : The eigenvectors are assumed normalized viac ic i ¼ 1. Of course these two eigenvectors are related, up to normalization, via the Gram matrix: As mentioned in Sec. III F, computing the full Gram matrix may be expensive, although we did find an indirect way to do it, described in Sec. IV III. 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 findsc 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 30 This also implies that the RG improvement discussed in Sec. V C is not of much use for this particular example: the mass parameter never appears on the right-hand side of the renormalization group equations (5.21), so their solution is straightforward and essentially given by (5.20), i.e. the unimproved equation. 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.

TRUNCATED CONFORMAL SPACE APPROACH IN D …
PHYSICAL REVIEW D 91, 025005 (2015) 025005-17 renormalization subtractions are suppressed by higher powers of 1=Λ UV . One last aspect we discuss here is an assumption implicit in the entire procedure of renormalization, namely that the contribution of high-energy states to low-energy observables is suppressed. It is possible to get a feel about this assumption by studying the distribution of eigenstate components in energy, defined as In Fig. 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.

A. 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 is 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 t ≡ m 2 =λ 2=ð4−dÞ : ð7:2Þ 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 jtj ≫ 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 second-order phase transition at t ¼ t c ; see Fig. 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 Fig. 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 Fig. 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.

B. 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 2even (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 Fig. 8). 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 Sec. V, with the theory specific details described in Sec. VII C 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 Fig. 9. Let us start with plots at a fixed R and varying λ. In Fig. 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  9 (color online). The range of R and λ explored in our study. Subsequent figures will show the spectrum dependence along the vertical and horizontal sections of this region, shown by the arrows. 32 It has been rigorously shown in the constructive field theory literature that for positive m 2 , perturbation theory is Borelsummable for all couplings in d ¼ 2 [17] and d ¼ 3 [18]. 33 See Sec. VII C for the RG equations for the Landau-Ginzburg flow. 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 Fig. 10.
Another interesting feature of the spectra in Fig. 10 is the absence of level crossing-eigenstates belonging to the same Z 2 sector do not cross. There are several values of λ when a pair of same Z 2 -parity eigenstates comes close to each other, but then repel. This should be 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 nonintegrable flows has long been noticed in the d ¼ 2 TCSA literature (see e.g. [20,21]), 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 Fig. 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. 35 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: 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.  34 Such tunneling effects were for example studied in TCSA in Ref. [19]. 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). 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 Fig. 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 [22]. 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 Sec. III C, 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 [22]. 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: λ c ≈ 0.55-0.56: ð7:7Þ

Spectrum for a fixed λ and varying R
We next present how the spectrum depends on R for a fixed value of λ (Fig. 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 on 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 would not 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. Focusing 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 Fig. 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

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 [16]: 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 do not 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 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 [23]. In Fig. 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 [23]. It would be interesting to work out the corresponding theory for the sphere. Finally, in Fig. 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.    16 (color online). Blue curves: E i − E 0 − 2M for the lowest two orbital excitations in the two-particle 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 fourparticle state at rest. Dashed black lines: 2ðΔM curv þ ΔM l Þ for l ¼ 1; 2. 2 × ΔM l ; As is clear from Fig. 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.
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 Sec. VII B 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 Fig. 17 the excitation energies times R. We vary λ in the range 0.53…0.56, roughly the range determined in Sec. VII B 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 [7,22]). 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 ϵ 0 , 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 abovementioned states in the spectrum in Fig. 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 do not 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 Sec. VII C).
We mention here that while TCSA should in principle be able to reproduce long-distance physics both in the gapped (massive) and the gapless (CFT) phases, it may not be the best approach from the point of view of numerical accuracy if one is interested only in the IR CFT. For example, the recently revived conformal bootstrap [28] is likely then to give better precision [see e.g. [27,[29][30][31] for the ongoing work about the Ising and OðNÞ models in d ¼ 3].
We also mention in this respect the recent work [32] describing a Monte Carlo simulation of the critical point of the three-dimensional Ising model not in the traditional R 3 geometry, but in the S 2 × R geometry, identical to the one used by us. In principle, this method could be used to simulate the full flow, not just the critical point, but one has to be careful about the approach to the continuum limit, making sure that the quartic coupling becomes small at the cutoff scale. Another issue faced by the lattice simulations on S 2 × R is that it is hard to lattice-discretize the theory on the two-sphere, because of its curvature. That is the reason why [32] uses a discretized icosahedron rather than the sphere. Preparatory work to find a true spherical discretization is ongoing [33,34], and we are looking forward to realistic simulations.
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 Sec. VII B 1-it means that the Z 2 symmetry is spontaneously broken. In the infinite R limit we expect a pair of degenerate vacua j0i L;R , each with its own tower of excitations jii 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 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. 38 It is 4 in d ¼ 2 and in d ¼ 4 − ϵ, and ≈3.83 in d ¼ 3 [24][25][26][27]. 39 It is interesting to compare the discussion below with Sec. VI of the contemporaneous work [5] devoted to Landau-Ginzburg flows in d ¼ 2; see also footnote 6 . of energies E i AE δ 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 Fig. 13-in each of the three approximately degenerate pairs, it is 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 is not what we see in the λ ¼ 0.9 plots in Fig. 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 is 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 re-expand 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 say a few words about the overall structure of the massive spectrum at large R. We identify the two neardegenerate pairs of even/odd states with two massive excitations, of mass M 1 ≈ 1.6; M 2 ≈ 2.5: ð7:13Þ Avery 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 [35]. 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 [36]. Their binding energy, exponentially small at weak coupling, is known in the leading and first subleading exponential approximation [36,37]. Apparently, here we are observing the same effect in d ¼ 2.5 dimensions and at strong coupling.

C. Renormalization details
The renormalization in the Landau-Ginzburg flow is determined through the leading OPEs of the deforming operators in (7.1), given by 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 nontrivial functions of the cutoff. We are in a position to use the RG-improved formalism of Sec. V C.
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 Eq. (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 Fig. 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 Sec. V C 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 Sec. V C, 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 Sec. V C 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 nonlocal terms given in (5.22 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 right-hand sides of the flow equations. 40 In Fig. 19 we present an example of the flow of several couplings. We see that g 2 receives substantial corrections, . Exact values are given in blue (isolated dots), a moving average in red (dashed line), and our continuum approximation in black (smooth curve). 40 For g ðH CFT Þ 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 wave function renormalization in perturbative RG. We implemented it, but found the numerical effect of this improvement to be very small. 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.
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: (i) ΔH loc which reflects the change in all local couplings; (ii) ΔH 1 and ΔH 2 which include the nonlocal terms proportional to H CFT :V c þ V c :H CFT and ðH − E r Þ:V c þ V c :ðH − E r Þ, respectively; (iii) ΔH nl . It is now time for numerical diagonalization. We add ΔH loc and ΔH 1 directly to the bare TCSA Hamiltonian, and diagonalize. Let us 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 is 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 is 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 Sec. VI D.
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: ðΔE n Þ nl ¼ c n :ΔH nl ðE n Þ:c n ¼ ðc n Þ i ðΔH nl ðE n ÞÞ i j ðc n Þ j : ð7:20Þ 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 Sec. VII B. 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 nonnegligible 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.

D. Nonunitarity and complex energy levels
As we observed in Sec. III C, 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 Fig. 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 negativenorm 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 Sec. III E, 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 negativenorm 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 are 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 seven negative-norm states. We find that the matrix of the ϕ 4 perturbation within this subspace has one pair of complex conjugate eigenvalues ≈1.85 AE 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 FIG. 21 (color online). The spectrum around Δ ¼ 11.5 for m 2 ¼ 1, λ ¼ 0.55, and 0 < R < 0.15 with a step of 10 −4 . 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 ¼ 12. 41 In principle, numerical spectrum for a finite Λ UV could contain eigenvalues with small imaginary parts, approaching zero as Λ UV → ∞. However, in our numerical studies at d ¼ 2.5 we observed that even the truncated spectrum was real for the purely massive perturbation. 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 [38].
Going back to d ¼ 2.5, the above argument is confirmed numerically in Fig. 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. Did we not say that the Landau-Ginzburg flow is not integrable, and that in nonintegrable flows energy levels do not 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-level-crossing rule does not apply. To see this, consider a toy-model 2 × 2 symmetric generalized eigenvalue problem where σ ¼ AE1 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 Fig. 22 we show the spectrum with RE 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 FIG. 22 (color online). The spectrum at RE 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. conjugate pair, which sometimes re-emerges as a pair of real eigenvalues for a slightly larger value of R. The most prominent such collision happens at R ≈ 0.7.
To resolve the multitude of eigenvalue curves in Figs. 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 do not 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 Sec. VI, whose extended Hilbert state included some null states, but only as a matter of convenience. In a nonunitary 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 nonunitary 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 negativenorm 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.

VIII. 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 Figs. 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 do we 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 three-dimensional gauge theories. The easiest problem in this class might be three-dimensional QEDthe 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 non-Abelian gauge theories in d ¼ 4. 42 42 The interaction in this case is only marginally relevant. In principle, nothing prevents using marginally relevant perturbations in TCSA, once renormalization issues are understood. In the d ¼ 2 case such flows were studied in [39].
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 Sec. VII B 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 nonsupersymmetric 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.
integrable SUSY sine-Gordon model was studied using TCSA methods in [45].
We also mention a nice paper [46] (see also [47]), 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,20], one also separated the states into quasiprimaries 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 ¼ hijji 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, hijϕð1Þjji. Here one uses the commutation relations with a Virasoro primary ½L n ; ϕðz;zÞ ¼ z n ½ðn þ 1Þh þ z∂ z ϕðz;zÞ; ðA2Þ 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 [20]. 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 [48] (described in more detail in [21]) 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 [12,13] (see also [49]). 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 [11] and more recently in [14] (see also [39], which in particular includes a study of marginally relevant perturbations). Our discussion of renormalization in Sec. V was inspired by [11], although, as we discussed in Sec. V D, our method is different in several aspects from [11,14].

APPENDIX B: OTHER HAMILTONIAN TRUNCATION TECHNIQUES
We 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 [50], which studied a two-dimensional 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 [51,52] (see also [53]) this approach was applied to the massive ϕ 4 theory in d ¼ 2 and d ¼ 3 dimensions. The reader should compare the spectra in Fig. 8 of [51] (d ¼ 2) and Fig. 1 of [52] (d ¼ 3) with our Fig. 10 in d ¼ 2.5. We believe that it is worth revisiting the approach of [51,52] and possibly improve on some of the implementation details. 44 As we mentioned in Sec. III E we are planning to do this at least in d ¼ 2 where TCSA is not directly applicable [9]. 44 For example, their "quasisparse eigenvector method" does not seem to take into account the high density of states present at high energy, which may compensate the smallness of individual contributions of each one of these states. Their "stochastic error correction" procedure [52] computes the sums of squares of highenergy matrix elements of the type encountered in Eq. (5.12) via a Monte Carlo procedure. As we observed in Sec. VII one should be able to compute such corrections analytically.