Truncated Hilbert space approach to the 2d ϕ4 theory

We apply the massive analogue of the truncated conformal space approach to study the two dimensional ϕ4 theory in finite volume. We focus on the broken phase and determine the finite size spectrum of the model numerically. We interpret the results in terms of the Bethe-Yang spectrum, from which we extract the infinite volume masses and scattering matrices for various couplings. We compare these results against semiclassical analysis and perturbation theory. We also analyze the critical point of the model and confirm that it is in the Ising universality class.


JHEP10(2016)050 1 Introduction
Interactive quantum field theories (QFT) cannot be solved exactly. There are some exceptions including two dimensional integrable QFTs, which can be investigated by the bootstrap method, but these theories are very special, as the scattering matrix is factorized and there is no particle production or decay. To study more realistic QFTs, one has to rely on alternative approximative methods such as perturbation theory, variational methods, lattice discretizations or the truncated conformal space approach (TCSA) just to mention a few.
The truncated conformal space approach was originally developed by Yurov and Zamolodchikov to study the relevant perturbation of the Lee-Yang conformal field theory [2]. The idea of the method is to put the system in a finite volume, which makes the spectrum discrete, then one can truncate the Hilbert space by keeping states below some energy cut and diagonalize the truncated Hamiltonian numerically to get an approximate spectrum as the function of the volume. This is basically the usual variational method of quantum mechanics with the unperturbed energy levels as basis states and the obtained results are not perturbative.
The TCSA method has been extended to analyze the perturbation of other conformal field theories [3][4][5], too. The advantage of the method is that it can be applied for non-integrable perturbations [6]. Indeed it has been successfully applied to describe such phenomena as phase transitions, decay rates, form factors or boundary entropies [6][7][8][9]. It also has a natural extension for models with boundaries and with a defect [10,11]. Recently the method has been extended even to higher dimensional theories [12].
The simplest non-integrable interacting quantum field theory is probably the two dimensional φ 4 theory. The TCSA method was used in [13] to study the spectrum by compactifying the zero mode of the field. This approach gives a good result for small couplings. Alternatively, in [1,14] the authors used the finite volume massive basis to build up the truncated Hilbert space and analyzed both the unbroken and the broken phases of the model for a wide range of couplings. There special emphasis was put on developing the method and on analyzing carefully the convergence properties in eliminating the truncation. The authors also searched for the critical point, beyond which the parity symmetry is broken and checked Chang duality, by which one can realize the same system with two different parameter sets, which are related to each other in a strong-weak fashion. The truncation dependence of the spectrum was further analyzed in [15]. The aim of our paper is to use the same truncation of the finite volume massive basis for the broken phase and obtain results for a wide range of couplings similarly to [1,14]. We focus on the volume dependence of the spectrum and try to extract the infinite volume parameters from finite size effects. This analysis requires precise data, which we achieved, first, by implementing an efficient C++ program to set up and diagonalize large sparse matrices, and second, by an extrapolation in the truncation level. We compare the extrapolated finite size spectrum with semiclasssical results, perturbative results and also with the Bethe-Yang (BY) spectrum. This BY spectrum incorporates the leading finite size corrections based on the infinite volume masses and scattering data and enables one to determine these infinite volume characteristics at moderate volumes when the truncation errors can be controlled.

JHEP10(2016)050
The paper is organized as follows: in section 2 we review the classical and semiclassical results. Section 3 contains the definition of the quantum model. In section 4 we explain our truncated Hilbert space method. Theoretical predictions are summarized in section 5, while numerical results are presented in section 6 and compared to both the semiclassical, the perturbative and to the Bethe-Yang spectrum. From this comparison we extract the infinite volume mass spectrum and scattering phases. We also analyze the spectrum near the critical point. We close the paper with conclusions in section 7. Appendix A explains volume dependent normal ordering and its relation to the ground-state energy. Appendix B introduces the parametrization of the scattering matrix from unitarity relations, while appendix C explains the numerical implementations.

Classical, semiclassical and perturbative results
In this section we recall the available results for our theory We follow [13], but see also [16,17] for more details. Sometimes we simplify formulas by introducing q as λ = 6q 2 . We take = c = 1. As in two dimensions the scalar field is dimensionless, both m and q have dimension of energy, i.e. 1. The parameter m can be chosen to determine the mass scale and then the physics depends on the dimensionless ratio: m q . Observe that the coefficient of φ 2 has a "wrong" sign, such that we analyze the theory when the Z 2 parity symmetry will be spontaneously broken.

Classical analysis
The potential has two minima at ± m q and its expansions around the minima reads as The mass of the small fluctuations is m 0 = √ 2m and the minimum of the potential is V 0 = − m 4 4q 2 . The two static solutions, called the kink and the antikink, interpolate between the vacua as and have masses The breather-like bound state of the kink and the antikink is not a stable solution as it radiates very slowly [18].

Semiclassical considerations
We can quantize the theory semiclassically. For this normal ordering is needed, which has to be defined w.r.t. a free theory. We recall for later reference that in the calculations [16,17] the reference theory was chosen to be the infinite volume free massive theory describing the fluctuations around any of the minima with mass m 0 . The relevant quantities at the quantum level are the masses and the scattering matrix.
The semiclassical corrections to the kink masses read as where, as we pointed out, the parameters m, q are understood as the coefficients of the normal ordered perturbations w.r.t. m 0 . The breathers can be considered as kink-antikink bound states with masses [19] m n = 2M sin n ζ 2 (2.6) where Although classically there is no stable breather solution due to the radiation into the small fluctuations of the field, at the quantum level the small fluctuations are identified with the breather itself thus can be stable. Higher breathers are bound states of the lightest one The existence of the n th bound state requires the scattering pole nζ to be in the physical strip: nζ < π, i.e. ζ < π n . The stability against the lightest particle further requires that m n < 2m 1 , which is always satisfied for m 2 but never for m 3 and for the higher ones. It was conjectured in [19] that this picture survives at the quantum level in the sense that we cannot have more than two breathers. Let us anticipate that in the deep quantum domain even the first breather is unstable against decaying into a kink and antikink pair.
The scattering matrix for the lightest breather satisfies unitarity and crossing symmetry below the two-particle threshold Here we introduced the rapidity variable, see appendix B for details. The simplest nontrivial solution of (2.9) is It implies a bound state with mass m 2 = 2m 1 cos κ 2 . Comparing this to the semiclassical mass formula (2.8), we conclude that κ = ζ.

JHEP10(2016)050
From the analysis of the small fluctuations over the kink solution the classical limit of the breather-kink scattering matrix can also be extracted [16]: The pole at θ = iπ/2 is related to the translational mode, while the pole at θ = iπ/6 signals an excited kink state with mass

Related integrable models and perturbative results
There are two integrable models, which are related to our theory. The sinh-Gordon theory is defined by the Lagrangian It contains only one type of particle with scattering matrix 14) The sine-Gordon theory can be obtained by an analytical continuation b → iβ, which admits kink and anti-kink solutions, too. The other interesting theory is the Bullough-Dodd theory Combining the perturbative expansions of the two models one can extract the leading order perturbative S-matrix of our model: Direct perturbative calculation for the mass correction gives [1]: We would like to compare all of the results of this section against the finite size spectrum we obtain from our truncated Hilbert space approach.

JHEP10(2016)050 3 Quantum theory
At the quantum level we analyze the spectrum of the Hamiltonian in a finite volume L: where φ(x) = ±φ(x + L), i.e. we consider both periodic (+) and antiperiodic (−) boundary conditions in order to accommodate one-particle kink states. As operators are multiplied in the same position we need to regulate them by normal ordering. Normal ordering is performed with respect to some prescribed mass scale µ and volume L, thus the Hamiltonian density takes the form: Here the field is expressed in terms of creation-annihilation operators as where [a n , a m ] = [a + n , a + m ] = 0 ; [a n , a + m ] = δ n,m (3.4) and The index n takes integer values for periodic, while half-integer values for antiperiodic boundary conditions. Normal ordering, denoted above by :: µ,L , is understood that all creation operators a + n are on the left of the annihilation operators a n . The vacuum is chosen to be annihilated by all a n . The free Hamiltonian is where the renormalized ground-state energy is chosen to be for periodic/antiperiodic boundary conditions. Since both the ground-state energy and the ground-state energy density is infinite for the non-renormalized Hamiltonian, we have to introduce proper renormalization. In appendix A, we show how E ± 0 (µ, L) appears from comparing different renormalization schemes. In our scheme the vacuum energy density is normalized to zero and the ground-state energy is chosen to vanish for infinite volume.
To relate (3.2) to another normal ordering at scale µ and volume L , we recall that [13,14] : φ 2 : µ ,L =: φ 2 : µ,L +∆Z ; Z(µ, L) itself is only meaningful with a UV cutoff N , but in the difference, ∆Z, the limit N → ∞ can be taken. Similarly : φ 4 : µ ,L =: φ 4 : µ,L +6∆Z : φ 2 : µ,L +3(∆Z) 2 (3.10) As we already mentioned, the scheme in which the semiclassical calculations were done is defined by normal ordering at L = ∞ and µ = m 0 . The parameters are chosen there as To relate this scheme to the one with mass m and L = ∞, we calculate such that the coefficients of the same physical system in this scheme are and As we work in the finite volume scheme we introduce [14]: For periodic boundary condition we have while for antiperiodic boundary condition we found The coefficient at volume L of our theory are given by: The correction to the vacuum energy density due to normal ordering is Since z ± (L) is exponentially small, it can be often neglected.

JHEP10(2016)050
It is tempting to introduce the normal ordering at µ = 0 and L = ∞, as then the scale would not appear as a dimensionful parameter [13]. Choosing µ = 0 implies, however, that Z(0, ∞) − Z(m 0 , ∞) is IR divergent, which can be regulated by an IR cutoff a, but then the theory will depend on the regulator a. Naively we could just use Z(0, L) (conformal scheme) instead of introducing Z(0, ∞) but we should keep in mind that L → ∞ is not a safe limit as it basically puts a to zero. As a result, the mass of the breather will diverge logarithmically in the L → ∞ limit in this conformal scheme.

Truncated Hilbert space approach
In this section we explain how we diagonalize the finite volume Hamiltonian. For numerical investigations we introduce dimensionless quantities such that the dimensionless Hamiltonian can be written as andẽ ± 0 (l) = E ± 0 (mL)/m. From now on all normal orderings are evaluated at scale m and volume L, unless otherwise indicated.
Technically we use the Hilbert space built up by the oscillators of the free Hamiltonian in volume L and diagonalize h at a certain truncation level. Clearly we have two dimensionless parameters: the dimensionless volume l and the quartic coupling parameter g. The dimensionful mass m sets purely the energy scale.

Mini Hilbert space approach
For periodic boundary condition we found it useful to adapt the "mini-superspace" approach. This is the method which was used many times in analyzing the Liouville theory and deals with the zero mode of the field separately [21][22][23]. The same approach was used also in [1] to handle the zero mode. Here we also found that the IR mode generated by a 0 and a + 0 plays a very special role in forming the low lying states in the spectrum. We thus decompose the elementary field as

Technical implementations
For periodic boundary condition we split the Hilbert space to the mini-Hilbert space containing the zero mode of the field and to the rest: The mini-Hilbert space H mini is generated from the vacuum by acting with a + 0 , whileH by acting with all other creation operators a + n =0 . Technically we determine the spectrum of h mini first. For this we choose 2000 eigenstates of h mini 0 and diagonalize h mini on these states. We found it enough to use only the first 6 states to build up the relevant product basis for (4.7).
For antiperiodic boundary condition we do not have the analogue of the mini Hilbert space, so we diagonalize the Hamiltonian on the full truncated Hilbert space.
We want to keep the dimension of the truncated space the same for each volume. Since by changing the volume the free energy levels are changing, too, we decided to generate the basis for the Hilbert space at mL = 1. At such a small volume the spectrum is close to the conformal spectrum and the energy can be approximately measured in units of 2π mL . Consequently, the mode (a + n ) jn contributes to this "conformal energy" as j n |n|. We introduce integer cuts in this unit, which we increase two by two (one by one in the antiperiodic sector), usually between 14 and 26. With each cut we keep the same oscillator content for any dimensionless volume, which typically takes integer values up to l = 30.
There are conserved charges in the theory, which help to decrease the dimension of the Hamiltonian below the energy cut in each sectors. Momentum is a conserved quantity and we focus through the whole paper on the zero-momentum sector. Field parity φ ↔ −φ is also a symmetry of the full theory. As a consequence we can diagonalize the Hamiltonian separately on the even and odd particle number sectors. Additionally, there is also space parity x ↔ −x, but in this paper we will not exploit this symmetry.

Description of the spectrum
According to Lüscher, the finite size spectrum of a quantum field theory can be described in terms of the infinite volume data [24,25]. These data consist of the masses of the excitations, their dispersion relation and the various scattering matrices. The leading finite size correction comes from momentum quantization, which is described by Bethe-Yang type equations and expresses the periodicity of the multiparticle wave function [25]. These corrections are polynomial in the inverse power of the volume. There are exponentially JHEP10(2016)050 small corrections, too, and they are related to vacuum polarization effects, or in some other words, virtual particles [24,26].

Infinite volume data
In infinite volume, the theory is parametrized by the masses of the excitations and their scattering matrices. We use unitarity and consistency with a generalized bootstrap idea to restrict the form of the scattering matrices and masses. The bootstrap idea states that the scattering matrix of the bound state is the product of the scattering matrices of its constituents. It is definitely satisfied in integrable theories and we expect that it is also satisfied for nonintegrable theories below the inelastic threshold.
Let us start with the breathers. According to the semiclassical analysis, we can have at most two breathers. Let us denote the lighter one by B 1 and its mass by m 1 , while the heavier is denoted by B 2 , whose mass is m 2 . As we are in a relativistic theory their dispersion relation is This motivates us to introduce the rapidity parameter: Lorentz transformation simply shifts the rapidity and, as a consequence, the scattering matrix depends only on their differences. Since B 2 is the bound state of two B 1 s, the B 1 B 1 → B 1 B 1 two-particle scattering matrix, which we denote simply by S 11 (θ), must have a pole whenever B 2 is in the spectrum. In appendix B we derive a parametrization of the scattering matrix with the existence of the required pole from unitarity: and the position of the pole iα and the density ρ 11 (θ) is an unknown function of the dimensionless parameter g. The threshold rapidity θ t corresponds to the kinematical point where the Semiclassically, for small g, according to (2.7), the parameter α starts as In the analysis we will assume that θ t is large and the exponential piece can be neglected. It would be very interesting, although not easy, to measure this piece and compare it to other scattering matrix elements. Since the scattering matrix has a pole at iα: Figure 1. The kink-anti-kink B 1 fusion process, together with the corresponding pole in the kink B 1 scattering process. the mass of the bound state m 2 can be expressed as Below the θ t threshold the scattering matrix behaves as if it were in an integrable theory. As a consequence one can try to use bootstrap to calculate the scattering matrix of however, this scattering could be intercepted by the B 2 B 1 → B 1 B 1 process, which is taken into account partially byS exp 21 (θ). Additionally the appearing pole at i 3α 2 should correspond to a resonance, since B 3 is not a stable particle, i.e. it should not lie on the imaginary line. All together the correct parametrization possibly is Let us focus on the kink and the antikink now. The Z 2 symmetry implies that these particles are degenerate in their mass, what we denote by M . Interestingly the kink-antikink scattering matrix is also diagonal and satisfies crossing symmetry and unitarity below the threshold [19]: where in the last line we used charge conjugation symmetry. A kink-antikink pair can fuse either to B 1 or to B 2 . The corresponding angles are β 1 and β 2 , such that the masses are m n = 2M cos βn 2 . The fusion process is shown on figure 1. This further implies that where we included S exp kk (θ) for later consistency, which has a similar form as (5.3). Here it is assumed that we are in the domain where the kink is the lightest particle. Note that kinematically it happens in a different domain of the coupling g , where we can trust in (5.3) without the exponential piece.

JHEP10(2016)050
The kink-B 1 sector has a topological charge which separates it from the neutral sector. Thus, even for m 1 < M , the scattering matrix S k1 also satisfies unitarity and crossing symmetry below the kink-B 2 threshold where again we used parity and charge conjugation symmetry. Due to the kink-anti-kink fusion process this scattering matrix must have a pole at i β 1 2 as demonstrated on figure 1. Taking into account that in the classical limit (2.11) we have a pole at θ = iπ/6 , which could get quantum correction γ, we arrive at the parametrization Demanding that fusing a kink with an antikink we obtain the breather-breather S-matrix the existence of the factor S(θ, iα)implies that β 1 = π − α (for α < π 2 ). This actually gives the expected mass relation m 1 = 2M sin α 2 and is also consistent with the fact that in the classical limit (2.11) β 1 /2 should be π/2.
Interestingly the fusion also gives two additional factors which indicates a pole in the physical strip at θ = 2π 3 − α 2 + γ. As this pole is already above the θ t threshold we shouldn't take it very seriously. Nevertheless, if we would take it, the only interpretation could be a self-fusing pole, which is also forced by having a cubic vertex in the perturbative definition (2.2). This implies the relation γ = α 2 and leads to the scattering matrix: This formal scattering matrix is very similar to the Bullough-Dodd scattering matrix (2.16) and satisfies the bootstrap equation S 11 (θ + i π 3 )S 11 (θ − i π 3 ) = S 11 (θ). Comparing, however, to the perturbative expansion (2.17) we can see a mismatch, thus we restrict our analysis only to the well-founded form (5.3). The γ = α 2 choice would lead to the following mass of the quantum excited kink: Using that m 2 = 2m 1 cos α 2 = 2M sin α , we can see that β 2 = π − 2α. From the fusion we also obtain Our analysis is very natural from an integrable point of view, moreover, a similar parametrization comes out from semiclassical investigations, too [19]. In the following we use these masses and scattering matrices to describe the finite size spectrum, which we compare to the numerical data. This will decide whether such bootstrap based proposals are adequate or not.

Polynomial finite size corrections
The leading finite size corrections are polynomial in the inverse of the volume and are due to momentum quantization. If the volume is large compared to the particle number then particles travel freely most of the time, except when they scatter on each other such that their wave functions pick up the scattering phases. Summing up all these phases when they move around the circle and adding the translational phase, the result should be an integer/half-integer multiple of 2π, ensuring that the wave function is periodic/antiperiodic: This equation is called the Bethe-Yang equation and is valid for integrable models, where particle numbers cannot change. We believe that in certain kinematical domains it gives a good description of the spectrum even for non-integrable models. Once the momenta are calculated, the energy is simply In the following, we elaborate this description for small particle numbers. For a one-particle state it implies that e im j L sinh θ = 1 ; j = 1, 2 (5.21) Taking its logarithm, we can see how the momentum is quantized in a volume L: This is identical to the free momentum quantization, thus the one-particle levels get only exponentially small, interaction-dependent finite size corrections. In the zero-momentum sector with periodic boundary condition we have zero momentum B 1 and B 2 particles, i.e. with n = 0. For antiperiodic boundary condition we have zero momentum kink states. For a state with two B 1 particles having opposite momenta p = ±m 1 sinh θ, the equation implies that e im 1 L sinh θ S 11 (2θ) = 1 (5.23)

JHEP10(2016)050
Taking its logarithm we arrive at Since the θ → −θ change brings (5.24) to the same equation but for n → −n, we assume that n ≥ 0. Once this equation is solved at a given L and n for θ, denoted by θ n (L), the energy of the corresponding two-particle state is In theories where there is a bound state pole in the scattering matrix, the equation (5.24) admits a complex solution. Indeed, let us assume that in eq. (5.23) the rapidity θ can be written as θ = iu with u being real. Clearly the factor e −m 1 L sin u goes to zero in the large volume limit. As the product is 1, the S-matrix contributions should diverge, which forces 2u = α + δu(L), where δu(L) goes to zero in the large L limit. We can even calculate how it vanishes by expanding eq. (5.23): The energy of this state is Clearly it describes a zero momentum B 2 particle together with its leading exponential finite size corrections (5.31). The equation (5.24) sums up all of these exponentially small corrections for n = 0. Similar analysis can be performed in the antiperiodic sector using the Bethe-Yang equations of the kink-breather scatterings. There the breathers are quantized as

Leading exponential corrections
The leading finite size correction for the vacuum is described by where m is the mass of the lightest particle in the spectrum. This correction is related to the virtual process in which a particle-antiparticle pair appears from the vacuum, travel around the world and annihilate when they meet again. The general form of the leading mass corrections was calculated by Lüscher [24]. It has two terms: the F -term corresponds to the process in which a virtual particle-antiparticle pair of type b appears from the vacuum and after traveling around the world and scattering on particle type a, they annihilate. The (−1) is the subtraction of a similar process for the vacuum. The µ-term is the residue of the F -term at the bound state pole: where µ c ab is the altitude of the mass triangle with base m c , and it is being understood that c appears as a bound state of a and b. Physically it corresponds to the process in which a particle of type a decays virtually into the particles b and c, which both travel forward in time. After traveling around the finite world they meet again on the other side and fuse back to particle a.

Numerical analysis and interpretations
In this section we summarize our numerical results. The paper [13] used the conformal normal ordering scheme. To test our numerical implementations, we first developed a program in this convention and successfully reproduced their weak coupling data. Then, to avoid infrared divergences for large volumes, we switched to the (m, L) normal ordering conventions (3.18). All of the data in this section are in this (m, L) convention. We plot the dimensionless energies E/m against the dimensionless volume l = mL and we analyze both the periodic and the antiperiodic sector.

Breather sector
We first show a typical spectrum at g = 1.2 on figure 3. Clearly the system has a ground state energy density such that for large volumes all energy levels contain a term E bulk L. This term can be analyzed already on the ground state, which goes as E 0 (L) = E bulk L + E exp (L). Here E exp (L) is exponentially small for large L. We normalized the energy density to zero for the free boson with mass m. But once the interaction is switched on, the ground state energy density becomes a measurable physical quantity.
In measuring every quantity we should keep in mind that we determine the spectrum at various truncation levels e cut . The dependence of the energy density on the cut-off is demonstrated on figure 4.
As truncation effects are very large, we have to extrapolate in the cut-off dependence in order to get a reasonable result. Based on the analysis in [14] we found it useful to fit the following function E(e cut ) = a + b log 2 e cut e 2 cut + c log e cut e 2 cut (6.1) In order to say something about the accuracy of the extrapolation, we also fitted a simpler form and used the difference of the two results for the error approximation on figure 4. We show on figure 5 the convergence of the energy density values along with the fitted function.
In the following we use (6.1) to extrapolate each energy level and show on the plots the extrapolated value only. In most of the cases we subtract the vacuum energy density, such that a typical spectrum looks like the one on figure 6.
We can easily identify the first and second breathers and even a resonance around 3m 1 . Let us start with the zero momentum bound-state particles. For the B 1 particle the energy is a constant function of the volume up to exponential finite size corrections,  however for B 2 we observed stronger finite size corrections. We blow up this part of the spectrum on figure 7. The reason is that B 2 is a weakly coupled bound state of two B 1 particles. In decreasing the volume we pump energy into the system, which makes the bound state dissolve and behave more like two B 1 states. That is why the equation (5.24) with imaginary rapidities, shown by a red solid line, provides a very precise description. This fitting is very sensitive to α, which turns out to be α = 0.06.
The other lines in the spectrum are described by the Bethe-Yang equations with real rapidities and quantization numbers n = 1, 2, 3. These curves are shown on figure 8.   Figure 8. Low lying even spectrum at g = 0.06: the first six energy levels (with the vacuum energy density subtracted) are plotted against the volume together with the solid Bethe-Yang lines for n = 0, 1, 2, 3. The resonance is well described by a three-particle imaginary BY line including a stationary breather.
The spectrum for small g can be very precisely described in a similar manner in terms of m 1 and α. On figure 9 we show a similar plot at g = 0.6.
Performing a similar analysis we demonstrated for g = 0.06 we can extract the coupling dependence of m 1 , m 2 and α. This is shown on figure 10.
We can compare the mass of the elementary excitation to perturbation theory [1]: Observe that the perturbative formula starts at linear order in g and gives a very good approximation, opposed to the semiclassical formula, which starts at quadratic order. Having a closer look at this quantum spectrum, we can see the following: for small couplings,  Figure 9. Low lying spectrum at g = 0.6: the first six even energy levels (with the vacuum energy density subtracted) are plotted against the volume together with the solid Bethe-Yang lines for n = 0, 1, 2, 3 and the three-particle BY line for the resonance. we have two breathers and a very heavy pair of kinks. As the coupling increases, the measured masses are smaller than the semiclassical values and the second breather disappears from the spectrum when it can decay into a kink-antikink pair, i.e. when m 2 (g) = 2M (g), this happens around g = 0.9. Similarly, the light breather, the elementary excitation of the quantum field disappears at around g = 1.5, where its mass equals twice the mass of the kink: m 1 (g) = 2M (g). For small couplings, the semiclassical calculations are reliable, which exclude the existence of more neutral particles. Indeed we can only see a resonance around 3m 1 . This confirms the conjecture of [19].  Figure 11. Low lying spectrum at g = 1.8: the low lying even (blue) and odd (red) energy levels (with the vacuum energy density subtracted) are plotted against the volume. Four-particle states are drawn with thick lines.
It is interesting to take a look at the spectrum at stronger couplings, when the breathers disappear from the spectrum.
Such a spectrum is shown on figure 11 for g = 1.8. Clearly we can see two and four particle scattering states. By identifying the points which correspond to the two particle Bethe-Yang lines with quantization numbers n and using (5.23), (5.25) we can measure the S kk scattering phase shift δ(θ) as follows: we first pick up a quantization number n and analyze the E 2 (n, L) curve. For a given L we transform the energy using (5.25) to rapidity E 2 (n, L) = 2M cosh (θ n (L)) −→ θ n (L) = cosh −1 E 2 (n, L) 2M (6.4) We then use (5.23) to transform the volume for the given θ to the phase shift as: Parameterizing the kink-anti-kink scattering matrix as S kk (θ) = −e iδ(θ) (6.6) the resulting phase-shifts are displayed for g = 1.8 and g = 2.4 in figure 12.
Observe that the phase shift is smaller and smaller once we increase the coupling and approach the critical point. We can fit a third order polynomial for the phase shift at g = 1.8 as follows: δ(2θ) = 2.82(3) θ − 0.77(10) θ 2 + 0.27(9) θ 3 (6.7) Clearly in the method above we needed the precise value of the kink mass. We measured this mass with two different methods: first from the splitting of the ground-state energy in subsection 6.1.2, then from the ground-state of the anti-periodic sector in section 6.2.
Both measurements indicated that the kink mass is close to its semiclassical value in the  Figure 12. Kink-anti-kink scattering phase shifts as functions of the rapidity for g = 1.8 (left) and for g = 2.4 (right). The various points were extracted from lines with different n in (5.24). For g = 1.8 we also presented a fitted third order polynomial in θ, while for g = 2.4 the errors of our data. whole range below the critical point, but the errors become large at stronger couplings. Therefore we used the semiclassical kink mass for these measurements. Since we do not consider the error of the kink mass, the errors are somewhat underestimated. So far we checked only the polynomial finite size corrections. Let us investigate now the leading exponential corrections. For g = 0.06, the mass of the lightest particle is m 1 = 1.40182 (10), which can be used to check the leading exponential correction of the vacuum energy (5.29). The comparison is shown on figure 13, providing a convincing evidence.
The similar comparison between the energy level of the first breather and (5.30) indicates that we are missing some interesting physics.
One may wonder if summing up all higher corrections coming from the elastic scattering matrix S 11 (θ) = S(θ, iα) together with the contributions of the second breather can save the situation. To test this idea, we can compare the spectrum of our model to the spectrum of the sine-Gordon theory, at the coupling where the first breather's scattering matrix agrees with S 11 (θ). We used the TCSA program from [27] and rescaled the energies and volume  such that the dimensionless quantities agree with this paper's convention. (In the TCSA program the soliton mass is normalized to 1.) The two spectra are plotted on the same figure on figure 15.
Clearly the spectra agree very precisely for smaller volumes, when the error of the TCSA analysis is small. (We did not perform any extrapolation in the TCSA energy cut which was chosen to be 20.) One surprising outcome is that the leading Lüscher correction (5.30) agrees very well with the TCSA data, but not with the φ 4 ones. It indicates that the inelastic processes play an important role in determining the finite size effects.  Figure 16. Difference of the ground-state energies at g = 0.6 of the even and the odd sector.

Kink mass from the vacuum splitting
So far we analyzed the even energy levels only. To gain information about the kink mass, one can analyze the gap between the even and odd ground-state energies [1,27]: The splitting is demonstrated on figure 16. In order to extract the kink mass, we first multiply the numerical data by √ L, then we measure the slope of the logarithm of the resulting curve by numerical differentiation. The result is shown on figure 17.
Using this analysis, we can obtain the kink mass as a function of the coupling. The curve on figure 17 seems to be monotonic. However, it actually has a little local minimum and a local maximum which can be used to estimate the error. The measured kink masses for different couplings are shown on figure 18.

Antiperiodic boundary condition
The lowest lying state in the antiperiodic sector is the kink state. After subtracting the periodic vacuum, we found the result presented on figure 19.
The low lying spectrum in the kink-sector is presented on figure 20. From the antiperiodic sector's ground state we can extract the kink mass. It is in complete agreement with the one we obtained from the ground state splitting and both are very close to the semiclassical result. These data are shown on figure 21.

Critical point and the spectrum of the Ising model
The vanishing of the kink mass indicates that there is a second order phase transition around g = 3. However, increasing the coupling above g = 1.8 renders both direct kink    Figure 20. The low lying spectrum in the antiperiodic sector at g = 0.6. We indicated the infinite volume masses of the semiclassical kink and excited kink by horizontal green and blue lines, respectively. mass measurement methods inaccurate. On the other hand, it is known that at the critical point the spectrum must be governed by a conformal field theory, which has the spectrum Here L 0 andL 0 are the Virasoro generators, which measure the left and right conformal dimensions, respectively. Their sum is the scaling dimension, while their difference is the conformal spin, which is proportional to the momentum. As we analyze the zero momentum sector of the theory we list the low lying scaling dimensions of the periodic Ising model in this sector in table 1. We compare this spectrum with l 2π times our spectrum.      figure 22. This indicates that the critical point is between g = 3.1 and g = 3.3.
We plotted the low-lying periodic spectrum at g = 3.24 in figure 23. Similarly to the critical behaviour of the two-frequency sine-Gordon model [6], we nicely recover the conformal Ising spectrum. Higher energy states approach their asymptotic values later, but all the multiplicities are correct. This result is a very strong confirmation that the critical behaviour is controlled by the conformal Ising model.

Comparison to [1]
In this subsection we compare our methods and results to those in [1]. We focus on the periodic sector only as [1] analyzed this sector. In order to facilitate the comparison we recall that the quartic coupling in [1] is denoted by g 4 ≡ g [14] , while the quadratic coupling is 1 4 M 2 [14] . Since both papers use the same scheme, it implies that M [14] = m 0 = √ 2m and g = 12 g [14] M 2 [14] . Consequently the dimensionless volumes are related as l [14] = M [14] L = √ 2l, while the energies as H m = √ 2 H [14] M [14] . In the periodical case both papers used the same mini-Hilbert space approach, but the numerical implementations of improving the cut-off dependence were different. The paper [1] used a smaller Hilbert space of size 10000 but decreased the truncation errors drastically by deriving an effective Hamiltonian, which took into account the cut-out Hilbert space in the leading non-vanishing order of the perturbation theory. (This was further improved in [15]). Contrary, we used efficient diagonalizing routines for large sparse matrices having up to 250 million nonzero matrix elements, and extrapolated in the cut dependence following their previous work [14]. We tested this method in the symmetric sector by comparing our results to the one in [14]: we used the code in [14] to produce raw and renormalized data, which we compared to our raw and extrapolated data. At g 4 = 1.5, m = 1, using an energy cut e cut = 16, the raw data agreed with high precision: the difference appearing in the fourth-fifth digit for the first eigenvalues (the truncation was slightly differently defined, so deviations were larger than the numerical precision). Our extrapolated data agreed with the renormalized data within their errors for all observables, except for the ground-state en-JHEP10(2016)050 ergy density, for which our data was a bit lower. Thus our method possibly slightly overextrapolates. For this reason, in the future we are planning to combine the two methods.
In comparing the results we note that we managed to analyze the model in a wider parameter range than [1]. For the overlapping parameters we found agreement within the numerical errors. For large couplings we identified the critical point g = 3.2 (1), which agrees with their findings using Chang duality. Since the numerical errors grow with the volume we put emphasis on the finite size spectrum including all polynomial and the leading exponential corrections. We did it in all sectors of the theory. In measuring the kink mass via the vacuum splitting we fitted their finite size formula. Having a good control of the finite size effects enabled us to extract physical information from the volume range l = [5,15] where both the exponential theoretical corrections and the numerical errors were small.

Conclusions
In this paper we investigated the two dimensional φ 4 theory in the broken phase by the massive generalization of the truncated conformal space approach. We called this method the truncated Hilbert space approach and we developed a code based on the diagonalization of large sparse matrices to diagonalize the truncated dimensionless Hamiltonian. We determined the spectrum numerically for various volumes and quartic couplings and both for periodic and antiperiodic boundary conditions. In the periodic case we found it useful to adapt the mini superspace approach. We compared the results with semiclassical calculations and also with finite size corrections. These finite size corrections are expressed in terms of the infinite volume masses and scattering matrices, for which we proposed parametrizations from the unitarity relations. We summarized the infinite volume mass spectrum and scattering parameter as functions of the coupling on figure 24.
We found a very good agreement with the semiclassical results for the breather spectrum for moderate couplings. We then mapped the full spectrum of the neutral and kink excitations. We found at most two neutral excitations and each particle disappeared from the spectrum when its mass exceeded twice the mass of the lightest particle. This confirms the conjecture of [19] that only two neutral excitations can exist in the spectrum.
In determining the kink mass we used both the splitting of the even and odd groundstate energies with periodical boundary conditions and also the lowest lying state in the antiperiodic sector. These two methods agreed very well and the results were close to the semiclassical spectrum. For large couplings we identified the critical point beyond which the broken Z 2 symmetry is restored. We also confirmed that the second order phase transition is within the Ising universality class. In computing the leading exponential mass correction of the lightest neutral excitation we found relevant deviations from the Lüscher formula (5.30). By comparing the spectrum with the analogue breather spectrum of the sine-Gordon theory we got convinced that this effect must be attributed to non-integrable, inelastic processes, which need to be understood further.
We believe that improving our data with the renormalization group method of [14] enables to measure other interesting quantities such as decay rates or the inelastic part of the scattering matrices. Another interesting direction of research is to describe the strongly coupled spectrum by perturbing the conformal Ising model and to compare the analytical results to the improved numerics.

Acknowledgments
We thank Zoltán Laczkó for a collaboration at an early stage. ZB thanks Joao Penedones for the discussions on dispersion relations and the parametrization of the scattering matrix, while ML would like to thank László Ujfalusi for useful discussions regarding the used numerical eigensolver package. We also thank Slava Rychkov and Lorenzo Vitale for suggestions and for comments. This work was supported by a Lendület Grant and by OTKA K116505.

A Ground-state energy from normal ordering
In this appendix we motivate the ground-state energy

JHEP10(2016)050
of the free massive boson from a normal ordering prescription. The free normal ordered Hamiltonian in volume L takes the form Normal ordering is understood that we expand φ as with ω n = µ 2 + k 2 n and k n = 2πn L and then put all the annihilation operators a n on the right of the creation operators a + n using the rule [a n , a m ] = [a + n , a + m ] = 0 ; [a n , a + m ] = δ n,m (A. 4) This implies that where we introduced a cut-off to make them finite.
In the paper we are interested in the following perturbation of the system where µ 2 = µ 2 + δµ 2 . We would like to understand in which way it is related to H µ ,L . They formally look the same -however, they are normal ordered at different mass scales. As a result, they differ by a constant term. To calculate this constant, we use : φ 2 : µ,L =: φ 2 : µ ,L +∆Z ; : (∂φ) 2 : µ,L =: (∂φ) 2 : µ ,L +∆Y ; where ω n = µ 2 + k 2 n and (∂φ) 2 = (∂ t φ) 2 + (∂ x φ) 2 . Observe that all the sums are finite, although this would not be true separately for (∂ t φ) 2 or for (∂ x φ) 2 . In the first term we can subtract the infinite volume limit

JHEP10(2016)050
to obtain Similarly, so that we get This implies that We can reformulate this result as 16) That is, if we normal order the interaction term at infinite volume and introduce the groundstate energy by E 0 (µ, L) and ground-state energy density by E bulk (µ) = µ 2 8π − µ 2 4π log µ then the perturbation will change these parameters in a consistent way. In the following, we do not keep track of the energy density and introduce the perturbing operator simply as : φ 2 : µ,∞ .
One can also wonder how the ground state energy changes by changing the volume, i.e. to relate the system with mass µ and volume L to the same system with volume L . Clearly the spectrum of the two systems are completely different. What we can compare is the ground state energy. We will demand that LH depends only on µL. This is satisfied for This implies that By this rule we can calculate the ground state energy at volume L and mass µ. Summarising, the vacuum energy at volume L and renormalization scale µ is E 0 (µ, L), which is independent of the path we arrive at this point, i.e. by adding perturbations or by changing the volume. But we do not keep track of the change in the vacuum energy density.

B Parametrizing the scattering matrix
Here we analyze the 2 → 2 scattering matrices of the theory. We first recall how the unitarity of the scattering matrix can be used to gain information on its analytic structure. The idea is to write the unitarity equation SS † = I for the interaction part of the matrix Taking matrix element between initial and final states and inserting a complete system we can write For diagonal matrix elements, this equation implies the optical theorem, namely the absolute square of the total cross section is related to the imaginary part of the forward scattering amplitude. This equation shows that S-matrix elements are not independent, they are connected analytically by a very intrinsic way. Let us focus on the 2 → 2 particle scattering process and denote the momenta of the initial particles as p 1 , p 2 and outgoing particles as p 3 , p 4 , which are all onshell p 2 i = ω 2 i − k 2 i = m 2 1 . Energy conservation can be factored out where due to relativistic invariance the scattering matrix depends on the Mandelstam variables On-shell condition implies that s + t + u = 4m 2 1 . In two dimensions the kinematics is further restricted since t = 4m 2 − s. The analytical S-matrix theory states that S(s) can have poles and cuts dictated by unitarity, 1 but otherwise it is an analytical function of s, which also satisfies crossing symmetry. Spelling out the equation (B.2) for the two-particle process we obtain: In factoring out momentum conservation we assume that k 1 > k 2 and k 3 > k 4 such that k 1 = k 3 and k 2 = k 4 must hold. Let us write out the first few terms. Clearly the vacuum does not contribute. The one-and two-particle contributions read as 2 eS(s, t)= dk n δ (2) (p 1 + p 2 − p n )Γ n 12 Γ 34 n + (B.6) + dk n 1 dk n 2 δ (2) (p 1 +p 2 −p n 1 −p n 2 )S(p 1 , p 2 , p n 1 , p n 2 )S(p n 1 , p n 2 , p 3 , p 4 ) + . . .
Here we factored out the total energy momentum conservation. However, on the r.h.s. other δ functions appeared, which eat up one integration leading to an energy delta JHEP10(2016)050 function in the one-particle and a constant piece in the two-particle contribution. Since these contributions appear only for certain analytically continued values of s, they tell us the non-analytical positions of the scattering matrix on the complex s-plane. Bound states correspond to poles, where energy conservation can be maintained, and the two-particle threshold gives a cut starting from s = 4m 2 1 . Let us go into the center of mass frame k 1 ≡ k = −k 2 . Then we have s = 4m 2 1 + 4k 2 and t = −4k 2 , such that relation t = 4m 2 1 − s is satisfied. Crossing symmetry implies that Switching from the i prescription to the −i prescription is equivalent to time reversal, which changes the scattering matrix to its inverse The analytical structure of the S-matrix is displayed on the following figure: The bound state poles are at m 2 2 and at 4m 2 1 − m 2 2 . There are cuts on the real line starting at the multiparticle thresholds (m n + m k + . . . ) 2 . The physical value of the S-matrix is at S(s + i ) where s > 4m 2 1 . In the following we parametrize the two-particle scattering matrix based on this analytical structure. We first switch to rapidity parametrization. In the center of mass frame k 1 = −k 2 = m sinh θ 2 , such that the rapidity difference is θ. The rapidity is related to the s variable as This resolves the cut starting at 4m 2 1 and maps the first sheet of the complex s-plane to the strip 0 < m(θ) < π. The bound state poles are located on the imaginary axes at θ = iζ and at θ = i(π − ζ), since m 2 2 = 4m 2 1 cos 2 ζ 2 (B.10) The first branch cut depends on the masses. Assuming the kinks are very heavy they appear at s = (m 1 + m 2 ) 2 or in the rapidity parameter at θ t :

JHEP10(2016)050
If the kinks are lighter, then even earlier. Now we turn this information into a parametrization of the scattering matrix. Recall that crossing symmetry translates to the rapidity parameter as S(iπ − θ) = S(θ) (B.12) while the S-matrices on the two sides of the cut are related as In order to have the right periodicity and good asymptotic behaviour let us consider the logarithmic derivative of the scattering matrix: 14) It has the properties Let us write it as where the contour surrounds infinitesimally θ. Now we blow up the contour and deform it to θ = 0 above the cuts and to θ = iπ below the cuts. In doing so we pick up the residues of the bound state poles: (B.17) where ρ(θ) is the jump of φ(θ) on the cut, which is related to the inelastic scattering processes through the unitarity relation. After integrating and exponentiating we obtain where we made the integral to converge for large θ . We note that similar parametrization should be valid for the kink-kink and kink-breather scattering matrices.

C A review of the numerical algorithm used
In this appendix we review the numerical implementation of our method.

C.1 Overview
The numerical results were obtained using an application written in C++, that is supplied with this document. We tested the simulation under SuSe and Ubuntu Linux platforms. As explained earlier,the matrix elements of the Hamiltonian are calculated in a truncated, finite-dimensional basis and then the lowest eigenvalues of the obtained large sparse matrix is determined. With the use of the PRIMME iterative eigensolver package [28] (of JHEP10(2016)050 which we currently use version 1.1), it was possible to work with matrices of up to 250 million nonzero matrix elements on a single personal computer, the computing time being a couple of minutes per measurement point. The eigensolver uses an improved version of the Jacobi-Davidson method (referred as JDQMR ETOL in PRIMME).
The algorithm itself is generally capable of handling a Hamiltonian composed of the sum of a finite number of Landau-Ginzburg interactions in 1+1 dimensions put into a finite volume. In this article, we focused on the parity-conserving ϕ 4 model, and thus the parameter corrections are only implemented to this case. In the periodic sector, it proved to be sufficient to use the minimal Hilbert-space approach, with the separate treatment of the zero mode and the construction of the overall Hamiltonian as a direct product of the minimal and non-minimal subspaces. An emphasis was also put on the antiperiodic sector, though obviously the minimal space approach is not applicable in this case.
A single run of the program usually generates several output files, each one corresponding to a different energy cutoff. A single file contains the relevant eigenvalues for different volumes between a minimal volume L min and maximal volume L max , having step L step , the interaction strength parameters being corrected corresponding to the change of volume, while all other parameters are kept fixed.
The parameters of the simulations can be set in a separate configuration file, in which the user can set the number of eigenvalues to be calculated, along with the interaction strengths, truncation dimensions, and output file name. It is possible there to decide the usage of the minimal space method, choose the boundary condition and select the desired parity and momentum sector. A basis mass and volume must also be given for the construction of the free Fock basis. The examined volume interval and step size are also to be given in the aforementioned file.
The source code is structured around C++ classes and is divided into several files. The entry point of the program is in the file main.cpp. Some auxiliary routines and numerical functions are declared and defined in utils.h and utils.cpp. A numerical integration routine from Numerical Recipes [29] is implemented in numint.cpp. The classes HilbertState and HilbertSpace, which are used to build the Fock basis, are defined in phi4.cpp. The algorithms regarding the calculation of a Landau-Ginzburg interaction term between two Fock states are implemented in phi n.cpp. The class Hamiltonian is realized in minspace.h and minspace.cpp. A class to store a sparse matrix in coordinate list (COO) format is defined in coo.cpp, along with a number of routines particularly useful for our considerations. The header primme.h provides the interface to PRIMME functions. The connection between the simulation and the eigensolver is realised in the utility functions of PRIMME int.h. Most notably, as PRIMME uses a different sparse matrix format (Compressed Sparse Row, CSR), there is a routine to convert the COO format sparse matrix to one that is operable by the eigensolver.

C.2 The basis
A state of the Fock space is represented by the class HilbertState in the algorithm. In fact, the particle content can be written in conventional representation |u = |(k 1 , n 1 ) , (k 2 , n 2 ) , . . . , (k r , n r ) where this notation stands for having a state |u that contains a number n 1 of particles with momentum k 1 , and a number n 2 of particles with momentum k 2 , and so on, until all particles are counted for in the state. Due to the finite volume, the momenta are quantized so that in the periodic sector: while in the antiperiodic sector, Any Fock vector is thus identified by a finite series of integer pairs. A series like this can be stored in a variable of type HilbertState.
The HilbertSpace class contains an array of HilbertState variables. Crucially, this array is filled by a member function in a determined way, so that the HilbertSpace contains all Fock vectors below a given energy cutoff. Starting from an auxiliary vacuum state (a HilbertState variable that is empty), the construction is performed recursively, each recursive step having two parts: first we add an additional particle of some momentum i to the auxiliary Fock state containing only particles with positive momenta. Then we distribute the signs of particle momenta every possible way so that the overall momentum lies in the appropriate sector. The construction routine is called with a starting auxiliary state, and has a loop in which the function calls itself with another auxiliary state having one additional particle. If the energy cut is reached, the loop breaks and the routine returns. If the embedded routine returns, the momentum i of the added particle gets increased. At the end, the members of the array get sorted according to their energy (with respect to the basis volume and mass).
The truncation of basis is also implemented in the HilbertSpace class; it corresponds to cutting off the highest energy states from the array.

C.3 Calculating interaction matrix elements
Due to the remarkable efficiency of the eigensolver routine, the construction of the sparse matrix and evaluation of matrix elements can also turn out to be the bottleneck of computation speed. The algorithm was designed to fit two criteria: firstly it should be able to decide very quickly if a matrix element is zero or not, and secondly, it should calculate the actual values of matrix elements as effectively as possible.
A Landau-Ginzburg interaction term of degree N reads

JHEP10(2016)050
where k i runs through the one-particle momenta allowed by the boundary condition, = ±1, a +1 n ≡ a + n , and a −1 n ≡ a n . First it should be noted that the matrix element a| O N |b can have a nonzero value only if the state |b can be transformed into a state equivalent of |a by a product of at most N creation-annihilation operators. Therefore the problem of fast discrimination of zero matrix elements reduces to using the particle content of |a to generate all states |b with energy E |b ≤ E |a and a| O N |b = 0 and looking them up in the array of HilbertStates. Secondly, the actual content difference of |a and |b determines which operator product terms from the sum in (C.1) contribute to the matrix element. (Since the simulation imposes normal ordering at the volume of computation, we don't have to consider terms of "tadpole" type, where new particles not contained in any of the states would pop up and annihilate each other.) The determination of the content difference is facilitated by the class StateDifference. This structure contains two arrays: the first array lists the momentum of all particles in the bra vector that are not present in the ket vector; the second contains the particles only present in the ket vector. The sum of array sizes gives the total number of creationannihilation operators needed to transform the states into each other. This total number will be called particle difference for simplicity. The comparison of the two HilbertSpace variables and the construction of these arrays is done by a constructor of StateDifference.
If the particle difference of two HilbertStates is exactly N , then there is only one type of contributing normal ordered operator product from (C.1). However, if we only count operator product terms that are different after normal ordering, we always have to remember that the action of normal ordering effectively produces a multinomial coefficient for each term. If the particle difference is a value N pd < N , then we may also annihilate and create N − N pd particles taken arbitrarily from the ket vector. Since the matrix is symmetric, we can interchange the bra and ket vector so that the number of annihilation operators is always greater or equal than the number of creation operators in the term under consideration. Thus every contributing term can be accounted for by examining the ket vector particle content and the particle difference. The contributions to the matrix element of O n are then calculated term by term. This computation is done for a single O N term using the class phi n interaction.

C.4 The Hamiltonian
As mentioned earlier, the matrix representation of the truncated Hamiltonian in the generated Fock basis is stored in a sparse matrix. The matrix is technically stored in a variable of type sparse matrix. A single matrix element is kept in a struct variable such that the row and column index is stored explicitly along with the value of the element. Tha class sparse matrix consists of an array of these structs, along with some useful member functions.
The interactive part of the Hamiltonian consists of a sum of terms O N in the form (C.1). The parameter corrections according to the scheme change are accounted for in the main function of the program.
The Hamilton operator is technically realized by a class Hamiltonian in the program. If the minimal space approach is not used, then at the beginning of the calculation, the JHEP10(2016)050 matrix is filled with the matrix elements of the interaction, the non-diagonal ones being only saved if they are nonzero. In the Fock basis, the free Hamiltonian is diagonal, and its matrix elements are also added to the matrix. After this initial construction, only the nonzero matrix elements are modified.

C.4.1 Minimal space approach
If this approach is used, the program initially calculates the normalized (i.e. g N = 1) matrix elements of O N between non-minimal basis states for every positive N smaller or equal than the maximal interaction exponent appearing in the Hamilton operator, and store the representations in an array of sparse matrices contained in the class Hamiltonian.
For each examined volume, the minimal space spanned by zero modes are accounted for first. This basically corresponds to solving an anharmonic oscillator with appropriate coefficients and obtaining the first eigenvalues and eigenvectors in the free basis. The ϕ 0 terms of the interaction are accounted for in the minimal space calculation. As a second step, the matrix elements of interaction terms O N are calculated between the numerically determined "exact" eigenvectors in the minimal space, and get stored as an array of (dense) matrices. Then the non-minimal (normalized) matrix elements of every O N get actualised in the stored array of sparse matrices. In the last step, the representation of the truncated Hamilton-operator is obtained as a sum of direct product terms according to (4.6).
If the computation is limited to a single parity sector, then the zero-mode oscillator is calculated separately in the even and odd particle number basis. The minimal-space O N matrices in the exact basis are constructed so that the basis vectors containing an even number of particles are ahead of those having an odd number of particles. Then the nonminimal representations of O N s are also transformed so that the even basis elements are ahead of the odd basis elements. Then the overall Hamiltonian is built up in the desired sector using the appropriate submatrices.

C.5 A quick comparison of sparse-matrix eigensolvers: PRIMME and Arpack
We compared Arpack (used through the C++ interface Arpackpp) to PRIMME focusing mainly on a sample problem of basis size around N = 13000, g 4 = 0.1, with the 7 smallest eigenvalues computed, but we explored other parameters as well. When the minimal space was used, its dimension was chosen to be N 0 = 1000, and the 6 smallest eigenvalues were kept. In all cases, the accuracy of the residual norm was chosen to be around 10 −9 .
The exact computing times depend strongly on the structure and spectrum of the matrices as well as the number of computed eigenvalues. As a matter of fact, Arpack tends to have a minimum in computation time against the number of calculated eigenvalues at a finite value of the latter (which sometimes turns out to be larger than 7). The JDQMR-ETol method generally seems to be more appropriate when a smaller number of eigenvalues is needed. In addition, computing times for JDQMR-ETol seemed to scale slightly better with the dimension of the matrix than Arpack did, although this may also depend on the matrix itself.
In the periodic sector, without the separation of the zero mode, JDQMR-ETol (PRIMME) won, being usually about 2−3 times faster than Arpack. Separating the zero JHEP10(2016)050 mode, PRIMME was significantly faster in the solution of the banded problems in the minimal space (the computing time being at least 4 − 5 times less than Arpack's). However, after the assembly of the full Hamiltonian, the resulting eigenvalue-problem was solved by Arpack around 3 times more efficiently. In the antiperiodic sector, the two packages had similar speeds for 5 eigenvalues. For a larger number of eigenvalues, however, Arpack got more efficient.
From our tests, we concluded that as a zeroth order approximation, the efficiency of the two packages is similar. We believe that the main strength of our program is that the whole algorithm was implemented in C++, which resulted in a much faster calculation of matrix elements, that enabled us to set up our matrices about 10 − 20 times faster than would be possible with a Python algorithm. This determined the overall speed of the calculations.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.