LSZ in Action: Extracting Form Factors from Correlators Nonperturbatively in 2d $\phi^4$ Theory

In this paper, we compute multiparticle form factors of local operators in 2d $\phi^4$ theory using a recently proposed method [1] for efficiently implementing the LSZ prescription with Hamiltonian Truncation methods, and we adopt Lightcone Conformal Truncation (LCT) in particular for our calculations. We perform various checks of our results at weak and strong coupling, and elucidate the parametric behavior of truncation errors. This opens up the possibility to compute S-matrix in various strongly coupled models using the LSZ method in LCT.


Introduction and Summary
The S-matrix is probably the most important observable in all of particle physics, but computing it at strong coupling poses several challenges.Often these challenges arise from the fact that the S-matrix is defined in terms of asymptotic states which strictly speaking require particles to be prepared in infinite volume, whereas most nonperturbative methods rely on finite volume.Locality implies that it should be possible to approximately extract the S-matrix scattering amplitudes given a sufficiently large volume, turning the problem into a practical issue of how quickly a given finitevolume approximation approaches the infinite volume limit. 1 One nonperturbative methods that does at least formally allow infinite volume is Hamiltonian Truncation2 evaluated in Lightcone Quantization, which exactly preserves a continuous set of boosts.However, even here there is a sense in which finite volume effects sneak in.The problem is that in order to numerically diagonalize the Hamiltonian, one restricts the space of states to a finite-dimensional subspace of the full Hilbert space, and this subspace usually does not contain particles at arbitrarily large separation.In more technical terms, this problem shows up in the following way.Given the full set of time-ordered correlators of local operators of a theory, scattering amplitudes can be extracted using the LSZ prescription, by taking the residues of poles that arise when the external momenta are taken on-shell.However, poles arising from the presence of multiple particles in the 'in' state or 'out' state' typically do not arise at finite truncation, but instead are essentially smeared out by truncation effects.
Recently, [1] proposed a clever prescription for circumventing this problem.The idea was to extract the LSZ on-shell poles by taking matrix elements of derivatives of operators of the form A(x) ≡ −(∂ 2 + m 2 p )ϕ(x), (1.1)   where m p is the physical pole mass and ϕ is a local operator that creates the particle of interest, and then to use the equations of motion to express A in terms of operators without time derivatives.For instance, in ϕ 4 theory, i.e., the theory of a real scalar field ϕ with a potential V (ϕ) = m 2 0 ϕ 2 /2 + λϕ 4 /24, the equations of motion allow one to replace A with as long as this expression is used inside Wightman correlation functions.This method was tested in the 3d O(N ) model at large N in [1], as a method for computing scattering amplitudes.Our goal in this paper will be to apply the same method in 2d ϕ 4 theory, without invoking any large N limit.We will focus on a similar but comparably simpler observable, namely the two-particle form factor for the stress tensor T µν between the vacuum and a two-particle states: We will test our results by comparing them to perturbative calculations, specifically at one-loop and two-loops, as well as to strong coupling calculations from [6] of this quantity computed using a different method.We will compute all correlators using Lightcone Conformal Truncation (LCT), a lightcone quantization Hamiltonian Truncation method [7, 8].Much of the new technical work involved will be to efficiently compute the matrix elements of ϕ and ϕ 3 between basis states when there is non-zero momentum exchange flowing through these operators.We will present three separate methods for these quantities, which provides a strong consistency check on the calculations.We also make example code available for implementations of all of them.Probably the most important new qualitative feature that arises away from the large N limit is that both terms ϕ and ϕ 3 in (1.2) contribute to the LSZ calculation.This fact is less innocent than it might at first appear.The point is that individually, the (Fourier transforms of ) correlators of ϕ and ϕ 3 diverge on-shell, due to the presence of the LSZ pole ∼ 1/(p 2 − m 2 p ).Therefore, it is not a priori obvious that in a numerical calculation, their sum will be sufficiently well-behaved to give even a finite answer, much less the correct one as the large truncation limit is taken.Nevertheless, we will indeed see that in our perturbative calculations the divergent pieces in ϕ and ϕ 3 do correctly cancel.In fact, we will explicitly see that they individually grow like O(∆ max ) at large truncation, but with equal and opposite contributions to the correlator of A, and the residual finite pieces after this cancellation match the perturbative results from Feynman diagrams.The convergence of the results is fastest at lower energies, and we will see the best results for the Taylor series expansion coefficients of the form factor as a function of s in an expansion around s = 0. Larger values of s converge more slowly, and in fact if we take the s = ∞ limit first then they never converge.Turning to the strong coupling regime, a slightly more complicated picture emerges.The divergent O(∆ max ) pieces in ϕ and ϕ 3 still cancel, but at strong coupling there are O(1/∆ max ) truncation errors in the value for m 2 0 − m 2 p in (1.2) that contaminate the residual finite terms in the form factor.We find that fixing the value of m 2 0 − m 2 p by matching the form factor at a single value of s leads to an accurate result over a very broad range of s.
This paper is organized as follows.In section 2, we discuss the LSZ prescription in the context of Lightcone Conformal Truncation, and review the proposal from [1].In section 3, we describe three different methods for computing the new matrix elements we need in order to evaluate the form factor.In section 4, we compare our results in perturbation theory to the one-and two-loop results computed using Feynman diagrams.In section 5, we compare our results at strong coupling to results from [6].Finally, in section 6, we conclude with a discussion of potential future directions.

Application of LSZ in truncation
2.1 Applying LSZ to two-particle form factor Our main object of interest is the two-particle form factor for an operator O, which can be extracted by applying LSZ to the following correlator: (2.1) We will work in lightcone coordinates for the momenta, p ± = p 0 ∓p 1 √ 2 , p 2 = 2p + p − .There are two cases we need to consider, one where p 1− and p 2− have the same sign, and one where they have different sign; without loss of generality, by Hermitian conjugation we can take p 1− > 0. Because there are no physical states with p − < 0, we have Take the case where both p 1− and p 2− are positive.Because of the time-ordering and the observation above, contributions to the integral vanish unless y + < 0 and z + < 0. We can focus on the case 0 > y + > z + and obtain the case 0 > z + > y + by symmetrizing p 1 ↔ p 2 .Moreover, we will insert a complete set of states as follows: where µ i labels the (in our truncated system, discrete) set of eigenvalues of the operator P + , through P + |µ i , p⟩ = (2.4) Now, since ϕ(x) = e iP •x ϕ(0)e −iP •x , we can perform the integrals over y, z, p and p ′ to obtain .
(2.5) To create a two-particle asymptotic state, LSZ instructs use to read off the residue of the pole at p 2 1 = m 2 p and p 2 2 = m 2 p .The pole at p 2 2 = m 2 p is manifest because the single-particle state is an eigenstate of the Hamiltonian, i.e., it is just the eigenstate with µ 2 j = m 2 p .We define the usual wavefunction overlap factor Z as ⟨Ω|ϕ(0 Reading off the residue of this pole and dividing by Z 1/2 , we obtain where |p 2− ⟩ refers to the one-particle state with lightcone momentum p 2− . The function F + O depends on three variables, but in our lightcone truncation setup there is an exact invariance under boosts which we can use to boost to a momentum frame where (p 1 + p 2 ) − = 1.We will do so implicitly from now on unless explicitly stated otherwise, and to reduce clutter we will denote |µ i , (p 1 + p 2 ) − ⟩ as |µ i ⟩.That leaves only two variables, which we will generally parameterize in terms of p 2− and s defined as In terms of these variables, we have There is a manifest pole in this expression when s passes through the value of any physical state, as we expect.However, what is more difficult to see is where the pole at p 2 1 = m 2 p comes from.At fixed s, all dependence on p 2 1 enters through its implicit dependence on p 2− , which appears in the matrix elements ⟨µ i |ϕ(0)|p 2− ⟩ in the numerator.To make this dependence on p 2 1 explicit, we can solve the following relation for p 2− in terms of s and p 2 1 : The ambiguity in the choice of which solution one takes to this equation for p 2− goes away after adding the symmetric combination G + (p 1 , p 2 ) + G + (p 2 , p 1 ) from (2.4).
In the large truncation limit, the pole in p 2 1 has to emerge from the sum over the amplitudes ⟨µ i |ϕ(0)|p 2− ⟩.We will see that these functions are all individually polynomials in p 2− , so a true pole cannot emerge at any finite value of truncation. 3nstead, we will be forced to apply the LSZ prescription indirectly.
Before moving on to strategies to overcome this issue, let us briefly revisit the second case referred to above, where p 1− > 0 and p 2− < 0. In this case, due to the time-ordering, the only integration region in y + and z + that contribute are y + < 0 and z + > 0. Inserting a complete set of states as before, (2.11) The analysis proceeds similarly to before, except that now both the poles in p 1 and p 2 arise from one-particle eigenstates and are manifestly visible even at finite truncation.Taking p 2 1 = p 2 2 = m 2 p and reading off the residue of the pole from the one-particle eigenstates, one obtains This quantity can therefore be easily computed in the truncation framework, as was done in [6].Its main disadvantage is that it requires taking one in state and one out state, and so does not directly probe the interesting kinematic region with a two-particle asymptotic state.

Introducing Equation of Motion
In the previous subsection, we discussed why it is difficult in truncation to directly read off the residue of poles in order to create multi-particle asymptotic states.The strategy proposed in [1] is to evaluate the differential operator i(∂ 2 + m 2 p ) acting on ϕ before taking the momenta on-shell.That is, we compute It is crucial that the mass m p is the physical mass, not the bare mass.Note that we have already used the manifest pole for the first particle to produce |p 2 ⟩ in the ket state.When the derivatives ∂ 2 y hit the time-ordered correlator, there is a contribution from the time-derivative hitting the time-ordering θ functions, and another contribution where the derivatives all act within the time-ordering brackets.Consider the former contribution.Writing ∂ 2 = 2∂ + ∂ − and the time-ordering as we have 15) The first term on the RHS can typically be simplified by using the canonical commutation relations or, if O is a conserved current, Ward identities.We will mainly focus on the case where O = T −− is the energy-momentum tensor, in which case the commutator can be evaluated for instance as follows.In lightcone quantization, the canonical commutation relations are and therefore (2.17) since T −− = (∂ϕ) 2 .Substituting this into the expression for G + O , it generates the following contribution: where we have used ⟨ϕ(x)|p⟩ = e −ip•x .This term is the tree-level contribution, which we will separate out from the rest of the form factor and concentrate on computing the remaining term, from ⟨Ω|T {O(0)∂ 2 ϕ(y)}|p 2 ⟩.When p −1 = −p 2 − , (i.e., s = 0), the tree-level contribution does not receive any higher order corrections, due to the Ward identity for T µν .The remaining terms, where ∂ 2 + m 2 p act within the time-ordering, are more complicated.Inserting a complete set of states and repeating the analysis that led to (2.9), the only change is that the terms where i(∂ 2 +m 2 p ) acts inside the time-ordering become The strategy proposed in [1] is that we can evaluate ∂ 2 ϕ directly, using the equations of motion in λ 4! ϕ4 theory: Note that m 0 in this equation is the bare Lagrangian mass, not the physical mass m p .This strategy explicitly relies on using a Lagrangian theory, so that one knows the equations of motion for the fields.Before we discuss how to implement this approach, we will first take a step back and describe a more general strategy that only relies on the Hamiltonian itself, which can be obtained in non-Lagrangian theories. 4The trick of using the equations of motion above will then emerge naturally as a special case.
To begin, we write :

.21)
To evaluate the commutator, we can use the fact that we know how to compute the matrix elements of P + and ∂ϕ explicitly in a basis of our choosing, say with basis states |k⟩: When we sum over the intermediate states |k⟩, we have the option to insert a complete set of basis states |k⟩ for a much larger subspace than the truncated basis we use when we diagonalize the Hamiltonian.In fact, it is not too hard to see that it is necessary to take a larger basis for the intermediate states. 5n a Lagrangian theory, there is a slick way to evaluate this commutator with an infinite truncation for the intermediate basis states |k⟩ in (2.22): we can use the canonical commutation relations.This way of evaluating the commutator is equivalent to using the Euler-Lagrange equations of motion for ϕ.To see this explicitly, take the P + Hamiltonian for λ 4! ϕ 4 theory: and use the canonical commutation relations to evaluate [P + , ∂ϕ]: (2.26) Substituting into (2.21), where we have defined the mass shift as The advantage of this expression is that we have essentially cancelled the " " factor by using the equations of motion, so that we should be able to take both p 1 and p 2 momenta on-shell from the beginning.Imposing this constraint, and adding in the symmetric combination F − O , we finally obtain the following expression for the on-shell form factor: over eigenstates |µ i ⟩ of P + .Then the matrix elements of the commutator reduce to where we recall we are using the frame p 1− + p 2− = 1.In this case, (2.21) becomes

.24)
Note that when both particles are on-shell, then s = where p 2−,os and (1 − p 2−,os ) are the two on-shell values of p 2− , i.e., the solutions to ) .In order to evaluate (2.29), we need to construct matrix elements for ϕ and ϕ 3 between arbitrary basis states.The new ingredient compared to previous work is that the bra state and ket state can have different p − momentum, so that there is some momentum flowing through the operators ϕ and ϕ 3 themselves. 6Computing ϕ and ϕ 3 matrix elements In this section, we describe three separate methods for calculating the matrix elements of ϕ and ϕ 3 in a basis of states created by Fourier transforms of chiral primary operators acting on the vacuum: Previous work [8] developed efficient methods for calculating the matrix elements of the Hamiltonian P + in this basis, and the methods presented here will be generalizations that allow momentum to flow through ϕ and ϕ 3 .The first method works directly with the Fock space of the free massless scalar field.All basis states are expressed in terms of their wavefunctions in the multiparticle Fock space, and matrix elements of ϕ and ϕ 3 reduce to integrals over all the individual particle momenta.We call this the 'Jacobi Polynomials' method because the wavefunctions of states (3.1) created by primary operators are products of Jacobi polynomials.The second and third method instead first evaluate correlators of operators in position space, and obtain all the required matrix elements by performing Fourier transforms of these correlators; unlike the first method, only two integrals (one for the bra state and one for the ket state) have to be performed in this case.The second ('Wick Contraction') method evaluates the position space correlators by performing Wick contractions, and does not use (or rely on) the conformal invariance of the free massive scalar theory.By contrast, the third ('Radial Quantization') method significantly simplifies the evaluation of the position space correlators by using the fact that conformal transformations can be used to map the bra and ket primary operators to 0 and ∞, where they can be represented in terms of radial quantization modes. 7The computational efficiency of the Radial Quantization method is significantly higher than the other two methods and we were able to produce results up to ∆ max = 40.Having all three methods is a useful consistency check and we present numeric code for them at https://github.com/andrewliamfitz/LCT. Below, we will briefly describe these approaches, focusing on the novel aspects required for generalizing previous techniques to the case of the ϕ and ϕ 3 form factors.

Jacobi Polynomial Method
The simplest representation of the basis states we use in our truncation is in terms of the multiparticle Fock space states.Each single particle state is labeled by its p − momentum and has the following inner product norm: (3.3) Here, the number of particles n depends on the primary operator O.Because the operators O are defined to be primaries, each wavefunction F O is an eigenfunction of the conformal Casimir, which implies that the F O s are linear combinations of products of Jacobi polynomials of the following form [8]: where l = (l 1 , l 2 , ...l n ), P is the Jacobi polynomial, and |p| i ≡ i j=1 p j .In addition, because of Bose symmetry, the physical states are linear combinations of these polynomials that are symmetric in all the p i s.Schematically, the corresponding primary operators are was simplified by the fact that T −− is itself a primary operator.By contrast, ϕ 3 is not a primary operator, and computing its matrix elements is consequently more involved.
Given a representation of the form (3.3) for the bra and ket states, any matrix element of powers of ϕ can be computed by integrating over all the p i momenta.Computing matrix elements of ϕ and ϕ 3 between bra and ket states with differing momentum in this method is straightforward, and we provide example code that does this explicitly.
The main disadvantage is that the integrals become very computationally expensive very quickly as the number of particles in the external states grows.

Wick Contraction
Instead of using the wavefunctions directly at an early stage, another attempt is to keep the operator expression until the very end.The primary operators in our basis can all be expressed as a linear combination of "monomials" [8] where a n are numbers, and ∂ kn ϕ are monomials.As an example, for k n = (1, 2, 2), For any two-point function of primary operators, we can then expand it in the form of Similarly, when we want to evaluate the three-point function of ϕ 3 , i.e. the ϕ 3 matrix element in the orthoprimary basis, we can expand it as Where we can always set ϕ 3 to be at the origin due to translational invariance.With known coefficients for each composite monomial of a primary operator, the question is now to compute the three-point function of monomials.We can do this by performing Wick contraction, and rewrite the three-point function of monomials as a product of a three-point function and a two-point function, which we have closed form expressions for where k/k i means removing the element k i from vector k, and and recursively we can also write Here we are contracting the ϕ 3 in the middle with one ϕ on the left and two ϕ's on the right.Note that we require that |k/k i | = |k ′ /k ′ j,k | after Wick contraction, therefore k ′ always has one more element than k.This means that the number of particles in state created by O ′ is always one more than that by O.We call this the n → n + 1 piece of ϕ 3 matrix.Intuitively, this means that ϕ 3 annihilated one particle on the left, and two particles on the right, so that the remaining has a non-zero overlap.Similarly, we can construct the n → n−1, n → n+3, and n → n−3 pieces of ϕ 3 , as well as the n → n+1 and n → n − 1 pieces of ϕ.All of them have similar closed form expressions, therefore we can compute the exact three-point function of monomials in position space.
The next step is to Fourier transform this expression into momentum space, so that we can use the LSZ prescription and Equation of Motion as presented in the previous section.Notice that all three-point functions of ϕ and ϕ 3 we evaluate take the form of const.
Since the constant will not affect the result of Fourier transform, we can simply focus on which has an analytic solution Now we have an expression for three-point functions with monomials in momentum space, one simply needs to sum over these results with appropriate coefficients to obtain the three-point functions with primary operators, i.e. the matrix elements of ϕ and ϕ 3 .

Radial Quantization
In section 3.1, we presented the Jacobi method which computed integrals of Jacobi polynomials in a brute-force way.This method works in any QFT, without restriction to CFTs.In section 3.2, the Wick contraction method was also general for any QFT, but we used conformal symmetries to obtain a closed-form expression for the two-point and three-point functions.In this section, we present the radial quantization method, which is a way of quantizing states in a CFT, allowing us to exploit the orthogonality between states generated by monomials, giving much higher coputational efficiency.
In radial quantization, we expand ∂ϕ as where a k , a † k are the annihilation and creation operators in radial quantization mode.Then we can write the monomials in radial quantization as where Now we can easily see how this speeds up the computation.Note that the in state and out states are completely orthogonal, i.e. where Here BC k is the bin count of k, or the number of times that k show up in k.Then we can obtain the Zamolodchikov metric g OO ′ in a simple form Notice that in this case, we only have a single sum over k due to orthogonality, instead of a double sum over k and k ′ as in (3.9).
To compute the matrix element of ϕ and ϕ 3 , we need to go through the following procedure.As an example, we would like to compute the ϕ 3 matrix element The first step is to make a substitution ϕ 3 (y) → ∂ϕ(y 1 )∂ϕ(y 2 )∂ϕ(y 3 ) (3.26) We made this substitution because ∂ϕ is a primary operator and a monomial, while ϕ is not.With the radial quantization decomposition, it is easy to compute a correlator of monomials, as we can see in the later steps; the fact that ∂ϕ is primary also allows us to perform a conformal transformation on the entire correlation function, which facilitates the computation.We can simply retrieve the original three-point function of ϕ 3 by integrating over y i : Then we can decompose the right of (3.27) into monomials where We then need to map it to the radial quantization matrix element through conformal transformation.Define then we have where ∆ is the dimension of an operator.We use the symbol ∼ = because this equation is only true when the linear combination of k and k ′ add up to primary operators, since individual monomials does not have to obey conformal symmetry.Now we know that we can obtain G (∂ϕ∂ϕ∂ϕ) kk ′ (x, y 1 , y 2 , y 3 , z) by conformal transforming the radial quantization matrix elements G (∂ϕ∂ϕ∂ϕ) kk ′ (y 1 , y 2 , y 3 ), we can compute these matrix elements by decomposing the monomials into creation and annihilation operators Then we need to decide whether to pull out a creation or an annihilation operator from ∂ϕ.For ϕ 3 , there are 4 cases: n → n + 3, n → n + 1, n → n − 1, and n → n − 3.If we consider the n → n + 1 case for example, it means that there is one more a † in a † k ′ than a in a k .To make the entire function non-vanishing when acting on the vacuum, we need two a and one a † from (∂ϕ) 3 to balance the left and the right.Then we can write Then we can move a † l 3 to the left and a l 1 , a l 2 to the right, and contract with all possible a k and a † k ′ .Note that we can ignore the delta function brought by [a l , a † l ] as it is removed by normal ordering of ϕ 3 by definition.After these contractions we are left with Notice that the resulting expression is only non-zero when k/k i = k ′ /k ′ j , k ′ l , which means that k and k ′ can only defer by 2 elements.One can show that G (∂ϕ∂ϕ∂ϕ) kk ′ (y 1 , y 2 , y 3 ) is indeed a sum over terms of the form of If we want to transform the radial quantization matrix element back to the five-point function with finite position, we can perform the conformal transformation on each individual term Now we have the expression for the five-point function with (∂ϕ) 3 , and the next step will be to integrate over y i to obtain the matrix elements for ϕ 3 .The integration evaluates as when we take y = y 1 = y 2 = y 3 .Now if we sum over possible k i , k ′ j , and k ′ l , we can retrieve the ϕ 3 matrix elements of monomials and therefore reconstruct the matrix element in primary basis.
The next thing is to Fourier transform these elements to momentum space.For convenience, define then the Fourier transform of G To evaluate integrals of the form (3.39), first let us write them as Making the following change of variables the integration becomes gives a pole at w = −iϵ and w = 1 + iϵ.When we were integrating over z ′ , there was a term e iz ′ qw .The z ′ integral is performed by wrapping the branch cut from z ′ = 0 to z ′ = ∞, therefore when we evaluate e iz ′ qw , we always take z ′ positive.An important assumption we made is that p ′ > p, which means that q < 0. Therefore, to ensure that e iz ′ qw converges, we have to wrap the w integral in the lower half plane, where w goes to −i∞.This translates to picking up the pole at w = −iϵ when performing the w integral.To evaluate the integral, we simply need to expand (3.43) and extract the residue at w = 0 (in the limit where ϵ → 0). 8 In our code, we used one more step to speed up the computation.Notice that F (v) is indeed a product of g k (v), which can be expanded as a sum in the form of Now instead of expanding a product of powers of v, we only need to consider the Fourier transform of a single power of v each time, which reduces the processing time of residue computation.
In summary, we now have the following complete recipe to compute matrix elements of ϕ 3 : First, we evaluate G (∂ϕ∂ϕ∂ϕ) k i k j k l (y 1 , y 2 , y 3 ) for possible contractions k i , k ′ j , and k ′ l between two monomials, then transform it to the five-point function in position space G (∂ϕ∂ϕ∂ϕ) k i k j k l (x, y 1 , y 2 , y 3 , z).Then we integrate over y i to obtain the ϕ 3 elements from (∂ϕ) 3 , and Fourier transform the position space correlator to momentum space.Lastly, we sum over possible contractions to get ⟨∂ k ϕ(p)ϕ 3 (0)∂ k ′ ϕ(p ′ )⟩, and then sum over the monomials with appropriate coefficients to get the matrix elements of ϕ 3 in primary basis.

Cross-checking with Feynman diagram results
In this section, we will present the truncation results in radial quantization, compared with first order and second order corrections from Feynman diagrams assumed weak coupling.We will first compare the functions of form factor in the regime where s < 4. In this regime, the function is smooth, so that we can avoid the comparison between poles shifted by truncation effects, and observe a better agreement between truncation results and Feynman diagram results.By analytic continuation, if the results from truncation agrees with the analytic solution in the negative s regime, they should also agree in the positive s regime.For the one-loop result, we will also compare the integrated results in the positive s regime.From now on, we will take units of m 0 = 1 to avoid clutter.
By old-fashioned time-independent perturbation theory, we can expand the form factor to leading orders 2) , (4.1) where µ (0) is the unperturbed eigenstate of the Hamiltonian, and µ (1) is the first order correction, and the same notation applies to p 1 .The definition of ), we can extract the first order and second order corrections to the form factor, as we will present in more detail in the subsections.
When we want to compare the form factor with negative s, there are two ways to do it: the first way is to read off the residue of the pole from one-particle eigenstates, as presented in Equation (2.12), which we will call the t-channel method.Another way is to use the dispersion relation to reconstruct the entire form factor from just the imaginary piece in the positive s regime.Taking the imaginary part of form factor gives a delta function and in truncation, taking the integral on the right side of Equation (4.2) gives an expression for F (s) that can be evaluated for any complex value of s: where µ 2 i are the eigenvalues of the Hamiltonian.By contrast, our original expression (2.29), which took (p 1 + p 2 ) 2 = s, was restricted to the physical regime s > 4m 2 .When we take the imaginary part, the δ function in (4.3) enforces the on-shell condition (p 1 + p 2 ) 2 = µ 2 i and allows us to fix the value of p 1 .For this reason, when we use equation (4.4) to compute the form factor at negative s, we have to also choose p 1 to take different values for each term in the sum.

One-loop check
From Equation (4.1), we can extract the first order correction to the form factor 1 6 In Fig. 1, we show the first order corrections with and without applying LSZ in truncation.Notice that in truncation, we only need to consider two-particle intermediate states for µ (0) , and only the n → n − 1 piece of ϕ 3 interaction will contribute.The oneloop Feynman diagram gives the O(λ) contribution F 1 to the form factor (F ⊃ λF 1 ):   and we can compare this function with the truncation result in the negative s regime, as shown in Fig. 2. We can see that the two lines are almost in exact agreement, with a difference as low as 10 −10 .Next, we would also want to check the agreement in positive s regime with this method.In truncation, the eigenvalues of the Hamiltonian are always discrete, therefore it is impossible to reconstruct the poles perfectly, and the position of the poles are so sensitive to truncation effects that we cannot compare them with the Feynman diagram result.Rather than looking at a function with poles in positive s, we can compare the integrated imaginary part of the form factor: where µ n means that we require the intermediate state to be on-shell.In truncation, we now have a sum over step functions which are less sensitive to truncation effects.
Indeed, we can exploit the fact that we only need to consider two-particle inserted states, and generate an analytic result for the one-loop correction.The closed-form expression for the two-particle eigenvalues in the free theory are [8] Moreover, although we do not have an analytic derivation, it is easy to verify empirically that the following relation holds approximately at large ∆ max , up to higher order in 1/∆ max corrections: At this order, the matrix elements ⟨µ i |ϕ 3 (q)|p 1 ⟩ do not depend on q.Using the large ∆ max expression for the eigenvalues µ j , we can analytically do the sum over eigenstates to evaluate the form factor: where we have used Taking the imaginary part of this expression, we see that it matches the Feynman diagram result A comparison between analytic and truncation results in the positive s regime is shown in Fig. 3.With 484 intermediate states involved, we can see that the truncation is successfully reproducing the analytic result, which is another check that applying LSZ in Hamiltonian truncation is a valid method.

Two-loop check
From Equation (4.1), the second order corrections to the numerator are 3 2  and their corresponding diagrams are shown in Fig. 4.
In addition to the quantities in (4.12), we also need to consider the correction from the denominator: (4.13)The complete second order correction reads 1 6 (4.14)There are two new qualitative features about the form factor at second order as compared to first order.The first is that now there are four-particle states that must be included in the sums.This is not a significant conceptual difference, but it does make it more difficult to obtain analytic calculations of the form factor in truncation since the four-particle states are more complicated than the two-particle ones, and it also means the maximum ∆ max we can take in practice will be lower since the number of four-particle states grows more rapidly with ∆ max .
More significantly, at second order there is no a contribution both from the ϕ and the ϕ 3 term in the equation of motion.This fact introduces an interesting subtlety.Because we evaluate the form factors on shell, LSZ tells us that in the exact second order form factors, the contributions from ϕ and ϕ 3 individually will be diverge due to the 1/(p 2 − m 2 ) pole.Instead, it is only the sum of these two terms as they appear in the equation of motion that will remain finite on-shell.Since we are evaluating the form factor with a finite ∆ max , in our case the contributions from ϕ and ϕ 3 will remain finite on-shell; indeed, this is our main purpose for using the equations of motion.As ∆ max is taken to ∞, we should see the ϕ and ϕ 3 terms individual grow without bound.It is important to actually check in practice whether or not this divergence, combined with other possible truncation effects, correctly cancels out in the large ∆ max limit.In fact, it is not enough for the divergence to merely cancel -it is crucial that the subleading finite terms also produce the correct continuum limit result and are not spoiled by subtle residual truncation errors due to the leading divergence. 9lotting the imaginary part of the second order result at positive s is challenging because in truncation, the imaginary part is a sum over δ functions and the second order piece therefore involves derivatives of δ functions.Instead, we will focus on the real part of the second order result in the negative s regime.In Fig. 6, we present the comparison between truncation results after assembling the ϕ and ϕ 3 pieces as presented in Equation (4.14).The agreement is not as good as the one-loop result, but we can still see the convergence towards the Feynman diagram result as we increase the truncation dimension.
In Fig. 7, we present the quotient of the first four Taylor series coefficients of truncation and Feynman diagram results, and their quotient converges to 1 as 1/∆ max .In Fig. 8, we present the ϕ and ϕ 3 contributions to the form factor, and both of them diverge as ∆ max .This is the aforementioned result of the cancellation between ϕ and ϕ 3 pieces, as each term individually is growing while their sum is converging.
In Fig. 7, we can see that the higher order terms from truncation are closer to the analytic results.To explain this, when we Taylor expand the form factor, it can be written as Notice that the term constant in s has the lowest power of µ i in the denominator, therefore is the most sensitive to the truncation effect.In constrast, µ i in the denom-  inator has higher power in the higher order terms, therefore the truncation effects are suppressed.In the same truncation dimension, we always see that higher order Taylor coefficients are closer to the expected value.Finally, let us look in more detail at the contribution from the ϕ(q) to the two-loop form factor.In the large truncation limit,10 A simple consequence is that after combining the p 1 ↔ p 2 contributions, Following a similar derivation to the one leading to (4.10), the contribution to the form factor is so in this case, we can explicitly see that the contribution from ϕ(q) alone diverges like O(∆ max ) at large ∆ max .

Strong Coupling Result
Finally, in this section we will check our results at strong coupling.The main comparison we will make is with results at s < 0 obtained previously [6] by evaluating the form factor using the alternate expression (2.12) that only requires single-particle eigenstates, which are easy to obtain with Hamiltonian truncation.We compared these results at two different couplings λ = 6/π and λ = 36/π.As shown in Fig. 10, the truncation result matches the result from [6] to about 80% accuracy; however, we actually find that this accuracy does not seem to approach 100% at large ∆ max .What might be going wrong?As shown in Fig. 8, we expect the individual ϕ and ϕ 3 contributions to grow with O(∆ max ), and therefore a slight deviation of δm 2 of order O(1/∆ max ) can result in a major shift in the final form factor.We did not have to face this issue at weak coupling due to the fact that in that case the mass shift δm 2 converged sufficiently quickly, like ∼ ∆ −2 max , to its exact value that these small truncation deviations in δm 2 did not contaminate the residual finite terms ∼ O(∆ 0 max ) in the form factor.However, truncation errors in δm 2 become significant at finite coupling because now it converges more slowly.As shown in Fig. 9, at small coupling, δm 2 converges quadratically, whereas at strong coupling δm 2 converges linearly.
More conceptually, because ∆ max is a UV cutoff, it causes shifts in the bare parameters which can be absorbed into counterterms.At strong coupling, these corrections are O(1/∆ max ), and because the ϕ form factor alone has a O(∆ max ) divergence, the O(1/∆ max ) corrections to δm 2 lead to finite corrections in the final form factor.
There are various ways we might attempt to fix the value of δm 2 .In our previous computations, we determined δm 2 from the eigenstate spectrum of the truncated Hamiltonian.In particular, we obtained δm 2 by setting m 2 p = e 1 , where e 1 is the lowest eigenvalue of the Hamiltonian.Alternatively, we could try using the two-particle eigenstate to compute m 2 p by setting m 2 p = e 2 /4, where e 2 is the second lowest eigenvalue of the Hamiltonian, i.e. the two-particle mass.The results are shown in Fig. 11.In this plot, the truncation result appears to be converging to the t-channel result, but has a much larger error than the previous attempt and the convergence is slower than we would like.In fact the issue with slow convergence becomes worse in this case at weak coupling.The reason is that at weak coupling, the perturbative corrections to e 2 are smaller, so the truncation errors (which are present in e 2 even in the free theory) swamp them and therefore dominate the value of δm 2 until one reaches much larger truncations.In fact at the smaller coupling λ = 6/π, within the range of values of ∆ max that we used, the two-particle mass is greater than m 0 and δm 2 becomes negative, whereas the physical mass shift is positive.So at best, using the second-lowest eigenvalue from truncation to set the value for δm 2 appears to converge slowly.
Finally, instead of using truncation eigenvalues, we can try to fix the value of δm 2 by matching the form factor result with our LSZ method to the form factor result from [6] at a single value of s.We choose to do this at s = −5, and the resulting form factor is plotted in green in Fig. 10.The fact that the two results still match within 3% over the full range of s shown, after tuning only a single parameter, is a promising sign for this approach.In principle, the value of δm 2 should be ∆ max -dependent but not observable dependent, and it would be good to check if the same value of δm 2 gives accurate results when used across several different observables rather than just the T −− form factor.

Future Directions
In this paper we have focused on the two-particle form factor ⟨Ω|T µν |p 1 , p 2 ⟩ in 2d ϕ 4 theory, as probably the simplest possible example in which to try out the LSZ proposal from [1] and explore potential issues.The method itself is designed to be more general and hopefully it can be applied across a fairly wide range of models and multiparticle observables.The next simplest case to study would be a two-totwo scattering amplitude in 2d ϕ 4 theory, especially in the inelastic regime.important direction to explore is models where the light particle spectrum is composed of bound states, such as 2d QCD in particular, so that LSZ requires the use of composite operators which may complicate the strategy of using the equations of motion.Even more generally, there are models (one of the simplest being 2d Ising Field Theory with a σ deformation) without a Lagrangian description and therefore no obvious equations of motion at all.In this case, one might still resort to the strategy mentioned in section 2.2 of evaluating (∂ 2 + m 2 )O by taking the commutator [P + , O] with a much larger basis of states inserted as intermediate states in the matrix product.Although we used LCT in this paper, the proposal itself is more general and it should be possible to apply it using equal-time truncation methods applied to 2d ϕ 4 theory [10, 11].Finally, it is important to test the method in higher dimensions than d = 2.The original proposal [1] studied the O(N ) model in d = 3 in the large N limit, and a natural next step is to do 3d ϕ 4 theory, which has been studied in LCT in [12] and in d > 2 in [13-16].
Returning to the case we studied in this paper, there are still a number of issues to investigate further.At strong coupling, we had to choose the mass-squared shift δm 2 that appears in the LSZ prescription formula by matching the form factor to an independent calculation.It would be good to check if this same value of δm 2 agrees across multiple observables, and to have a deeper understanding of what fixes its value in the truncation framework.We also found that the accuracy of the form factor at a fixed truncation parameter ∆ max is better at small s and worse at larger s, and in fact if the limit s → ∞ is taken before the limit ∆ max → ∞, then the form factor apparently does not converge to the correct answer (see Fig. 6).It would be useful if there were a way to improve the method to obtain more accurate results at large s, or at least to have a better understanding of what range of s is expected to be valid for a given ∆ max .

µ 2 i
2p − |µ i , p − ⟩, and p − labels the (continuous) set of eigenvalues of the operator P − .So we have

Figure 2 .
Figure 2. Plot of one-loop result from Feynman diagram (black solid line) and truncation (red dashed line, at ∆ max = 20) in negative s regime.The residual plot shows the difference between the two lines.

Figure 3 .
Figure 3. Plot of s 0 ds ′ F 1 (s ′ ) comparing the analytic solution (black solid line) and truncation result from Jacobi method (red dashed line, included 484 intermediate states) in positive s regime.The residual plot shows the difference between the two lines.

Figure 4 .
Figure 4. Second order corrections in truncation.

Figure 6 .
Figure 6.Comparison between truncation and Feynman diagram 2-loop results, in negative s regime.The inset plot shows the same function at different scale, and one can see that for finite ∆ max , at sufficiently large s the truncation approximation eventually breaks down.

Figure 7 .Figure 8 .
Figure 7. Ratio between Taylor series coefficients of truncation result and Feynman diagram results, expanded at s = 0.The dashed lines are linear fits to Taylor series coefficients of different orders.

0 th order coefficients 1 st order coefficients 2 nd order coefficients 3 rd order coefficients
11Another