Solving the Bars-Green equation for moving mesons in two-dimensional QCD

The two-dimensional QCD in the large $N$ limit, generally referred to as the 't Hooft model, is numerically investigated in the axial gauge in a comprehensive manner. The corresponding Bethe-Salpeter equation for a bound $q\bar{q}$ pair, originally derived by Bars and Green in 1978, was first numerically tackled by Li and collaborators in late 1980s, yet only for the {\it stationary} mesons. In this paper, we make further progress by numerically solving the Bars-Green equation for {\it moving} mesons, ranging from the chiral pion to charmonium. By choosing several different quark masses, we computed the corresponding quark condensates, meson spectra and their decay constants for a variety of meson momenta, and found satisfactory agreement with their counterparts obtained using light-cone gauge, thus numerically verified the gauge and Poincar\'{e} invariance of the 't Hooft model. Moreover, we have explicitly confirmed that, as the meson gets more and more boosted, the large component of the Bars-Green wave function indeed approaches the corresponding 't Hooft light-cone wave function, while the small component of the wave function rapidly fades away.


I. INTRODUCTION
Two-dimensional QCD (hereafter QCD 2 ) has long served as a valuable toy model to mimic some essential dynamical features of strong interaction in the real world. A gratifying feature of this theory is that, due to the absence of transverse degree of freedom, the gluon is no longer a dynamical degree of freedom (at least in non-covariant gauges), but merely provides a linear color Coulomb potential, with the quark confinement as an almost trivial outcome. Despite this great simplification, QCD 2 still constitutes a rather nontrivial quantum field theory, which contains rich hadron phenomenology for mesons and baryons. There have been some numerical explorations of the QCD 2 with finite N based on the first-principle approaches, e.g., from the lattice Monte Carlo simulations [1,2] and from discretized lightcone quantization [3].
The 1/N expansion is a powerful and indispensable arsenal to tackle the nonperturbative dynamics of QCD [4][5][6]. As first exemplified in a 1974 seminal paper by 't Hooft [7], thanks to the dominance of the planar diagrams, QCD 2 in the large N limit (hereafter abbreviated as the 't Hooft model) indeed becomes much more tractable. The limit of infinite number of colors is in the following sense: where g is the strong coupling constant in QCD 2 , which carries the mass dimension one, and λ is often referred to as 't Hooft coupling constant. The first two conditions are standard large N assumptions. The last requirement, that the quark masses, m f , are much greater than the gauge coupling, is usually referred to as the weak coupling regime of the QCD 2 . The bulk of investigation on 't Hooft model has been mainly concentrating on this particular regime.
Employing the light-cone gauge A a + (x) = 0 and invoking the large-N limit, 't Hooft showed that mass spectra of an infinite tower of mesonic states can be inferred from the following integral equation: which is nothing but the light-cone Bethe-Salpeter equation for a relativistic qq bound state. m 1,2 are quark (antiquark) bare masses, M n denotes the mass of the n-th excited mesonic state (n is the principal quantum number to characterize a meson living on a string), and φ (n) (x) signifies the corresponding light-cone wave function, with x ∈ [0, 1] representing the fraction of the light-cone momentum carried by the quark with respect to that by the meson. The symbol− implies that a principle-value prescription is exerted to eliminate the infrared divergence that occurs at y → x. In general, Eq. (2) is not admissible to an analytic solution, yet can only be solved using numerical recipes. Intriguingly, the discrete meson spectrum determined from (2) exhibits the Regge trajectory. Specifically speaking, for highly excited states (n ≫ 1), the squared meson masses are well described by M 2 n = 2π 2 λn + m 2 1 + m 2 2 − 4λ ln n + O(n 0 ).
In 't Hooft's pioneering work, some essential properties that closely resemble the ordinary QCD, such as color confinement and Regge trajectory, have already been revealed. Afterwards there have been extensive investigations on various "phenomenological" aspects of the large-N QCD 2 , e.g., hadron decay/scattering amplitudes, current correlators, form factors, (naive) asymptotic freedom and parton model, fragmentation functions, Pomeron, (generalized) parton distribution functions, quark-hadron duality, and many more else [8][9][10][11][12][13]. Apart from these work, a very remarkable feature of this model has also been uncovered in the mid 1980s: the non-vanishing quark condensate, the spontaneous breaking of chiral symmetry, and the clarification of (quasi-)Goldstone mode [14][15][16].
Historically, most aforementioned features of the 't Hooft model have been deduced by utilizing the light-cone gauge (often peppered with the light front (LF) quantization). The greatest advantage of this procedure is that, it allows to yield the manifest boost-invariant bound-state equation, (2), and generates compact expressions for various physical quantities. Nevertheless, there also exist some disadvantages inherent to this approach, e.g., the mechanism initiating chiral symmetry breaking becomes obscure, due to the perturbative nature of the vacuum in LF quantization.
Interestingly, there also exists an alternative perspective to tackle QCD 2 , that is, by imposing the axial gauge (A a 1 (x) = 0) condition in the ordinary equal-time quantization. Despite its technical complication, this approach does possess some notable virtues. For example, unlike the light front quantization, the equal-time quantization accommodates a nontrivial vacuum state, which makes the study of the spontaneous chiral symmetry breaking much more transparent. Moreover, this approach has a natural connection to the familiar constituent quark model, in analogy with the intimate connection between the light-cone gauge and the parton model.
In 1978, by quantizing 't Hooft model in the axial gauge, Bars and Green presented a formal proof that Poincaré algebra does close in the color-singlet channel [17]. They further derived the axial gauge Bethe-Salpeter equations for mesons in an arbitrary frame. The resulting relativistic bound-state equation (hereafter dubbed Bars-Green equation) does look much more sophisticated than its light-cone counterpart, Eq. (2). As is widely known, it is highly nontrivial to conduct the Lorentz boost for a bound-state wave function constructed in the equal-time quantization procedure [18]. It was anticipated that Bars-Green equation must preserve Poincaré invariance, i.e., the meson mass should not rely at all on which Lorentz frame one is carrying out the measurement, which is clearly a rudimentary requirement for any sensible theory. A special and gratifying situation is when a meson is viewed in the infinite momentum frame (IMF), the Bars-Green equation can be proven to exactly reduce to the 't Hooft equation [17]. Nevertheless, it remains an analytic challenge to prove that Bars-Green equation does preserve Poincaré invariance in any finite momentum reference frame.
In general, it is impossible to solve the Bars-Green equations in an analytic fashion. The numerical investigation of these equations were pioneered by Li and collaborators in late 1980s [19], but only for the mesons in the zero-momentum frame. They indeed confirmed that the calculated meson spectra using axial gauge agreed with what were found by solving the 't Hooft equation.
The aim of this work is to extend the earlier investigation in [19], by numerically solving the Bars-Green equations for a generic moving meson, with hadron species ranging from the chiral pion to heavy quarkonium. Our primary goal is to numerically validate the Poincaré invariance of the Bars-Green equations. Moreover, we wish to quantitatively assert that, to which extent when a meson gets boosted, the Bars-Green wave function would resemble the corresponding 't Hooft wave function to a decent degree.
It is worth mentioning that, the Bars-Green wave function is intimately related to the so-called quasi-distributions in QCD 4 , which have received lots of attention in past few years. The quasi-distributions, a set of instantaneous yet spatially non-local correlators, was recently introduced by Ji as a proxy to help extract the light-cone distributions, which can be directly computed by lattice simulation [20]. One of the key properties of the quasidistributions is that, when boosted to the IMF, they will reduce to their light-cone cousins, e.g., parton distribution functions and light-cone distribution amplitudes. We wish that our comprehensive numerical study of Bars-Green wave functions will lend some guidance on quantitatively understanding of the properties of quasi-distributions in realistic QCD. The paper is structured as the following. In Section II, we present a relatively succinct, yet self-contained review on the course of arriving at the Bars-Green equation in the 't Hooft model. We also add more details in illustrating how does the Bars-Green equation reduce to 't Hooft equation in the IMF. In Section III, we investigate the renormalized chiral condensates with a variety of quark masses in the axial gauge, and compare with the respective values obtained in the light-cone gauge. We also present the analytical formulas for the decay constants of the even-and odd-parity mesons. In Section IV, we briefly describe our numerical strategies in solving the 't Hooft equation and Bars-Green equations. In Section V, we then present comprehensive numerical studies of the mass spectra and Bars-Green wave functions for a variety of meson species: chiral π, physical pion, a fictitious strangeonium, and charmonium, for each of which several different meson momenta are chosen. We also examine how fast the Bars-Green wave function for a highly boosted meson converges to the respective 't Hooft wave function. We conclude this Section by examining the frame-independence of the decay constants. Finally we summarize in Section VI.

II. REVIEW OF THE BARS-GREEN FORMALISM
In this section, our main goal is to sketch some key intermediate steps in deriving the Bars-Green equation in axial gauge. Nothing in this section is really new, and the purpose of including this section is mainly for the sake of completeness. We will follow the Feynman diagramatic approach to derive the BG equation. It is worth noting that, there also exists an elegant alternative way to derive the same equations from the Hamiltonian approach.
We start from the Lagrangian of the 1+1-dimensional QCD with the SU(N) color gauge symmetry: where the gluon field strength F a with T a being the color SU(N) generator in fundamental representation and a running from 1 to N 2 −1. We adopt the Dirac-Pauli representation for the γ-matrices: Throughout this work, we are imposing the axial gauge condition A a 1 (x) = 0 1 . Like in the light-cone gauge, the nonlinear term in F a 01 , the major characteristic complication of QCD, simply drops out in the axial gauge. Moreover, A a 0 is no longer a dynamical variable, instead can be expressed in term of a quark current through the Euler-Lagrange equation. Hence, in the canonical Hamiltonian form, the gluon field A a 0 (x) has been completely eliminated, whose effects are fully encoded in the instantaneous, yet spatially-nonlocal current-current interaction. As a common practice, the current-current interaction can often be simulated by a gluon propagator in the axial gauge: with the only survivor from the 00-component, and x µ = (x 0 , x 1 ). It can be immediately identified with the instantaneous linear Coulomb potential. It is instructive to rewrite it as a Fourier integral: The momentum-space gluon propagator only depends upon the spatial component of k µ , reflecting that A a 0 in the axial gauge is a non-propagating degree of freedom. Moreover, due to the singular behavior of the integrand near k 1 → 0, one must introduce a proper prescription to make the above Fourier integral well-defined. This is the origin of the ubiquitous occurrence of the principle value prescription in two-dimensional gauge theory. Let S(p µ ) denote the full (dressed) quark propagator, Σ(p µ ) signify the 1PI quark selfenergy. In the large N (planarity) limit, the standard rainbow approximation in the Dyson-Schwinger equation for the quark self-energy, as pictorially depicted in Fig. 1, becomes a rigorous procedure:

A. Mass-Gap Equation
where− implies a principal-value prescription. It turns out that the quark self-energy Σ only depends on the spatial component of the two-vector p µ , and can be parameterized as For notational brevity, from now on we often use the symbol p to represent p 1 , unless otherwise explicitly stated. It is convenient to introduce two new variables E(p) and θ(p), in replacement of the the functions A(p) and B(p): where E(p) characterizes the energy dispersion of the dressed quark, and θ(p) is usually referred to as the Bogoliubov (chiral) angle. Integrating (7b) over k 0 , and after some algebra, one arrives at the so-called mass-gap equation [17]: θ(p) is an odd, and, monotonically rising function in p, which approaches ± π 2 as p → ±∞, respectively. It turns out to be an impossible mission to express θ(p) in terms of the known special functions. As matter of fact, this nonlinear integral equation can only be solved numerically, even in the chiral limit. Notice in the free theory (λ = 0) limit, θ(p) = tan −1 (p/m), recovers the familiar Foldy-Wouthuysen angle in the free Dirac theory.
Interestingly, Eq. (10) can also be derived from an alternative perspective, viz, by the requirement of minimizing the vacuum energy 2 .
Once the θ(p) is determined, one can proceed to determine the dispersive law of the dressed quark [17]: which is clearly an even function of p. Note this dispersive relation is not even Lorentz covariant. This can be attributed to the fact that, since the Poincaré algebra does not close in colored sector, so the Lorentz covariance is scarified in a single quark sector, though it must hold in color-singlet channel. As will be elaborated in Section IV, we solve the mass gap equation (10) numerically using Newton method. In Fig. 2, we plot the profiles of the chiral angle and dispersive law as function of quark momentum, for several different quark mass. We see that the Bogoliubov angle in chiral limit indeed assumes a nontrivial shape. For small bare quark mass, when the quark momentum is very soft, the dressed quark energy may even become negative. This pathological behavior can be readily seen from Fig. 2, which can also be understood from the approximate formula

B. The Bars-Green Equation
As a confining theory, QCD 2 admits an infinite tower of stable color-singlet mesons in the large N limit. We are interested in inferring the bound-state equation from the familiar Bethe-Salpeter approach, although the alternative approach, i.e., the Hamiltionian method, may be particulary illuminating in certain aspects. For simplicity, throughout this work we have focused on the flavor-neutral quarkonium state, though the extension to the flavored mesons are straightforward. The meson-quark-antiquark vertex, denoted by Γ(p µ , P µ ), obeys the homogenous Bethe-Salpeter equation: where P µ is the two-momentum of the meson, and p µ (P µ − p µ ) is the momentum of the external quark(antiquark) leg. This equation is pictorially represented in Fig. 3, where the ladder approximation also becomes justified, thanks to the planarity condition. It is a standard practice to introduce the Bethe-Salpeter wave function: Before proceeding, it is convenient to rewrite the full quark propagator in (7a) as It is also useful to decompose the Bethe-Salpeter matrix wave function Φ into a pair of wave functions φ ± : Substituting (14) into (12), and integrating both sides over p 0 by employing the method of residues, one ends up with two coupled equations for each mesonic state with mass M n (of the n th mesonic level) [17]: where P µ P µ = M 2 n , the dressed quark energy E(p) is given in (11), and where the Bogoliubov angle θ(p) is deduced from solving the gap equation (10). Eq. (16) is the mesonic bound-state equation in axial gauge with equal-time quantization, hereafter referred to as the Bars-Green (BG) equation. It is the instant-form counterpart of the 't Hooft equation, (2). Unlike a single meson wave function φ in light-front formalism, here one must introduce a pair of meson wave functions φ ± , in order to warrant the Lorentz covariance in the equal-time quantization. The much more sophisticated form of the BG equation with respect to the 't Hooft equation, simply reflects the widely-spread tenet, that boosting the equal-time bound-state wave function is a highly nontrivial mission, in sharp contrast to the boost-invariant LF formulation.
The wave functions φ ± , representing the large (small) component of bound-state solution, respectively, characterize the probability amplitude for the qq pair moving forward (backward) in time. Their physical meanings become even more transparent in the bosonization approach, where the φ ± can be directly interpreted as the coefficient functions associated with a Bogoliubov transformation from the composite quark-antiquark creation operator to the mesonic creation operator [21]. Specifically speaking, there are two ways to produce a mesonic state from the vacuum in the equal-time quantization. One can always produce a meson state by creating a qq pair out of the vacuum, regardless of the (non)trivial nature of the vacuum, with the probability amplitude characterized by φ + . Nevertheless, when the vacuum is nontrivial, e.g., which may accommodate a nonzero quark condensate, as what is happening in the equal-time formulation for large N QCD 2 , one can also create a meson by removing a redundant pair of qq from some correlated quark-antiquark pairs constantly fluctuating out of the nontrivial vacuum. It is intuitively conceivable that, the relative magnitude of φ − with respect to φ + gets more and more suppressed with the increasing quark mass/meson momentum/principal quantum number.
The meson wave functions φ ± obey the following orthogonality and completeness condi- The ubiquitous minus signs are reminiscent of the nature of Bogoliubov transformation, which are manifest using the Hamilton approach 4 . It is instructive to examine the properties of the wave functions under discrete symmetry transformation. Since QCD 2 is symmetric under space inversion, charge conjugation, the flavor-neutral quarkonium wave functions must be subject to the following relations: where η n = (−) n+1 signals the intrinsic parity of each meson 5 . Thus all the mesonic levels simply alternate in parity: parity-odd states (n = 0, 2, 4, . . .) and parity-even states (n = 1, 3, 5, . . .). These symmetry relations signal a notable virtue of Bars-Green formalism versus LF formulation, since it is far from straightforward to realize the parity transformation in the latter setup.
A very special case is the ground-state meson in the chiral limit, which turns out to be parity-odd and exactly massless. For this reason, it is often dubbed chiral pion, π χ . Its BG wave functions are known in a semi-analytical fashion [21]: where θ(p) is the corresponding Bogoliubov angle in the chiral limit. Applying (19a) to (20), we can obtain the BG wave functions of the chiral pion that moves toward the negative x-axis: Since one cannot boost a massless particle to its rest frame, some sort of irregularity is anticipated to occur in the P → 0 limit. For a fixed p, the BG wave functions for a chiral pion, φ πχ ± (p, P ), are continuous yet nonanalytic across the point P = 0, as indicated in (20) and (21).

C. Bars-Green equation in IMF
Examining the the coupled integral equations (16), it is by no means transparent to prove the meson mass spectra are frame-independent. Nevertheless, Bars and Green have argued that, in the IMF, the backward-moving wave function φ − must diminish, so the BG functions must reduce to the celebrated 't Hooft equation, consequently the forward-moving wave function φ + can be identified with the 't Hooft wave function φ(x). Bars and Green have already outlined all the necessary clues for the proof. Nevertheless, for the sake of completeness and clarity, we decide to supplement more technical details in intermediate steps, together with some pictorial evidences, to validate Bars and Green's claim.
Let a flavor-neutral meson carry the nonzero momentum P > 0. Let us first introduce a pair of dimensionless ratios x, y, by x = p/P and y = k/P , where p and k represent the quark momenta appearing in (16). We subsequently reexpress the Bars-Green wave function as φ ± (x, P ) ≡ φ ± (p = xP, P ). At this stage, the range of x and y remains unbounded. Let us temporarily assume x, y are not in proximity to 0. In the IMF limit P → ∞, the Bogoliubov angle θ(xP ) is then approaching its asymptotic values, π 2 ǫ(x). Consequently, one finds in the P → +∞ limit, where Θ designates the Heaviside step function. Therefore, the C function equals 1 if 0 < x, y < 1, or x, y < 0, or x, y > 1, and vanishes in all other cases; The S function always vanishes except when x < 0, y > 1 or x > 1, y < 0, in which cases it equals -1.
For the sake of clarity, we take the chiral limit as a concrete example, with the respective Bogoliubov angle θ(p) shown in Fig. 2. We then generate the corresponding C(x, y, P ) and S(x, y, P ) functions, with several different choice of P . From Fig. 4, we clearly see the trend that, in the IMF limit, C(x, y, P ) and S(x, y, P ) indeed exhibit the behavior as dictated in (22). Following Ref. [17], it is straightforward to see that the right-hand sides of Bars-Green equations (16) must scale as O(1/P ) in the IMF limit. Therefore, any term in the left-hand side which scales as P 1 or P 0 must cancel, and the O(1/P ) terms in both sides must be matched. By examining the asymptotic behavior for the factor E(p) + E(P − p) ± P 0 , we thereby find that φ − (x, P ) must vanish for all x, and φ + (x) is non-vanishing only when 0 ≤ x ≤ 1. Notice than, for 0 ≤ x ≤ 1, Eq. (22) then implies that C(x, y, P → ∞) → Θ(y(1 − y)), and, S(x, y, P → ∞) → 0.
From the gap equation (10), it is easy to see that tan θ(xP ) = xP m +O(1/P ) in the P → ∞ limit. Therefore, the dispersive law in (11) in the IMF simplifies into Employing the aforementioned simplifications in the IMF, the Bars-Green equation (16) then reduces to Approximating P 0 = √ P 2 + M 2 by P + M 2 /2P , and matching both sides of (24) through the linear order in 1/P , one finds that As promised, the BG equation (16) for φ + (x, P ) in IMF does reduce to the 't Hooft equation (2) (with m 1 = m 2 = m), while φ − (x, P ) dies away with a pace ∝ 1 P 2 . In the following sections, we will numerically examine the tendency of the BG wave functions φ ± with an ever-increasing meson momentum.

III. SOME LORENTZ-INVARIANT QUANTITIES IN AXIAL GAUGE
There are some basic yet important nonperturative quantities, exemplified by the quark vacuum condensate (for arbitrary quark mass) and meson decay constants, which have been extensively studied in the LF formulation of QCD 2 . In this section, we revisit these quantities in the axial gauge in equal-time quantization. To our knowledge, the studies from the perspective of Bars-Green formalism are novel. The purpose of this Section is to make a nontrivial examination of the gauge and Lorentz invariance (frame independence) of these simple QCD matrix elements.

A. Quark condensate
Since the mid-80s, it became widely known that the 1 + 1-dimensional QCD in the large N limit actually accommodates spontaneous chiral symmetry breaking (SCSB), signalled by the non-zero quark condensate [15,16].
In passing, it is worthwhile to elaborate on the possible phases in the 't Hooft model. The massless QCD 2 (N → ∞) can be classified in two distinct regimes, depending on the order of taking the N → ∞ and the chiral limit, which turn out not to commute [15]: 1) In the weak coupling regime, one assumes m q ≫ g ∼ 1 √ N , and the N → ∞ limit is taken prior to ultimately sending m q → 0. This phase corresponds to the familiar mass spectrum from solving 't Hooft equation, where the spontaneous chiral symmetry breaking occurs. 2) In the strong coupling regime, where m ≪ g ∼ 1 √ N is instead assumed, and one first takes the chiral limit, then followed by sending N → ∞. Chiral symmetry remains unbroken in this phase, and the corresponding spectrum is rather different, where there appear massless composite fermion rather than the massless meson [22][23][24][25][26]. In this work, we have tacitly assumed to exclusively consider the weak-coupling phase.
At first sight, the occurrence of SCSB in QCD 2 (with N → ∞) appears to contradict Coleman's theorem [27], which seems to rule out the possibility of spontaneous breakdown of any continuous symmetry in two dimensional field theory. This puzzle was first resolved by Witten [14] in the context of the SU(N) Thirring model in the N → ∞ limit (this model is also commonly referred to as the Gross-Neveu model). He pointed out that the Berezinskii-Kosterlitz-Thouless (BKT) phenomenon actually occurs in this case [28,29], so that the chiral symmetry is "almost" spontaneously broken. Later Zhitnitsky realized that, in the weak coupling regime, the 't Hooft Model also exhibits exactly the same BKT effect [15], so that the SCSB also occurs in the N → ∞ limit. The spontaneous chiral symmetry breaking is consistent with the spectrum of 't Hooft model that the mesonic states with opposite P -parities are non-degenerate in the masses.
Specifically speaking, one can show that the following two-point correlator in the 't Hooft model possesses the following large-|x| behavior [15]: In the N → ∞ limit, the correlator approach a non-vanishing constant, thus exhibiting the true long-range order, and heralding the occurrence of the massless boson mode; for any large but finite N, the correlator falls off very slowly with x, and there does not arise massless meson. Hence there is no contradiction with Coleman's theorem. Despite the notion of pertubative vacuum in the light-front quantization, a nonvanishing chiral condensate was first discovered from this formalism. Using the operator expansion technique, Zhitnistsky has found an analytic result for the vacuum quark condensate in the chiral limit in 't Hooft model [15]: Since the condensate is nonanalytic in the 't Hooft coupling λ, it characterizes a type of nonperturbative effect that cannot be captured by summing perturbation series in λ n (n being a non-negative integer) to all orders. Later, Burkardt has extended Zhinitstsky's analysis, and presented an analytic formula also for massive quark [30]. After subtracting the logarithmic UV divergence, he obtains the renormalized quark condensate for an arbitrary value of m: where α = 2λ/m 2 , γ E = 0.5772 . . . is the Euler constant, and .
A nontrivial vacuum state naturally emerges in QCD 2 if equal-time quantization is taken. It was first by Li [16] who first reported a nonzero quark condensate in the chiral limit in the axial gauge: Tr Substituting the numerical solution of θ(p) from (10), into this equation, one readily verifies (27) obtained from light-front formalism, to a high numerical accuracy. For nonzero quark mass, the integral in (30) becomes logarithmically UV divergent. Subtracting the analogous term arising from cosine of the Foldy-Wouthysen angle of a free quark, which amounts to performing an additive renormalization, one finds that the renormalized quark condensate in axial gauge is In Fig. 5, we plot the the quark condensate as function of quark mass, stemming from the light cone gauge, (28), and the axial gauge, (31). Quite satisfactory agreement is achieved, firmly establishing the gauge invariance of the quark condensate.

B. Decay constants
One can define the meson decay constant f (n) as where ǫ µν is the antisymmetric Levi-Civita tensor in two dimensions.
In the light-cone gauge, Callan, Gross and Coote [8] were able to identify the decay constant for the n-th mesonic level simply with the integral of the 't Hooft wave function: For the chiral pion, π χ (the massless parity-odd state affiliated with m = 0) 6 , the 't Hooft wave function possesses a peculiar form: φ πχ (x) = Θ(x(1 − x)), so the decay constant simply is The decay constant in the axial gauge can be most conveniently worked out following the Hamiltonian method [21]. With the aid of the bosonization technique, one can reexpress the axial vector current in term of meson's creation and annihilation operators, and readily ascertain the intended decay constant.
We separately discuss the decay constants of mesons with odd and even parity, as designated in (32). First we consider the mesonic level with even n (odd parity): These two expressions for the decay constant are obtained by utilizing the different axial vector Lorentz index in (32). Although both analytical expressions superficially differ, they are doomed to be equal by Lorentz invariance. Furthermore, although these expressions explicitly depend on the meson momentum P , the frame-independence of the decay constant enforces some identities that θ(p) must obey. In Section V, we shall present explicit numerical evidences for the frame/Lorentz-index independence of the meson decay constants. Note it is quite delicate to extract the decay constant for a stationary (P → 0) meson from (35). We emphasize that the inclusion of the small component of the BG wave function, φ − , is crucial to warrant the frame-independence of the decay constant. It is interesting to examine the decay constant of the chiral pion in the axial gauge. Substituting the analytic BG wave functions (20) into (35), for a pion carrying arbitrary positive momentum P , we find Though far from obvious to see why the integral is exactly equal to P , it has to be so to match the LF result for chiral pion, Eq. (34). Next we turn to the mesonic levels with odd n (even parity). For flavor-neutral mesons, such states have odd C parity , so the corresponding 't Hooft wave functions are odd in exchanging x and 1 − x. As a result, the decay constants simply vanish in line with the prediction from the light-cone gauge, (33).
Notwithstanding this trivially looking result, it is still instructive to examine these decay constants from the angle of axial gauge. The respective decay constants in this case read Again we show the expressions extracted from (32) by utilizing two different axial vector Lorentz indices. Making use of the fact that θ(p) is an odd function of p, and the odd C-parity of the BG wave functions for the odd-n states as encoded in (19b), one can prove that the integrals in (41) indeed vanish, for all possible meson momentum.

IV. NUMERICAL RECIPES FOR SOLVING BOUND-STATE EQUATION
Numerically solving 't Hooft equation has gained a mature status, so here we just briefly describe the numerical strategies adopted in this work. This type of equation is usually solved by the spectrum method, with the solution presumed to be a linear combination of a set of basis functions. For massive quark (m ≫ √ 2λ), it proves convenient to invoke the so-called Multhopp method, which utilize the trigonometric basis functions [10] [31]. For light quark (m ≤ √ 2λ), yet it is more advantageous to follow 't Hooft's original method [7], that adopts a set of basis functions such as Ψ n (x) = Ax β 1 (1 − x) 2−β 1 + Bx β 2 (1 − x) 2−β 2 + n C n sin(nπx), in which parameters β 1,2 are determined by the boundary conditions πβ 1,2 cot(πβ 1,2 ) = 1 − m 2 /2λ. Empirically, n ∼ O (10 1 ) is sufficient to yield stable first three digits.
Prior to solving the Bars-Green equations, one has to first determine the chiral angle θ(p) to a decent accuracy. Here we follow Ref. [19] to use the generalized Newton method. It is convenient to first change the variable from p to ξ using p = √ 2λ tan ξ, so ξ ∈ (− π 2 , π 2 ), within a finite interval. The mass gap equation in (10) is then discretized to a set of matrix equations: Suppose N is a prescribed large positive integer, designating the size of a grid. Both variables ξ j , η j = jπ 2N are evenly partitioned on the grid, with j being an integer obeying −N ≤ j ≤ N. In addition, the boundary conditions θ(± π 2 ) = ± π 2 must be imposed. Greater N will generally decrease the discretization errors, nevertheless render the computation more expensive. Practically, N = 100 works well for the case of light flavors, i.e. π and ss mesons. A finer grid with N = 300 or higher is required for heavy mesons, i.e., for cc quarkonium. The numerical results of the chiral angle θ(p) and dispersive function E(p) have been shown in Fig. 2.
Analogous to the 't Hooft equation, Bars-Green equation can be solved by means of the spectrum method as well. The major complication is due to the emergence of the additional φ − (p, P ) component, so one inevitably confronts coupled integral equations. In the late 80s, Li et al. solved the BG equations for a variety of stationary flavor-neutral mesons, choosing the basis functions as the quantum harmonic oscillator's eigen-functions [19].
As described in Section II C, it is convenient to introduce a momentum fraction variable x = p/P for a flavor-neutral meson, with meson momentum denoted by P and quark momentum represented by p. The Bars-Green wave function is then effectively expressed as φ ± (x, P ) ≡ φ ± (p = xP, P ). In contrast to the light-front momentum fraction x ∈ [0, 1] in the 't Hooft wave function, the range of x in BG wave functions is completely unbounded. The normalization condition in (18a) for the BG wave functions can be rewritten aŝ To apply the spectrum method to a moving meson, we generalize Li et al.'s orthogonal basis functions as follows 7 : where H m represents the m-th Hermite polynomial and α is a variational parameter that can be tuned to minimize the mass of the ground state.
φ ± (x, P ) = N −1 m=0 a ± m Ψ 2m (α, x, P ), n even; Solving the original coupled integral equations are then transformed into the matrix eigenvalue problem. After diagonalization of the N × N matrix, one can determine the mass spectra of the first N same-parity mesonic states from the discrete energy eigenstates, M n = (P 0 n ) 2 − P 2 , as well as the corresponding BG wave functions φ (n) ± (x, P ). Practically, for most cases taking N ≈ 20 appears to be adequate.
Before concluding this section, we describe the principal-value prescription employed in this work in solving 't Hooft and Bars-Green equations. To tame severe infrared divergences, two distinct strategies are implemented: where f (y) is a smooth test function which is regular at y = x. The first recipe is the subtraction method utilized in [19], while the second is the Hadamard regularization for hypersingular integral [33]. In practice, both prescriptions yield stable and convergent results.

V. NUMERICAL RESULTS
Being a super-renormalizable theory, QCD 2 bears the gauge coupling g with unit mass dimension. In the large N limit, we set the absolute mass scale following the ansatz in Ref. [12], e.g., choosing the value of the 't Hooft coupling λ such that πλ = 0.18 GeV 2 , in conformity to the value of string tension in the realistic QCD 4 . For notational brevity, we will express any dimensionful quantity in units of √ 2λ = 340 MeV in the rest of this section. In the hypothetic 1+1-dimensional world, we attempt to mimic realistic mesons in QCD 4 as much as possible. In the √ 2λ unit, the masses of physical π and J/ψ mesons are M π = 0.41, and M J/ψ = 9.03, respectively. Solving the 't Hooft function for ground state, the corresponding quark mass are found to be m u = 0.045 8 , m c = 4.23.
We have also intentionally invented a fictitious strange quark with m s = 0.749. As can be seen in Fig. 2, it corresponds to a peculiar threshold θ(ξ) (with ξ = tan −1 p), through which the profile of θ(ξ) passes from the convex to concave with the increasing quark mass. This particular strange quark mass is determined by minimizing the relative distance between the θ(ξ) and the straight line θ(ξ) = ξ. The lowest-lying strangeonium state has M ss = 2.18. Naively speaking, the strange quark with m s = 0.749 might be thought of severing a threshold, below which is called light flavor, and above which is called heavy flavor. For completeness, we also consider the chiral limit with m u = 0, which can host a lowestlying massless state named chiral pion. The mass spectra of the first low-lying mesonic levels, for each quark mass, are listed in Table I as well as in Fig. 6. As is well known, the excited states fit into the linear Regge trajectories to a good degree. We stress that, the mass spectra found by solving the BG equation appear to be frame-independent, and always agree with what are obtained by solving 't Hooft equation. Thus these numerical studies constitute a nontrivial validation of Poincaré invariance in the Bars-Green formalism.  In passing, one might like to take a closer look at the spontaneous chiral symmetry breaking in QCD 2 . The celebrated Gell-Mann-Oakes-Renner relation states that With the aid of (34) and (27)  (126 MeV), which is quite close to the input pion mass M π = 0.41 (139 MeV). Thus, the pseudo-Goldstone nature of the "physical" pion is explicitly validated.
We proceed to show the profiles of various bound-state wave functions. In Fig. 7, we first plot a number of 't Hooft (LF) wave functions for the ground state and the first excited state in the mass spectra, affiliated with the different quark species. As dictated by charge conjugation symmetry, the LF wave functions with even/odd n are symmetric/antisymmetric under the exchange x ↔ 1 − x. The LF wave functions always vanish in both end points x = 0, 1. For lighter quark, the LF wave functions for ground states exhibit a very steep rising/falling behavior when x approaches the boundaries, and a stable plateau in the majority of range in x (Note the slope in the boundaries becomes infinite in the chiral limit!). For heavy quark, the LF wave function possesses a much milder rising/falling shape near the end points, and the plateau disappears.
In Fig. 8, 9, 10 and 11, we plot various Bars-Green wave functions of the ground state and the first excited state associated with several quark flavors: massless u, m u = 0.045, m s = 0.75, and m c = 4.23, respectively. For the sake of comparison, we also juxtapose the LF wave functions of the corresponding meson in each figure. The BG wave functions with even/odd n are symmetric/antisymmetric under the exchange x ↔ 1 − x, as required by the charge conjugation symmetry for flavor-neutral states.
Obviously, the BG wave functions are spatially much more spread in x-axis than the 't Hooft wave functions, which are confined within the interval x ∈ [0, 1]. From these figures, one clearly observes that, for all species of quark flavors, when the mesons are boosted with higher and higher momentum, the φ + component of the BG wave functions will approach the 't Hooft wave functions, while the φ − components rapidly dies off. These behaviors are completely compatible with the anticipated IMF limit of the BG wave functions in (25).
In most cases, the "backward-motion" wave functions φ − (x, P ) are always much less significant in magnitude than the "forward-motion" components φ + (x, P ). The only exception is a "wee" (very low-momentum) lowest-lying meson, as exemplified by chiral and physical pions in Fig. 8 and Fig. 9. As mentioned before, the comparable magnitude between φ + and φ − is expected for the soft Goldstone boson, which can be intuitively attributed to the nontrivial vacuum structure, characterized by the nonzero quark condensate. The φ − Figure 8: Bars-Green wave functions for the low-lying uū states in the chiral limit: chiral pion π χ (first row), and the first excited state (second row). component becomes quickly suppressed with respect to φ + , provided that the meson momentum increases, or going to higher excited states, or increases the quark mass, which can be attributed to the rapidly decreasing energy denominator 1/(E(p) + E(P − p) + P 0 ) in (16). This is somewhat analogous to the case of the Dirac equation, where the disparity between the large component and small component of the Dirac spinor becomes substantial when going to nonrelativistic/ultra-relativistic limit.
We would also like to mention a technical nuisance. As can be seen in Fig. 11, some wiggles have emerged in φ − (x, P ) for the lowest-lying and first excited charmonium state. This should be regarded as the calculational artifact, which presumably arises from the truncation error due to the insufficient number of our basis functions. In principle, these wiggles would vanish if we include an infinite number of orthogonal basis function. Typically in this work we choose about 20 harmonic oscillator basis functions. Perhaps we should seek a smarter set of basis functions that allows for a faster convergence behavior. On the other hand, we note that, whenever the wiggles appear, the corresponding φ − component is always typically about 3 or 4 orders-of-magnitude smaller than the φ + component, thus completely negligible in a practical sense.
With the BG wave functions φ ± (x, P ) available, we can employ the formulas derived in Section III B to compute the meson decay constant. In Fig. 12, we plot the decay constant of the ground state meson as functions of the quark mass and meson momentum. One finds an overall satisfactory agreement between the light-cone-gauge and axial-gauge predictions. In the Bars-Green formalism, we have explicitly examined that the meson decay constant is indeed frame-independent, as it must be. This can be viewed as another nontrivial verification of the Poincaré invariance of the Bars-Green formalism.

VI. SUMMARY
The 't Hooft model has constantly served a fruitful theoretical laboratory to sharpen our understanding about certain aspects of the realistic QCD. In contrast with the widelystudied light-front quantization of QCD 2 , much less work has been conducted in the equaltime quantization. The most notable formalism in this category is based on the axial gauge quantization, with the corresponding bound-state equations first developed by Bars and Green in late 1970s [17]. It was formally proved that when the meson is boosted to the IMF, the large component of the BG wave equation would exactly reduce to the 't Hooft wave function. Moreover, it is believed that Poincaré invariance should be preserved for the color-singlet meson wave function with arbitrary finite meson momentum. Unfortunately, until now this important feature has never been explicitly verified in a numerical fashion.
To date, the most comprehensive numerical solutions of the BG equation were those works done by Li and companions more than three decades ago, yet only for the stationary mesons [16,19]. In this paper, we have moved an important step forward, by numerically solving the Bars-Green equation for arbitrarily moving mesons, with meson species ranging from the chiral pion to heavy quarkonium. We are able to numerically establish the validity of Poincaré invariance of the 't Hooft model. Moreover, we have explicitly confirmed the tendency that, as the meson gets more and more boosted, the large component of the Bars-Green wave function is indeed approaching the corresponding LF wave function obtained in the light-cone gauge. We also computed the quark condensates and meson decay constants with a variety of meson momentum, and explicitly verified the frame-independence and gauge invariance of these physical quantities. As a topical application, the t' Hooft model may serve as a concrete toy model to extract some general features of the recently proposed quasi parton distributions [20]. We note that the relation between the 't Hooft light-cone gauge formulation and the Bar-Green axialgauge formulation for the two-dimensional QCD, is very similar to that between the LF parton distributions and the quasi parton distributions. Right now, the lattice simulation of the quasi distributions in the QCD 4 is still in its infancy. Therefore, we hope that our comprehensive understanding of the Bars-Green wave functions may shed some important light on the nature of quasi-distributions in realistic QCD [34].