{'fraction_non_alphanumeric': 0.032962759630341294, 'fraction_numerical': 0.009698966053618812, 'mean_word_length': 4.477634381656434, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Small x divergences in the Similarity RG approach to LF QCD

We study small x divergences in boost invariant similarity renormalization group approach to light-front QCD in a heavy quark-antiquark state. With the boost invariance maintained, the infrared divergences do not cancel out in the physical states, contrary to previous studies where boost invariance was violated by a choice of a renormalization scale. This may be an indication that the zero mode, or nontrivial light-cone vacuum structure, might be important for recovering full Lorentz invariance.


I. INTRODUCTION
Light-cone or light-front [1] description of field theories has long attracted attention for various reasons, often related to its peculiar kinematics [2]. Light-cone coordinates are used in vastly different applications usually in Hamiltonian formulation (see [3] for an extensive list of references), ranging from deep inelastic scattering and phenomenology to attempts to formulate the theory of everything. Our own interest in light-front field theories is motivated by an assumption that in light-front coordinates it is possible to derive from first principles a self-consistent and systematic constituent approximation to QCD [3]. We will give more details below.
As a consequence of the light-front kinematics, one component of the three-momentum, so called longitudinal momentum k + , is positively definite and can be interpreted as playing role of a Newtonian mass [4,5]. Thorn has long advocated that for this reason, light-front formulation of string theory is one of the best hopes for a truly fundamental, non-perturbative description of strings [6] and to quantitatively realize the conjectured correspondence between string and field theory [7].
Because of the above mentioned properties of the longitudinal momentum k + , light-front vacuum was thought to be trivial. This notion turned out to be naive, but the light-front vacuum indeed can contain only particles with k + = 0 known as the zero mode. Therefore, the light-front vacuum can be made trivial by imposing appropriate cutoffs. This is the last point that the light-front community agrees on. Opinions differ when it comes to the question how to proceed further.
Basically, there are two approaches to the problem of light-front vacuum [8]. One is to address it head-on [9]. Typically, this is done for various, often lower dimensional, field theories (for review see [10], for extensive list of references see also [3]) in the context of discreet light cone quantization (DLCQ) (for a recent review of the method as well as an expansive list of references see [11]). The theory is formulated in a box of a finite longitudinal size which makes k + discreet and non-zero. The zero mode [12] is then determined as a solution to constraint equations. Constraint and dynamical zero modes are often discussed. Critics complain that so far, it is not clear how much of the zero mode determined in this manner is an artifact of the method (in other words, a counterterm specific to the DLCQ) and whether the continuum limit exists 1 .
As an untypical example of direct studies of the zero mode, we would like to mention a recent work by Tomaras, Tsamis and Woodard [14] on back reaction in light-cone QED. Though motivated by the back-reaction in quantum gravity occurring on an inflating background, their work addresses some issues of the light-front vacuum without having to evoke DLCQ. They had constructed a full operator solution to a free QED coupled to a constant external electric field in continuum (3+1) dimensions. In this set up, all modes are forced to go through the zero mode at which point particle pairs are created. The zero mode of the constraint components of the fermionic field is shown to be crucial for unitarity.
The other approach to the problem of zero mode is more pragmatic. Instead of trying to solve for the zero mode, it is simply cutoff, be it with DLCQ [15] or an explicit infrared cutoff in a continuous formulation [3]. Physics associated with this mode can then be put in form of counterterms, if needed, for example, to restore symmetries or account for phenomena associated with the vacuum. Traditionally, spontaneous symmetry breaking was viewed as an example of such a phenomenon [3,12,13]. However, Rozowsky and Thorn [15] have argued recently that, while conceding that the inclusion of a fundamental zero mode is a valid theoretical option, it is not necessary to describe spontaneous symmetry breaking where its presence seems to be most needed. Indeed, in scalar quantum field theory in (1+1) dimensions DLCQ the physics of spontaneous symmetry breaking is completely and accurately described without the zero modes [15].
The need to put in new counterterms associated with the infrared (IR) regulator in our continuum formulation [3] was anticipated [13] but was not encountered yet in the applications to hadronic physics so far [16]. Perry [17] has shown that even though the one body and two-body effective operators are each separately divergent as k + goes to zero, the divergences exactly cancel in any color singlet state. The cancelation does not occur for non-singlet states, leaving them with an infinite mass. This together with a naturally generated confining potential (imprecisely referred to as "logarithmic") is a plausible feature of the approach. Note that both the effective confining potential and the infinite mass of the color non-singlets originate from small k + regions.
A shortcoming of the above described result is that it was obtained in a similarity renormalization group limited to matrix elements that required a specific infrared regulator (theta function) and an introduction of an arbitrary scale P + which violates an explicit kinematic symmetry of light cone. In a bound state calculation one can argue that there is a preferred scale, i.e. one associated with the typical longitudinal momentum of the state, or, the total center of mass P + ; however, consequences of the violation of the kinematic symmetry are not known.
Since then, the similarity renormalization group approach has advanced so that it is no longer necessary to violate the kinematic boost invariance, and to generate counterterms dependent on the total center of mass P + . It is also possible to use an arbitrary form of the small x regulator [18] (x being a dimensionless, boost invariant fraction of the total P + ). Thus, for the first time we are able to study with some degree of generality the issues related to the light-cone zero mode. Recently, G lazek has found [19] that even though infrared divergent terms cancel out in the running coupling, there is a residual finite dependence on the functional form of the infrared regulator. In this paper, we wish to study the issue of small x (or infrared) divergences in color singlet states consisting of a heavy quark and antiquark of the same flavor, for simplicity.
The paper is organized as follows. Section II is devoted to an overview of the similarity renormalization group approach light-front QCD. We start with a general description of the approach. For the sake of making this paper self-contained, the reader is reminded of the structure of the classical light-front Hamiltonian. The corresponding quantum Hamiltonian cannot be defined without regularization. Details of regularization determine structure and form of counterterms. In this case, details of regularization are particularly important, because only one of the regulators can be properly removed by renormalization. After introducing regularization, we present a technical description of the boost invariant similarity renormalization group for particles. Remaining sections deal specifically with quarkantiquark color singlets to second order in g. Section III contains the effective one body operators to second order in g; section IV the effective two body operators to the same order. In section V we address the issue of infrared divergences in those operators. We conclude with a short summary and conclusions.

II. SIMILARITY RENORMALIZATION GROUP APPROACH TO QCD
In this section we briefly review the similarity renormalization group (RG) approach to light-front (LF) QCD, introduced in ref. [3]. We show the unregulated canonical lightfront Hamiltonian. Then the regularization that we use is briefly explained, and finally, the similarity renormalization procedure is outlined.
The basic assumption behind the approach is that it is possible to derive a constituent picture of hadrons from QCD. To separate vacuum fluctuations, it is convenient to use light-cone coordinates with cutoffs preventing zero longitudinal momentum. Then vacuum is forced to be trivial. In addition, since such a cutoff introduces a nonzero minimum for the (kinematic, positive) longitudinal momentum, it also restricts number of particles in any state with a fixed longitudinal momentum P + . For massive particles, due to light-cone free energy increasing with decreasing longitudinal momentum, many-body states tend to have higher free energy than few body states. Mixing of low energy few particle states with many-body states can be expected to be naturally small, at least at weak coupling. These features make for useful prerequisites toward the constituent picture of hadrons.
The apparent difficulty with renormalization of light-front Hamiltonians (compared to Lagrangians) is turned into an advantage by using similarity renormalization [20]. The basic idea of the similarity renormalization group is simple. The regulated bare Hamiltonian that mixes all energy scales is transformed via a unitary transformation to a Hamiltonian that contains direct couplings only between neighboring scales. Such a Hamiltonian is referred to as "band-diagonal" in the sense that it does not allow for direct couplings between arbitrary scales. At any finite order of perturbation theory, a band-diagonal Hamiltonian cannot produce any ultraviolet divergences providing its matrix elements are finite. Therefore, by requiring that the band-diagonal Hamiltonian be independent of the regulator, counterterms that need to be added to the bare Hamiltonian can be identified.
Upon completion of renormalization, the effective band diagonal Hamiltonian is subject to diagonalization. We refer the reader to [16] for more details regarding this step of the procedure.

A. Classical light-front QCD Hamiltonian
The starting point is the canonical light-front QCD Hamiltonian in light-cone gauge, A + a = 0. We will not explicitly show terms that are not important for the specific calculations presented in the next sections. For a detailed discussion of the light-front Hamiltonian see [3,21]. Purely gluonic terms are studied in [19].
The part of the classical unregulated canonical Hamiltonian that is relevant for our study of the bound state of a quark and antiquark of the same flavor is, where H free is the free light-front Hamiltonian: and contains the standard order g quark-gluon coupling. Here ψ and A µ ≡ a A µ a T a are free light-front fields: and The constrained fields contain, in addition to ψ − and A − shown above, O(g) terms that are fully determined functions of the physical degrees of freedom. These produce new terms in the canonical Hamiltonian, among which is the so-called instantaneous gluon exchange between two fermions. As the next step, fields are expanded in mode functions consisting of (light-cone) spinors u pσ , v pσ [22], (p denotes momentum, σ denotes spin), SU(3) color spinors χ c and plane waves, The Hamiltonian is expressed in terms of what will become, after quantization, creation and annihilation operators of fermions (b † , d † , b, d) and gluons a † , a. For example, the kinetic part H f ree becomes and E k ≡ k ⊥2 k + are free light-front energies for massive and massless particle, respectively.
As another example, consider the interaction V 2 that contains (among other terms), where x i ≡ p + i /P + is the fraction of the total P + carried by the particle i. Momentum conservation is assumed implicitly. Similarly, V 1 is also rewritten.

B. Regularization
As will be illustrated in Section III, the ultraviolet (UV) and IR divergences in the lightcone QCD are interwined (for additional discussion see Apendix E of ref. [19]). For this reason and because at present only UV renormalization is achieved in a systematic manner, details of regularization are important. We regulate the canonical quantum Hamiltonian with cutoffs on the changes in transverse momenta and on longitudinal momentum fractions as in ref. [19], term by term in expansion of creation/annihilation operators. For example, an operator consisting of one creation and two annihilation operators would contain product of two regulators, one for each of the annihilation operators. Regulation of operators containing more than one creation operator (and their Hermitian conjugates) is slightly more involved. For details see [19]. The instantaneous interaction (9) is viewed for regulation purposes as an exchange of a virtual particle between two vertices. These general rules can be satisfied by various cutoffs. Specifically, we use a regulating factor of the form (arguments are generic for the moment) that is particularly convenient for analytical calculations.
Here and in what follows r δ is an infrared regulator that prevents its argument to become zero. We wish to study dependence on this regulator, hence we leave its functional form unspecified. ∆ is an ultraviolet cutoff. The limit with the cutoffs removed is achieved by ∆ → ∞, δ → 0, in which case r ∆δ → 1.
In order to preserve kinematic symmetries, we choose the arguments of the regulators to be kinematically invariant quantities, in particular, Jacobi momenta defined with respect to the vertex in question. As an example of a regulated interaction, the standard quark-gluon coupling has the following form in the regulated quantum theory, For clarity, summation indices are implicit here, sum is over spins and colors of all 3 particles, integration runs over all 3-momenta with a momentum conserving δ function, and index i is and the transverse Jacobi momentum κ ⊥ 12 is defined by With the specific choice of the functional form of the UV regulator in (10), the entire regulating factor in (11) is This concludes our definition of the bare regulated Hamiltonian.

C. Boost invariant similarity renormalization
The dependence of the Hamiltonian on the ultraviolet (UV) regulator ∆ is removed via similarity renormalization scheme whose basic idea has been outlined in the introductory part of this section. Similarity renormalization can be performed in various ways, in terms of matrix elements (the original G lazek-Wilson formulation) or in terms of coefficients of creations and annihilations operators (G lazek's similarity for particles), and for various definitions of the similarity formfactors in the effective Hamiltonian [18,20]. For a review and comparison of numerical efficiency of various formulations for Hamiltonian in QCD in (2+1) dimensions see [23]. The focus in this paper is restricted to the similarity renormalization group for particles (i.e. for creation and annihilation operators) that we use [19]. For introduction to the similarity renormalization group for particles see the latest reference in [18].
Let us denote the regulated canonical bare light-front Hamiltonian, H can ∆δ , together with the as of yet undetermined counterterms, X ∆ , in accordance with [19] as H ∆δ . H ∆δ is written in terms of creation/annihilation operators of bare, or current, quarks and gluons. For the purpose of this discussion, we denote all creation/annihilation operators, bosonic and fermionic, generically a † , a. If this current Hamiltonian were to be used to describe mesons, physical states would have to be very complicated superpositions of current quarks and gluons. We wish to find an effective Hamiltonian that provides a simpler description at hadronic scales λ. To arrive to a simpler picture of mesons, the effective Hamiltonian should rather be expressed in terms of constituent quarks and gluons. Yet, the two Hamiltonians are equal, meaning that they have the same eigenvalues. λ is a parameter of dimension of mass that will be explained in more detail below. As the Hamiltonians are related by the unitary transformation U λ , so are the current and constituent particles, The creation/annihilation operators a † , a satisfy the usual fermionic or bosonic commutation relations. We use Lorentz invariant normalization 2(2π) 3 p + δ 3 (p ′ − p) for the commutation relations.
Applying the transformation U λ to both sides of (15) and owing to its unitarity, one obtains From (16) follows that all dependence on λ in (17) is in the coefficients of the various combinations of creation/annihilation operators. So differentiating with respect to λ affects only the coefficients, not the operators, and we obtain where we have denoted in accordance with [19], and τ λ ≡ U † λ U ′ λ generates infinitesimal transformations. Equation (18) is an operator equation, but apart from c-numbers it contains only known, current, creation/annihilation operators and their commutators. For this reason, it can be solved in an operator form, in contrast to earlier formulations in terms of matrix elements.
The transformation is constructed from the requirement that the Hamiltonian H λ be band diagonal, which means, as we stated in the introductory part of this section, that it does not contain direct couplings between arbitrary scales. Parameter λ is a measure of width of momentum space formfactors f λ that appear in vertices of the renormalized H λ . For this reason, λ is often referred to as band width. Schematically, we require that which indicates that we want H λ to vanish in the part of the phase space were the similarity formfactor f λ is zero. This equation also defines G λ 2 . It is easier to work with G λ than with H λ 3 . G λ is split into the free part G 0λ = G 0 that is, in present work, independent of interactions and does not change with λ, and the interacting part G Iλ . Rewriting the differential eqn. (18) gives The left side of the equation is independent of G 0 because G 0 does not change with λ and because f λ does not change on the diagonal. G Iλ consists of the canonical interactions as well as new effective interactions ensuring that the effective Hamiltonian is equivalent to the original current one. Therefore, (21) describes two unknowns, τ λ and G Iλ . We have freedom to arbitrarily split this equation into two. Without loss of generality, we can assume that when the interactions vanish, so does the generator τ λ . If we choose and the generator of the transformation is of order G Iλ . Changes in G Iλ , corresponding to new effective terms, are of second order in interactions. Equation (23) can be solved for τ λ , for example order by order in the interactions. In fact, with G 0 being the free Hamiltonian (8), the solution is simple. It is the same operator as the right hand side of the equation (23) with an additional factor [ a E i − c E j ] −1 , where a, c stands for annihilated/created, and E's are free light-front energies created/annihilated in the vertex (see (8)). The additional factor arises from the free energy in G 0 and contracting one creation/annihilation operator with the conjugate in G 0 . It is, however, not necessary to explicitly find τ λ . If we denote the solution to (23) then the solution to (22) subject to the boundary condition where H ∆δ includes both known canonical terms and unknown counterterms, can be written as Counterterms X ∆ should remove any dependence on ∆ arising from the second term. This criteria is sufficient to determine X ∆ . Finite parts of the second term constitute new effective interactions and/or modifications to canonical terms. Finally, the Hamiltonian is found by substitution to (20). The equation (26) is an operator equation but the creation/annihilation operators are not affected by differentiation. They merely play a role for "counting" purposes. If we are interested in writing down a differential equation for a coefficient of a specific combination of creation/annihilation operators, to any given order in g we can see what operators in the interaction Hamiltonian need to be included on the right-side of the equation, and which of the creation/annihilation operators on the right side need to be contracted. Note that since the right side of (26) is proportional to a commutator, the coefficients in the effective Hamiltonian arise only from connected terms.
The preceding discussion was valid for an arbitrary choice of the similarity formfactor, and we have not specified how the formfactors are imposed. Specifically, we require that any combination of creation/annihilation operators in the Hamiltonian is accompanied by a formfactor f λ ; a particularly useful choice of the similarity formfactor f λ has the following form where M 2 c ≡ ( i created p i ) 2 is the square of sum of all free four-momenta created in the vertex and similarly, M 2 a ≡ j annih. p ′ j 2 is the square of all free four-momenta annihilated in the vertex. For example, for any finite λ Hamiltonian, the integrand of the operator (11) is modified by insertion of a factor The corresponding term in G λ , i.e. the Hamiltonian without the vertex formfactor, is to the lowest order in interactions just equal to (11) and independent of λ.
In this work, we expand the Hamiltonian relevant for a quark-antiquark state in powers of g, and find effective interactions up to g 2 . The perturbation theory for G λ is well explained in section 3.c of reference [19]. As it is often the case, in practice the most convenient way to determine the factors arising from the similarity transformation from infinity down to the finite scale λ in eqn. (26) is with the help of diagrammatic rules. Unfortunately, these seem to be the best kept secret in the field, and cannot be found listed in any reference, perhaps for fear of confusion with the ordinary perturbation theory for the S-matrix [18]. The calculation of effective interactions to O(g 2 ) is much simpler than at higher orders, because it involves only G λ at O(g) which is independent of the similarity scale. Consequently, all dependence on the scale in eqn. (26) is in f σ and d f σ /dσ. This is not true in general, when care must be taken to keep track of all σ dependence. Therefore, the second order calculation does not require "full-blown" diagrammatic rules and we will not develop them here. We will, however, show the simplified rules sufficient for the present calculation, with the understanding that they should and can be generalized for higher order applications.
To draw a diagram for the effective coefficient of any operator O λ in G λ , • draw an external line for each creation and each annihilation operator in O λ . Draw all connected diagrams connecting the external lines. (At second order, one-body operator has only one diagram, Fig. 1, reminiscent of a self-energy diagram in the ordinary perturbation theory. The two-body operator has two.) • Let us denote total momentum annihilated/created in a vertex i as a i , c i . (At second order, i = 1, 2.) Let p + i , p ⊥ i be the momentum conserved in the vertex i, and light-cone At each vertex i, include the standard factors, as well as regulators, already present in the G σ , such as, for example, the familiarū i ǫ µ γ µ u,with a momentum conserving delta function, at each qqg vertex i.
• The novel element comes from the similarity formfactors and scale. The entire diagram is to be integrated λ ∞ d σ. At second order, assign to the diagram. This factor contains all the dependence on the scale, and, for f σ given in (27), it can simply be integrated, resulting in • Diagrams can acquire an overall sign and/or combinatorial factors so that the same ordering of uncontracted creation/annihilation operators on both sides of the eqn. (26) is achieved [18]. At second order, however, these overall factors are obvious.
• Generalization of the procedure so that it is applicable for higher order calculations is straightforward.
As a specific example, we present in the next section the renormalization of the one body operator. Diagram for the quark one body operator at second order in g, not to be confused with the ordinary perturbation theory for S matrix. Thick lines correspond to uncontracted creation/annihilation operators. Double line is a contracted gluon, single thin line is a contracted quark. Vertical dashed lines show which momenta add up to a i , c i , respectively. All other notation is explained in the text.

III. RENORMALIZATION OF THE ONE-BODY OPERATOR
The one-body operator was previously derived, but not further studied, in ref. [18] with slightly different choices.
The one-body operator when λ = ∞ is where m 2 ∞ is the bare mass, ∆ is the UV regulator and δm 2 ∆ is the counterterm to be determined so that renormalized Hamiltonian is independent of the regulator.
When λ is finite, there are corrections (with coefficient δm 2 λ /p + ) to the one-body operator 4 which start at O (g 2 ).
To visualize the diagrammatic approach outlined at the end of previous section, we present in Fig. 1 a diagram for the lowest order correction to the one body operator with the notation used in the list of rules. For example, momentum a 1 is the free momentum of the uncontracted quark p; momentum c 1 in this diagram is equal to sum of the free momentum of the virtual gluon k g and the virtual quark k q . Three momenta are conserved, so only one of the virtual particles momenta is independent.
One factor of (11) (with appropriate momenta and without integration over the external momentum) stands for each vertex. Sums of the polarisation vector, spin and colors of the virtual particles lead to standard expressions involving color factors and matrix elements of D µν (k g )γ µ (/ k q + m) γ ν in the spinors of the uncontracted quark. It is convenient to express the virtual momenta in terms of Jacobi momenta κ ⊥ , x. This loop momentum is integrated over.
The similarity factor (30) further simplifies. Using that p + 1 = p + 2 = p + , and M 2 In order to write the one body operator with the same ordering as the original mass term, the uncontracted b, b † have to be commuted which gives an overall (-1). After this straightforward algebra, we find that the correction to the one body operator is indeed of the expected form constant/p + b † b (independent of p ⊥ ) with the coefficient independent of p + . Here m ≡ m ∞ ; z arises from integration over transverse momenta and carries a dimension of mass squared. The coefficient (33) consists of terms µ 2 ∆ independent of λ, containing UV/IR divergences, and λ-dependent terms, both finite µ 2 λ and IR divergent µ 2 λ, δ , where: Note that the coefficient of the quadratic divergence c δ ≡ 1 0 dy[r 2 δ (1 − y)/(1 − y) + 2r 2 δ (y)/y 2 ] diverges as δ → 0.
The one-body operator can be written as where The UV divergent µ 2 ∆ is combined with the counterterm δm 2 ∆ so that the dependence on the regulator ∆ is removed. Nevertheless, the mass m 2 λ is not finite in the limit δ → 0. In what follows, we are concerned only with the behavior in this limit, rather than details of the finite terms. For this reason and for simplicity, we grouped all finite terms contributing to m 2 λ into m 2 (λ).

IV. TWO-BODY OPERATORS
We are seeking a coefficient V of an operator consisting of two creation and two annihilation operators, which we expressed in terms of relative (Jacobi) momenta and set the total transverse momentum of the initial and final two-body states to zero because of kinematical boost invariance. We also factored out the total P + because all regulators and form factors are expressed in terms of boost invariant quantities. The conservation of the three momentum is implicit here and in what follows. The coefficient V λ depends on momenta, and contains terms of O(g 2 ) and higher. 5 We are interested only in the color singlet and most IR divergent part of this operator. (This part 5 We factored out the omnipresent 4 √ x 1 x 2 x 3 x 4 because it cancels with the similar factor in the definition of the bound state wavefunction. is diagonal in spins). Depending on its degree of divergence, it can be concluded whether or not the subleading terms are IR convergent. The canonical Hamiltonian already contains a two-body operator of this form, the instantaneous gluon exchange V 2 (9). With the infrared regularization as in [19] its coefficient is where r δ is the (unspecified) infrared regulator.
With the choice of f λ (27) and regularization as in ref. [19], the spin independent, most infrared divergent, color singlet part of the coefficient of the effective two-body operator to lowest (i.e. second) order in coupling is It arises from contracting one pair of creation and annihilation operators of a gluon with where x 5 = |x 3 − x 1 | corresponding to the two time-orderings. Invariant masses M are defined as at each vertex and r δ are the infrared regulators, in accordance with [19]. V 1 corresponds to particle 1 emitting the gluon that is absorbed by particle 2; in V 2 particle 2 emits gluon, particle 1 absorbs it. The second line in both expressions comes from integrating the width parameter in the similarity RG procedure from infinity down to its value λ.
In terms of the Jacobi momenta, the arguments of the formfactors are The coefficients of the effective one-body and two-body operators diverge when the infrared regulator is removed, i.e. r δ → 1. In previous calculations [17] the infrared divergences were found to cancel in the color singlet qq states, leaving however, the nonsinglet states with infinite mass. The remaining part of the two-body potential was found to be confining. A similar observation was made recently regarding QCD in (2+1) dimensions [23]. We want to check whether the cancelation occurs here.
It is easier to address the infrared divergence issue in the framework of a bound state equation, or equivalently, in terms of expectation values of the effective Hamiltonian. In that case, the creation and annihilation operators are contracted, and the integrands of the momentum integrals are just c-numbers consisting of the coefficients of the operators and wavefunctions of the bound state. The bound state eqn. reads where (4m 2 + 4m E) is the eigenvalue 6 , and Φ ij = Φ(x i , κ ⊥ ij ) is the bound state wave function of a quark-antiquark state of total momentum P ≡ P + , P ⊥ = 0 : The first line of the right side of the bound state eqn. (56) comes from one-body operators, the integral part arises from the two-body operator. The infrared structure of the expectation value of the two-body operator depends on the similarity formfactors in the various coefficients V. The coefficients also contain two different sets of infrared regulators corresponding to the two time orderings. To make the dependence on all formfactors explicit, let us introduce v 1 , v 2 as V inst does not contain any formfactor, but it does contain the two sets of infrared regulators for x 3 < x 1 and x 3 > x 1 . It is convienient to introduce v inst 1 , v inst 2 in analogy with v 1 , v 2 : This defines v's. Note that the infrared regulators in v 1 and v inst 1 are the same; so are in v 2 and v inst 2 .
To extract the infrared structure of the integral in (56), we first note that the IR divergence in the coefficients (45) and (47), (48) occurs when x 5 → 0 and that the small x 5 behavior of v inst i has opposite sign compared to v i . For this reason we first regroup the integrand in (56) by adding and subtracting v inst i multiplied by the same formfactors f λ as those multiplying v i in (58), (59). This splits the integral into two parts, where and Assuming that the wave function Φ 34 is bounded, it is straightforward to show that I (62) (schematically, the integral containing (1 − f 2 )(V inst − V 1 − V 2 )) is bounded and thus, has to be infrared convergent.
The remaining integral, I (63) is divergent as x 5 → 0. This means that x 3 → x 1 and x 4 → x 2 . Upon inspection of the integrand one can see that the formfactors in (63) restrict k ⊥ 5 more and more severely as x 5 → 0. It appears that, consequently, Φ 34 → Φ 12 as the divergence x 5 = 0 is approached. So we add and subtract Φ 12 to Φ 34 , i.e. replace Φ 34 in (63) by ([Φ 34 − Φ 12 ] + Φ 12 ). Then by expanding Φ 34 in Fourier series around Φ 12 it can be shown that the part of the integral with [Φ 34 − Φ 12 ] is, indeed, convergent as expected. Therefore, the infrared divergence is contained in multiplying Φ 12 . This is similar to the divergent term from the one-body operator.
To compare the divergence in the two-body operator with that of the one-body operator, we need to isolate it. It is sufficient to consider just series expansion for x 5 → 0. Keeping leading and subleading terms, arguments of the formfactors (55) reduce to The overall similarity formfactor is in this limit The formfactors contain factors k ⊥2 5 /x 5 which make the limit x 5 → 0 somewhat obscure. This can be evaded by a change of variables k ⊥2 5 /x 5 = K ⊥2 . In these variables, The subleading correction gives dx 5 ( √ x 5 ) −1 which is finite. The integral over K ⊥ can be done analytically. Out of all infrared regulators, only those regulating x 5 → 0 are necessary, because the integral only diverges at its lower limit. The remaining infrared regulators can be set to 1. The result is or, If this term multiplied by Φ 12 is subtracted from the integral in (56), it becomes finite in the limit δ → 0.

A. Cancelation of IR divergences?
Let us summarize the infrared divergent terms in the bound state equation (56). We have rewritten the bound state eqn. (56) as follows If the last line of (73) vanished, the bound state equation would be independent of δ, just as in previous formulation of the similarity renormalization [17], or, as recently reported, in (2+1) dimensions [23]. It is obvious that this does not happen here, except for the leading order in the nonrelativistic expansion. Indeed, in the nonrelativistic limit x i = 1/2, and the two different functions of xes, constituting the coefficients of the infrared divergence from the one-body and from the two-body operators, have the same value (4/ √ 2) at this point. However, an arbitrarily small deviation from x = 1/2 introduces a positive divergent constant into the bound state equation in the limit δ → 0. The bound state equation therefore is not defined in this limit.
For any finite δ, the bound state equation is well-defined and can be solved. If δ is small, the IR divergence forces the wavefunction to be peaked at x = 1/2. Detailed behavior of the wavefunction is likely to depend on the specific form of r δ . In any case, the wavefunction width determines the allowed range of x 5 . There is no such direct mechanism to restrict the transverse momenta, but they are pushed to small values by the formfactors in [V 1 + V 2 + V inst ] as x 5 is pushed to zero by r δ . It is peculiar that while the similarity formfactor f λ (M 2 12 −M 2 34 ) also restricts the range of transverse momenta, for x 5 = 0 this restriction is less severe than for any other possible x 5 . This means that the implicit formfactors are more important when x 5 → 0.
If we assume that with |ǫ| << 1/2, (consequently, x 5 << 1/2) and that other momenta, compared to mass, are also small, the bound state equation (73) can be expanded. Care must be taken, however, to keep track of the infrared divergence. We find HereẼ is the eigenvalue shifted by a constant from the one body operator,Ẽ ≡ E − (m 2 − m 2 (λ))/m. On the right hand side of the equation, the third term is a correction to nonrelativistic kinetic energy. It is likely that this term is small. The difference between the current mass and mass including finite λ-dependent corrections can be expected to be negligible for heavy quarks. The last expression in curly parenthesis is simply the subtraction of the IR divergence in the two body operator to this order in momentum expansion, so that the entire potential is final. The last line contains the residual divergence for x = 1/2 to its lowest nontrivial order, O (ǫ 2 ) . The structure of the right hand side becomes more clear upon substitution In these variables, The bound state equation (75) reads .
The last line in (81) is the uncanceled divergence which prevents us to take δ = 0 limit at the Hamiltonian level. Without numerically solving (81) we cannot make any quantitative statements, but some conclusions can be drawn just from the structure of the bound state equation. If not for the Coulomb potential, i.e. the first integral term in (81), it is obvious that with decreasing δ, the wavefunction would get squeezed more and more. In the limit, the wavefunction would basically collapse to the point x = 1/2 and zero relative transverse momentum. This might allow the eigenvalue to stay finite, and possibly roughly the same. With the Coulomb potential in place there could be some interplay between the positive δ divergence and the Coulomb. Note, however, that the Coulomb potential in our bound state equation is incomplete due to the presence of the exponential. The exponential cuts off the Coulomb potential at small q (and nonzero q z ). It is still possible (although not guaranteed) that the mass of the eigenstate is finite and convergent in the δ → 0 limit, but at the expense of an unphysical wavefunction.

VI. SUMMARY AND CONCLUSIONS
Owing to nonzero mass of quark/antiquark, the only infrared divergence in the one-and two-body operators in a quark plus antiquark state is due to x 5 → 0, i.e. the virtual gluon longitudinal fraction going to zero.
The infrared divergent part of the one-body operator obeys the one-body dispersion relation by construction of the similarity transformation. Specifically, it goes like a constant multiplying x −1 i , and is independent of the spectator. The expectation value in a quark/antiquark state is a simple sum of the quark and the antiquark contributions. But the IR divergence from the two-body operator is a different, more complicated function of xes. In addition to an overall explicit function, multiplying the divergent integral, whose form can be traced to the choice of f λ , there is an implicit dependence through the arguments of the regulators. Thus a different x-dependence can be generated for different regulators, contrary to the one-body operator. This occurs because boost invariance dictates the form of the one-body operator, and at the same time the two-body operator cannot lead to a function of the same form.
In the leading order nonrelativistic reduction, all dependence on infrared regulator drops out. Beyond nonrelativistic leading order the divergences do not cancel out in general, and it is not obvious that a regulator such that the divergence would cancel out for all x, exists. Caution is in place when the nonrelativistic limit is considered, because the infrared divergence, even though related to the most infrared divergent part of the effective terms, enters at the subleading order in the momentum expansion.
What are physical implications of this result and how to evade the problem? This is a subject of an ongoing study. Here we only outline some possibilities. We wish to emphasize that at this point they are just speculations.
The first question one needs to address is whether the IR divergences are to be cured at the level of Hamiltonian. In previous section we argued that it is possible to solve for the bound state with any fixed value of δ and at the end look how the solution behaves in δ → 0. The eigenvalue can be expected to be finite and convergent, however, the wavefunction in the limit must collapse to a point. Such a wavefunction is counterintuitive. It is not necessarily a no-go, however. Wavefunction by itself is observable only through expectation values of other operators. These operators have to be regularized consistently with the Hamiltonian. Then the behavior of observables need to be studied in the limit. This approach was adopted in a recent work on QCD (2+1) by Chakrabarti and Harindranath [23]. They solved the eigenvalue problem for a few formulations of similarity renormalization for Hamiltonian matrix elements. To avoid any confusion, we wish to emphasize that we use similarity renormalization for particles, and these two frameworks can have very different challenges. With this in mind, we note with interest that they find convergent eigenvalues, and some structure around x = 1/2 which they interprete as a sign of slow convergence.
It may turn out that removing the divergences at the level of Hamiltonian is necessary, or preferable. There are some compelling reasons for this. Leaving aside the open question of whether the limit of bound state solutions with finite δ converges in the limit δ → 0, improving our understanding (or lack of thereof) of the vicinity of the zero mode could be very useful in constructing more convenient similarity transformations. In present work we used as our diagonal part in the similarity transformation the current mass operator. Consequently, all formfactors are function of momenta compared to current mass. Even though in the bound state equation the kinetic energy is expressed in terms of effective, or constituent, λ-dependent mass, any expansions in momenta have to be performed in reference to current mass. This is certainly not an issue in case of heavy quarks, but for light quarks it practically rules out a simple nonrelativistic description of the bound state. It is not necessary, though, to use a λ-independent diagonal part to drive the similarity transformation. It could be devised around a λ-dependent, constituent mass, providing we have a resolution to the issue of the IR divergences at the Hamiltonian level.
One possibility is that understanding of the light-cone zero mode is indeed crucial for finding the hadronic spectra, and that the requirement that the physical observables are in compliance with Lorentz invariance is not sufficient to fully define the IR cutoff theory and render it regularization independent. If this is the case, then we are in trouble.
Another option that comes to mind is more down-to-earth. In a fully covariant theory, the physical state surely contains all Fock components. To maintain/restore symmetries order by order in coupling, it might be necessary to include in the state at least those Fock components that can mix with qq at the given order. To second order the only additional Fock component to include would be qqg. It is possible that the infrared divergence in the qq cancels with an infrared divergence in qqg. To support this view, recall, that in the standard perturbation theory the infrared divergence cancels out between instantaneous exchange and one gluon exchange. Even though our framework differs from the usual perturbation theory, a similar cancelation might occur. The reason is that the formfactor structure of the mixing between qq and qqg is the same as that of the effective qq terms which are plagued by the infrared divergence. The main difference is that the mixing depends on a three-body wavefunction. Nevertheless, it is conceivable that both the qq effective interaction and mixing with the three-body Fock component can be made finite by adding and subtracting the same IR counterterm, and it is possible that any ambiguities in connection with the choice of a specific IR regulator can be used to improve manifest rotational symmetry in the qq sector.