Convergence of the light-front coupled-cluster method in quenched scalar Yukawa theory

We explore the convergence of the light-front coupled-cluster (LFCC) method in the context of two-dimensional quenched scalar Yukawa theory. This theory is simple enough for higher-order LFCC calculations to be relatively straightforward. The quenching is to maintain stability; the spectrum of the full theory with pair creation and annihilation is unbounded from below. The basic interaction in the quenched theory is only emission and absorption of a neutral scalar by the complex scalar. The LFCC method builds the eigenstate with one complex scalar and a cloud of neutrals from a valence state that is just the complex scalar and the action of an exponentiated operator that creates neutrals. The lowest order LFCC operator creates one; we add the next order, a term that creates two. At this order there is a direct contribution to the wave function for two neutrals and one complex scalar and additional contributions to all higher Fock wave functions from the exponentiation. Results for the lowest order and this new second-order approximation are compared with those obtained with standard Fock-state expansions. The LFCC approach is found to allow representation of the eigenstate with far fewer functions than the number of wave functions required in a converged Fock-state expansion.


I. INTRODUCTION
The calculation of the bound states for a given quantum field theory is an inherently nonperturbative problem. Various methods can be applied, the best known being, of course, lattice (gauge) theory [1]. Here we consider a method based on a Hamiltonian formulation in light-front coordinates [2,3]. The fundamental bound-state problem is then the eigenvalue problem where P − is the light-front Hamiltonian, M is the mass of the eigenstate, and P = (P + , P ⊥ ) is the light-front momentum. 1 The Hamiltonian is constructed from the Lagrangian L for a generic field φ as − L : . (1. 2) The eigenstate |ψ(P ) has definite momentum P , and, once known, can be used to compute properties of the state This formulation is particularly convenient for the computation of form factors, because |ψ(P ) is boost invariant. The standard approach to the solution of the eigenvalue problem is to write the eigenstate as a Fock-state expansion, which leads to a coupled system of equations for the Fock wave functions. This coupled system is then converted into a matrix eigenvalue problem, either by direct discretization, as in discrete light-cone quantization (DLCQ) [4], or by basis function expansion, as in basis light-front quantization (BLFQ) [5]. However, a finite matrix representation requires a truncation of the Fock space.
This truncation has serious consequences. In particular, there can be uncanceled divergences, and self-energy corrections become dependent on the Fock sector and on the presence of spectator constituents. These are the nonperturbative analog of what would happen to the contribution from a Feynman diagram if the diagram were decomposed into the various time orderings, with the removal of the time orderings that involve too many intermediates. These difficulties led to the idea of sector-dependent renormalization [6][7][8], which has its own difficulties [9].
As an alternative, we have developed the light-front coupled-cluster (LFCC) method [10]. No Fock-space truncation is invoked. Instead, the eigenstate is written as coming from the action of an exponentiated operator T acting on a valence state |φ(P ) with √ Z a normalization factor. 2 The valence state is something simple that carries all the appropriate quantum numbers, in addition to the total momentum; for a proton in QCD it would be the three-quark state. The operator T increases particle number in various ways and conserves all the quantum numbers of the valence state; in QCD, T would include gluon emission from a quark or gluon and pair creation from a gluon.
The original eigenvalue problem is converted into two parts, through multiplication by e −T and projection onto the valence sector and its complement. To express this, we define 1 We define light-front coordinates [2] and momenta as x ± = t ± z, x ⊥ = (x, y), p ± = E ± p z , p ⊥ = (p x , p y ).
The mass-shell condition for the total momentum is then M 2 = P + P − − P 2 ⊥ . 2 This construction was inspired by the coupled-cluster method used in many-body problems of nuclear physics and quantum chemistry [11]. the effective Hamiltonian P − ≡ e −T P − e T and the projection P v onto the valence sector. We then have Roughly speaking, the first equation determines M and any wave functions in |φ , while the second determines the functions that define the structure of T . In reality, of course, they are a coupled system, unless the valence state has a single constituent and therefore no wave functions. All of this is obviously more complicated than the original eigenvalue problem, but it is exact. The power of the approach comes from the approximation step: Rather than truncate Fock space, we truncate T . Even for the simplest T operator, its exponentiation allows the eigenstate to span an infinite Fock space, and, without much difficulty, one can arrange the approximate e T |φ to fully explore all Fock sectors relevant for the quantum numbers of the valence state. In terms of a Fock-state expansion, what we have done is to force the wave functions of the higher Fock sectors to be directly dependent on those of the lower sectors rather than setting these higher wave functions to zero, as would happen in a Fock-space truncation. Yet another way to interpret the LFCC approximation is that the eigenstate is represented by a generalized coherent state. In any case, the avoidance of a Fock-space truncation eliminates the sector dependence and spectator dependence of self-energy corrections and potentially controls the uncanceled divergences.
The LFCC equations themselves are also truncated. The complement projection 1 − P v is restricted to the lowest set of Fock sectors necessary to have enough equations to solve for the functions that define T . This means that the LFCC method is not variational; the effective Hamiltonian P − is not Hermitian, and the truncated projections are not equivalent to minimization of the expectation value ψ|P − |ψ .
One price to be paid for the gains of the LFCC method is that the LFCC equations are nonlinear. The existence of a solution can be difficult to guarantee. However, a linearized perturbative solution shows that the LFCC equations re-sum perturbation theory to all orders for a restricted set of diagrams. (The restriction arises because of the truncation of T .) This implies that, for weak coupling, a physical solution must exist. Depending on the structure chosen for T and |φ , the physical solution may disappear as the coupling is increased. An explicit example of this appears in an application to φ 4 theory [12], where the solution for the lowest-order approximation for T does not extend beyond a certain coupling strength. This is likely due to the restriction of the valence state to a single constituent in a regime near the critical coupling where all Fock sectors should contribute strongly.
One question that immediately arises has to do with the convergence of the method, in the sense that as one relaxes the truncations of T and 1 − P v , how does the solution improve? The present work answers this question in a particular context, with an application to quenched scalar Yukawa theory [13] in two dimensions. 3 In general, the correspondence between perturbation theory and the LFCC resummation at weak coupling shows that the convergence of the LFCC method is closely related to the convergence of perturbation theory at weak coupling. To get beyond weak coupling, we compare a nonperturbative Fock-state expansion calculation to LFCC calculations done with T operators of increasing complexity.
The quenching of the theory eliminates potential concerns about the vacuum. Recent work [15][16][17][18][19] has emphasized the need for care in considering the vacuum on the light front, but here no vacuum bubbles can occur. This also means that the Fock wave functions of a massive state do not include vacuum contributions and therefore have a direct physical interpretation. This is not generally true in equal-time quantization, where one must compute the vacuum state as well as massive states; as an example, see the work on φ 4 theory by Rychkov and Vitale [20].
The Lagrangian, Hamiltonian, and Fock-state expansion for quenched scalar Yukawa theory are given in Sec. II. The formulation of the LFCC method for this theory is developed in Sec. III. The results for both the Fock-state expansion method and the LFCC method are presented and compared in Sec. IV, with a brief summary provided in Sec. V. Details of numerical methods and diagrammatic rules are left to appendices.

II. QUENCHED SCALAR YUKAWA THEORY
The Lagrangian for scalar Yukawa theory [13] is where χ is a complex scalar field with mass m and φ is a real scalar field with mass µ. The two fields are coupled by a Yukawa term with strength g. In two dimensions, the light-front Hamiltonian density is The mode expansions for the fields are 4 The nonzero commutation relations of the creation and annihilation operators are In terms of these operators, the quenched light-front Hamiltonian and Pair creation and annihilation terms are suppressed, to stabilize the spectrum. 4 Beginning here and for the remainder of the paper the + superscript of the light-front momentum is suppressed.
We seek eigenstates of P − , for which the two-dimensional light-front mass eigenvalue problem is (2.8) We limit this to the charge-one sector. In the next section, we consider the LFCC approach to the solution of this eigenvalue problem, but here we develop the standard Fock-state expansion approach, to use as a basis for comparison. We write the Fock-state expansion of the eigenstate as (2.9) Projection of the eigenvalue problem onto n ′ j a † (y j P )c † + ((1 − n ′ j y i )P )|0 , and division by µ 2 , yields coupled equations for the Fock-state wave functions ψ n Herem ≡ m/µ is a dimensionless relative mass and λ ≡ g/( √ 4πµ 2 ) is a dimensionless coupling strength. We solve this system numerically by first truncating the Fock space at n = n max neutrals and expanding the wave functions in a symmetrized monomial basis. The details are discussed in Appendix A, and the results in Sec. IV.
The structure of the eigenstate is studied by considering the relative probabilities for Fock sectors with different numbers of neutrals. These are formed as the ratio Results for these ratios are shown in Sec. IV.

III. LIGHT-FRONT COUPLED-CLUSTER METHOD
The LFCC method constructs the charge-one eigenstate in the form where c † + |0 is the single-particle valence state. The T operator is expanded in a sequence T = n T n , with The factor p n/2 is included to keep T n dimensionless; p is the natural scale, being the momentum flowing through the operator.
The action of T n is to increase the number of neutrals by n, and the exponentiation of T provides for generation of all possible (quenched) charge-one Fock states, even if T is truncated to only T 1 . Without truncation, the functions t n provide for an exact solution, with a duality between the t n and the Fock-state wave functions ψ n . However, without truncation the eigenvalue problem is equivalent to an infinite coupled system of nonlinear equations for these t n .
We can then study the convergence of the LFCC method as the number of terms in T is increased. Here we consider the first two terms, T 1 and T 2 , and compare results with those from the truncated Fock-state expansion.
The LFCC form of the eigenvalue problem is Independent of the level of truncation for T , the first equation becomes The contributions to this equation come from the P − 0 and P − int T 1 terms in P − , as represented diagrammatically in Fig. 1. On division by µ 2 , the projected valence equation reduces to 1 1 1 the following expression for the eigenmass M: The self-energy term is then specified by The function t 1 is to be obtained by solving the remaining LFCC equations.
With each truncation of T there is a matching truncation of the projector 1 − P v to include only enough Fock sectors to determine the unknown functions in the retained terms of T . Given the truncation to T = T 1 + T 2 , the equations for t 1 and t 2 take the form of two projections, onto the sectors with one and two neutrals. The contributions to the first projection come from Diagrammatic representation of the projection onto the one-neutral Fock sector.
These terms are represented in Fig. 2 and yield the following equation for t 1 : The equation (3.8) for t 1 can be obtained either by explicitly carrying out the contractions of annihilation and creation operators or by diagrammatic rules listed in Appendix B. The first term in the second line of Eq. (3.8) can be simplified by rescaling the integration variable x by 1 − y; this shows the integral to be equal to ∆. The same self-energy integral appears in the last term of the first line. The terms proportional to ∆ can then be collected with them terms, to introduce M 2 with use of (3.5) For the truncation T = T 1 , this equation, with the t 2 term removed, is all that need be solved.
The appearance of the physical mass M in the invariant-mass terms on the left of (3.9) is typical of the LFCC method, where self-energy corrections are independent of the Fock sector and independent of spectators. This avoids the use of the sector-dependent bare masses that are frequently introduced in truncated Fock-state-expansion calculations [6][7][8][9], where self-energy corrections are sector and spectator dependent.
The contributions to the second projection, onto the two-neutral sector, come from the following terms in P − : We solve these equations numerically, as discussed in Appendix A, both for t 1 alone and the coupled system, for t 1 and t 2 .    6. Same as Fig. 3 but for the P − terms 1 The relative probabilities for different Fock sectors can be computed from the expansion of the exponential form of the LFCC approximation The Fock state wave functions can be extracted by comparison with the Fock state expansion in (2.9), after the actions of the operators T 1 and T 2 are taken into account. We find The relative probabilities for the one and two-neutral sectors can then be computed as before, using (2.11). The necessary integrals can be done analytically for the basis function expansions introduced in Appendix A; however, for the cross term between the second and third terms of ψ 2 , the analytic result is the value of a hypergeometric function and that term is instead integrated numerically with Gauss-Legendre quadrature. The overall normalization Z is not computable in a finite sum, which is the motivation for considering relative probabilities, rather than absolutes. Fock sectors higher than the two-neutral sector can be considered, but the wave functions become much more complicated.

IV. RESULTS
The results for the mass M in the Fock-state expansion method are shown in Figs. 8-10. Both the basis size and the Fock-space limit are increased to achieve convergence for the lowest eigenstate; however, for the ultrarelativistic case ofm = m/µ = 10, convergence of the Fock-space expansion is not achieved for stronger coupling values, as can be seen in Fig. 10. On the other hand, convergence for the nonrelativistic case ofm = 0.1 is almost immediate.
From the solutions to the LFCC equations, we compute the mass eigenvalues M and the relative probabilities of the one and two-neutral Fock sectors. The masses are shown in Figs. 8-10, where we plot results for both T 1 alone and T 1 + T 2 .
Results for relative probabilities are plotted in Figs. 11-13. These show that as the neutral constituents become lighter, makingm larger, the importance of the higher Fock sectors increases dramatically. The LFCC approximation for the one-neutral Fock wave function yields a nearly exact match to the one-neutral relative probability; this is seen in Figs. 11-13, where the solid line representing the LFCC result passes through the points from the converged Fock-state-expansion results for the one-neutral probabilities. We interpret this agreement to mean that the effect of the higher Fock sectors on the one-neutral wave function is well represented by the LFCC approximation to these higher sectors.
The results show that the LFCC truncation to T 1 + T 2 is sufficient to replicate the converged Fock-state expansion results, with T 1 alone just as good as a two or three-neutral Fock-sector truncation. Thus the LFCC approximation, using only the two functions t 1 (y) and t 2 (y 1 , y 2 ) of one and two variables, respectively, is sufficient to represent information that the Fock-space expansion encodes in many more wave functions. In addition, the number of basis functions required to represent the Fock wave functions is significantly greater than the number required for the LFCC functions. Thus, the matrix representation is much smaller for the LFCC approximation, which is ample compensation for its nonlinearity. The failure of the nonlinear solver to converge 5 for strong coupling in the ultrarelativistic case occurs in the same coupling range where the Fock-state expansion fails to converge. This is near where M tends to zero and may be indicative of the incompleteness of theory. Quenching may have stabilized the spectrum, but the theory is no longer a complete quantum theory. As discussed in the Introduction, a similar lack of solution convergence has been observed in φ 4 theory [12].

V. SUMMARY
We have shown that the LFCC approximation provides an efficient representation of a massive eigenstate in quenched scalar Yukawa theory. We have also found that the LFCC approximation converges quickly as more terms are added to the T operator. From a numerical standpoint, there is also an efficiency in the basis size required for a matrix representation of the fundamental equations; the LFCC functions are fewer in number than the Fock wave functions, depend on fewer variables, and need fewer basis functions for their accurate representation.
In doing these calculations, we have developed diagrammatic methods for the construction of the LFCC equations. These significantly reduce the effort involved, compared to literally 0 Ɖ ƈ FIG. 9. Same as Fig. 8 but for a constituent mass ratio ofm = 1 and with both LFCC approximations T = T 1 and T = T 1 + T 2 . The basis sets in each Fock sector were limited to orders N = 2, 6, 12, 10, and 8 for n max = 1, 2, 3, 4, and 5, respectively. The basis sets for the LFCC results have maximum orders of N 1 = 5 and N 2 = 5.
carrying out contractions of creation and annihilation operators in matrix elements of the effective LFCC Hamiltonian. Extension to other theories should be straightforward.

ACKNOWLEDGMENTS
This work was supported in part by the Minnesota Supercomputing Institute through grants of computing time and benefited from participation in the workshop on Hamiltonian methods in strongly coupled quantum field theory supported by the Simons Collaboration on the Nonperturbative Bootstrap. The operator diagrams were drawn with JaxoDraw [21].
Appendix A: Numerical methods

Fock-state expansion
We solve the coupled system (2.10) for the Fock-state wave functions ψ n in (2.9) by first expanding the wave functions in a simple polynomial basis ψ n (y 1 , . . . , y n ) = y 1 · · · y n (1 − where m is the order of the polynomial P (n) mj , j is an index that differentiates distinct polynomials of the same order (which is nontrivial for multivariate polynomials), N is the maximum  Fig. 9 but for the mass ratio of the constituents ism = 10. The basis sets in each Fock sector were limited to orders N = 5, 5, 4, 3, and 3 for n max = 1, 2, 3, 4, and 5, respectively, and to N = 2 for all higher Fock sectors. The basis sets for the LFCC results have maximum orders of N 1 = 5 and N 2 = 3. In this case, the Fock-space expansion has not converged near M = 0. Also, the nonlinear system solver failed to converge for the LFCC approximation with λ beyond 10.6 when T 2 was included.
order included, and the c (n) mi are unknown coefficients to be obtained. The polynomials are chosen to be simple monomials, suitably symmetrized but not orthogonal. They take the form P (n) mj (y 1 , . . . , y n ) = y j 1 1 y j 2 2 . . . y jn n + · · · , (A2) with n i j i = m. The truncation of the basis to the order N is, of course, an approximation necessary for a finite matrix representation; we study convergence with the respect to this truncation, allowing N to be different for each Fock sector.
Projection of the nth equation onto each basis function, y 1 · · · y n (1 − i y i )P (n) m ′ j ′ (y 1 , . . . , y n ), yields a matrix representation of the original coupled system The individual matrices are T (n) m ′ j ′ ,mj = n i dy i m 2 y 1 · · · y n + ny 2 · · · y n P (n) m ′ j ′ (y 1 , . . . , y n )P The integrals can be done analytically in terms of the generalized β function dx 1 · · · dx n x k 1 1 · · · x kn n (1 − x 1 − · · · − x n ) = k 1 ! · · · k n ! (k 1 + · · · + k n + n + 2)! (A8) This allows for efficient calculation of all the integrals, with the different β-function evaluations done recursively and stored for use. If the basis functions were orthogonal, S (n) would be diagonal, of course. However, we implicitly orthogonalize the basis by performing a singular-value decomposition S (n) = U (n) W (n) U (n)T . The columns of the matrix U (n) are the eigenvectors of S (n) , and W (n) is a diagonal matrix of the eigenvalues. The U matrices then define an orthogonal transformation to new vectors of coefficients c (n)′ = (W (n) ) 1/2 U (n)T c (n) and new matrices, such as T (n)′ = (W (n) ) −1/2 U (n)T T (n) U (n) (W (n) ) −1/2 . The new matrix problem is no longer of the generalized type, but simply The lowest eigenvalue is extracted by standard procedures for symmetric matrices.
The convergence of such a calculation, with respect to the basis size, is illustrated in Fig. 14. Convergence is quite rapid in general; for stronger coupling values, near where M becomes zero, larger basis sizes are needed.
For the D matrix, a change of variables has been shown explicitly; similar rescalings are done for many of the other matrices. These rescalings arrange for the arguments of the polynomials to be polynomials and for all integration ranges to be from 0 to 1. The integrals are then linear combinations of simple integrals of monomials. The nonlinear matrix equations obtained in this way are then solved by a modification of the Powell hybrid method [22] as implemented in the general nonlinear equation solver 'hybrj' of the MINPACK set of subroutines [23]. The method is recursive; the initial guess for the unknown coefficients is taken to be zero for the lowest coupling strength and, as an increasing series of coupling strengths is considered, the next initial guess is the solution for the previous coupling strength.
For the case where only T 1 is included and we solve only (3.9) for t 1 with t 2 = 0, convergence with respect to basis size is very rapid whenm = 1. The results for N 1 = 1 and N 1 = 2 are indistinguishable on a graph. For the full solution, with both T 1 and T 2 included, the dependence on N 2 , the maximum order for the t 2 basis, is shown in Fig. 15. Convergence is again quite rapid, except for stronger coupling where M 2 approaches zero. For smaller and larger values ofm, convergence is slower for t 1 , requiring N 1 = 9 form = 0.1 and N 1 = 5 form = 10. Convergence for t 2 is quicker, using N 2 = 3, except for strong coupling in the case ofm = 10 where the nonlinear equation solver was unable to converge to a solution.  x − x 1 x 1 x − x 1 − x 2 x 1 FIG. 17. Diagrammatic representation of the T 1 and T 2 operators, and their corresponding expressions, acting to the right and creating one or two neutrals, respectively, by first annihilating a charged scalar with momentum fraction x.
the one-neutral projection, shown in Fig. 2. Except for the 1 2 P − int T 2 1 term in (3.7), there is only one diagram for each term in P − ; for 1 2 P − int T 2 1 there are two. The rules then yield (3.8).