text:
# Relativistic Many-Body Hamiltonian Approach to Mesons

## I Introduction

In a series of publications  an ambitious QCD program has been initiated to comprehensively investigate hadron structure. The theoretical formulation entails renormalization and utilizes established many-body techniques to approximately diagonalize an effective confining Hamiltonian. This paper, a detailed continuation of our recent letter , focuses upon the quark sector and reports numerical results for mesons complementing our previous gluon study .

Over the years there have been many meson investigations, from the early, simple non-relativistic constituent quark model calculations to more involved relativistic, field theoretical approaches implementing current quarks and spontaneous chiral symmetry breaking. A common shortcoming of these analyses is an inability to consistently reproduce the physical mass spectrum of the scalar and pseudoscalar mesons. Our paper addresses this issue and significantly extends the pioneering work of the Orsay group , Adler and Davis , and the Lisbon investigators . In our approach the exact QCD Hamiltonian in the Coulomb gauge is modeled by an effective, confining Hamiltonian, $`H`$, that is fully relativistic with quark field operators and current quark masses. However, before approximately diagonalizing $`H`$, a similarity transformation is implemented to a new quasiparticle basis having a dressed, but unknown constituent mass. As described in Sec. II, this transformation entails a rotation which mixes the bare quark creation and annihilation operators. By then performing a variational calculation to minimize the ground state (vacuum) energy, a specific angle and corresponding quasiparticle mass is selected. In this fashion chiral symmetry is dynamically broken and a non-trivial vacuum with quark condensates emerges. This treatment is precisely analogous to the Bardeen, Cooper, and Schrieffer (BCS) description of a superconducting metal as a coherent vacuum state of interacting quasiparticles combining to form condensates (Cooper pairs). Excited states (mesons) can then be represented as quasiparticle excitations using standard many-body techniques which in this work will be the Tamm-Dancoff (TDA) and random phase approximation (RPA) methods. The two treatments are truncated at the one quasiparticle, one quasihole level and then numerically compared. Our RPA analysis confirms and extends the early work of Ref.  which utilized an extended Nambu-Jona-Lasinio mean field approach.

Two other comments are in order before proceeding. First, there are several reasons for choosing the Coulomb gauge framework. As discussed by Zwanziger , the Hamiltonian is renormalizable in this gauge and, equally as important, the Gribov problem ($`𝐀=0`$ does not uniquely specify the gauge) can be resolved (see Refs.  for further discussion). Related, there are no spurious gluon degrees of freedom since only transverse gluons enter. This ensures all Hilbert vectors have positive normalizations which is essential for using variational techniques that have been widely successful in atomic, molecular and condensed matter physics. Second, due to Fock space truncations our analysis is not Lorentz invariant. However, we only plan to use one frame and do not compute hadron form factors with this method. Interestingly, violating Lorentz non-invariance implies a prefered reference frame, which, as selected by chiral symmetry breaking, is the condensate rest frame.

This paper consists of six sections and two appendices. The next section introduces our effective, QCD inspired Hamiltionian and developes the BCS vacuum treatment leading to the quasiparticle mass gap equation. We also compare our approach to the classic Nambu-Jona-Lasinio model. In Sec. III we detail our numerical, supercomputer solution of the gap equation along with the quark condensate and constituent mass values. Sections IV A and IV B describe the TDA and RPA, respectively, while Sec. IV C addresses weak decays and Sec IV D presents a derivation of the Gell-Mann-Oakes-Renner relation. The TDA and RPA meson spectra are compared and discussed in Sec. V. This section also includes results from a simple $`SU_f(3)`$ flavor mixing analysis for the $`\eta `$-$`\eta ^{}`$ system and our predictions for the charmed mesons. Conclusions and future work are summarized in Sec. VI. Finally, Appendix A provides further details regarding the BCS transformation and vacuum state while Appendix B presents the most general TDA equation for arbitrary angular momentum.

## II Hamiltonian and mass gap equation

### A Effective Hamiltonian

By introducing a phenomenological confining potential, $`V_L`$, the QCD Coulomb gauge Hamiltonian  for the quark sector can be replaced by an effective Hamiltonian

$$H=𝑑\stackrel{}{x}\mathrm{\Psi }^{}(\stackrel{}{x})(i\stackrel{}{\alpha }\stackrel{}{}+\beta m)\mathrm{\Psi }(\stackrel{}{x})\frac{1}{2}𝑑\stackrel{}{x}𝑑\stackrel{}{y}\rho ^a(\stackrel{}{x})V_L(|\stackrel{}{x}\stackrel{}{y}|)\rho ^a(\stackrel{}{y})$$
(1)
where $`\mathrm{\Psi }`$, $`m`$ and $`\rho ^a(\stackrel{}{x})=\mathrm{\Psi }^{}(\stackrel{}{x})T^a\mathrm{\Psi }(\stackrel{}{x})`$ are the current (bare) quark field, mass and color density, respectively (for a more complete discussion, especially for the heavy quark sector, consult Refs. ). For notational ease the flavor subscript is omitted (same $`H`$ for each flavor) and the color index runs $`a=1\mathrm{}8`$. Motivated by lattice gauge studies we adopt a linear confining interaction, $`V_L=\sigma |\stackrel{}{x}\stackrel{}{y}|`$, with slope $`\sigma =.18`$ $`GeV^2`$ also specified by lattice and Regge phenomenology. In our analysis we have also performed calculations with and without the leading QCD canonical or Coulomb (one-gluon exchange) interaction, $`V_C=\frac{\alpha _s}{|\stackrel{}{x}\stackrel{}{y}|}`$, with $`\alpha _s=\frac{g^2}{4\pi }`$ $``$ $`.4`$. For most observables, especially the meson mass spectrum, the Coulomb interaction is not important and can be omitted. This can be understood by noting that in momentum space, where we perform all calculations, the two interactions have the same sign, i.e. $`V(r=|\stackrel{}{x}\stackrel{}{y}|)=V_C+V_L`$

$`\widehat{V}(k)`$ $`=`$ $`{\displaystyle 𝑑\stackrel{}{r}V(r)e^{i\stackrel{}{k}\stackrel{}{r}}}`$ (2)
$`=`$ $`4\pi {\displaystyle \frac{\alpha _s}{k^2}}8\pi {\displaystyle \frac{\sigma }{k^4}}.`$ (3)
Because the meson wavefunctions have a finite momentum distribution, most static meson properties are predominantly governed by the infrared (IR), or low, momentum region where the confining potential dominates. Including the Coulomb interaction is then roughly equivalent to using a slightly larger string tension, $`\sigma `$. There are certain observables, and in particular the gap equation detailed below, for which the Coulomb interaction is ultra-violet (UV) divergent. In such cases we regularize with a cut-off parameter and then could renormalize to remove cut-off sensitivity using one of our renormalization procedures detailed in Refs.  for the gluon sector. In this paper we only present unrenormalized results since this program is still in progress  and has not yet been completed for the quark sector. This is an additional reason for omitting the Coulomb interaction. Hence, with the exception of the current quark masses (we use $`m_u=m_d=5`$ $`MeV`$, $`m_s=150`$ $`MeV`$, $`m_c=1200`$ $`MeV`$), our approach entails only one pre-determined parameter which also sets the hadronic scale, $`\sqrt{\sigma }=424`$ $`MeV`$.

We further note that even though the confining potential is IR divergent, this singularity is cancelled (see Ref. ) in both the mass gap equation and all calculations for associated observables. Hence, the problem is the delicate numerical evaluation of this integrable singularity which we discuss in Sec. III.

Finally, we stress that in constituent quark models free quarks can exist which requires imposing color confinement. However, as demonstrated in Refs.  the Lorentz structure of our Coulomb gauge density-density confining interaction only permits stable solutions for color singlet states. Therefore, confinement naturally emerges in our approach.

### B BCS transformation and gap equation

We now wish to solve $`H\mathrm{\Psi }=E\mathrm{\Psi }`$ as accurately as possible. In this subsection we focus on the ground state and introduce the Bogoliubov-Valatin, or BCS, transformation. We begin by recalling the plane wave, spinor expansion for the quark field operator

$$\mathrm{\Psi }(\stackrel{}{x})=\underset{c\lambda }{}\frac{d\stackrel{}{k}}{(2\pi )^3}\left[u_{c\lambda }(\stackrel{}{k})b_{c\lambda }(\stackrel{}{k})+v_{c\lambda }(\stackrel{}{k})d_{c\lambda }^{}(\stackrel{}{k})\right]e^{i\stackrel{}{k}\stackrel{}{x}}$$
(4)
with free particle, anti-particle spinors $`u_{c\lambda },v_{c\lambda }`$ and bare creation, annihilation operators $`b_{c\lambda },d_{c\lambda }`$ for current quarks, respectively. Here the spin state (helicity) is denoted by $`\lambda `$ and color index by $`c=1,2,3`$ (which is hereafter suppressed). Because we can expand $`\mathrm{\Psi }`$ in terms of any complete basis we may equally well use a new quasiparticle basis

$$\mathrm{\Psi }(\stackrel{}{x})=\underset{\lambda }{}\frac{d\stackrel{}{k}}{(2\pi )^3}\left[U_\lambda (\stackrel{}{k})B_\lambda (\stackrel{}{k})+V_\lambda (\stackrel{}{k})D_\lambda ^{}(\stackrel{}{k})\right]e^{i\stackrel{}{k}\stackrel{}{x}}$$
(5)
entailing quasiparticle spinors $`U_\lambda ,V_\lambda `$ and operators $`B_\lambda ,D_\lambda `$. The Hamiltonian is equivalent in either basis and the two are related by a similarity (Bogoliubov-Valatin or BCS) transformation. The transformation between operators is given by the rotation

$$B_\lambda (\stackrel{}{k})=\mathrm{cos}\frac{\theta _k}{2}b_\lambda (\stackrel{}{k})\lambda \mathrm{sin}\frac{\theta _k}{2}d_\lambda ^{}(\stackrel{}{k})$$
(6)
$$D_\lambda (\stackrel{}{k})=\mathrm{cos}\frac{\theta _k}{2}d_\lambda (\stackrel{}{k})+\lambda \mathrm{sin}\frac{\theta _k}{2}b_\lambda ^{}(\stackrel{}{k})$$

involving the BCS angle $`\theta _k=\theta (k)`$. Similarly the rotated quasiparticle spinors are

$$U_\lambda (\stackrel{}{k})=\mathrm{cos}\frac{\theta _k}{2}u_\lambda (\stackrel{}{k})\lambda \mathrm{sin}\frac{\theta _k}{2}v_\lambda (\stackrel{}{k})=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\sqrt{1+\mathrm{sin}\varphi (k)}\chi _\lambda \\ \sqrt{1\mathrm{sin}\varphi (k)}\stackrel{}{\sigma }\widehat{k}\chi _\lambda \end{array}\right]$$
(7)
$$V_\lambda (\stackrel{}{k})=\mathrm{cos}\frac{\theta _k}{2}v_\lambda (\stackrel{}{k})+\lambda \mathrm{sin}\frac{\theta _k}{2}u_\lambda (\stackrel{}{k})=\frac{1}{\sqrt{2}}\left[\begin{array}{c}\sqrt{1\mathrm{sin}\varphi (k)}\stackrel{}{\sigma }\widehat{k}\chi _\lambda \\ \sqrt{1+\mathrm{sin}\varphi (k)}\chi _\lambda \end{array}\right]$$

where $`\chi _\lambda `$ is the standard two-dimensional Pauli spinor. We have also introduced the gap angle, $`\varphi _k=\varphi (k)`$, which is related to the BCS angle, $`\theta /2`$, by $`\varphi =\theta +\alpha `$ where $`\alpha `$ is the current, or perturbative, mass angle satisfying $`\mathrm{sin}\alpha =m/E_k`$ with $`E_k=\sqrt{m^2+k^2}`$. Hence

$$\mathrm{sin}\varphi _k=\frac{m}{E_k}\mathrm{cos}\theta _k+\frac{k}{E_k}\mathrm{sin}\theta _k$$

$$\mathrm{cos}\varphi _k=\frac{k}{E_k}\mathrm{cos}\theta _k\frac{m}{E_k}\mathrm{sin}\theta _k.$$

Similarly, the perturbative, trivial vacuum, defined by $`b_\lambda |0=d_\lambda |0=0`$, is related to the quasiparticle vacuum, $`B_\lambda |\mathrm{\Omega }=D_\lambda |\mathrm{\Omega }=0`$, by the transformation

$$|\mathrm{\Omega }=\mathrm{exp}\left(\underset{\lambda }{}\frac{d\stackrel{}{k}}{(2\pi )^3}\lambda \mathrm{tan}\frac{\theta _k}{2}b_\lambda ^{}(\stackrel{}{k})d_\lambda ^{}(\stackrel{}{k})\right)|0.$$
(8)
In this paper we will denote the BCS vacuum by $`|\mathrm{\Omega }`$ (in Sec. IV we introduce the RPA vacuum labeled $`|\mathrm{\Omega }_{RPA}`$). Expanding the exponential and noting that the form of the operator $`b^{}d^{}`$ is designed to create a current quark/antiquark pair with the vacuum quantum numbers, clearly exhibits the BCS vacuum as a coherent state of quark/antiquark excitations (Cooper pairs) representing $`{}_{}{}^{2S+1}L_{J}^{}=`$3P0 condensates. One can regard $`\mathrm{tan}\frac{\theta _k}{2}`$ as the momentum wavefunction of the pair in the center of momentum system.

We now seek an approximate ground state for our effective Hamiltonian by minimizing the BCS vacuum expectation, $`\mathrm{\Omega }|H|\mathrm{\Omega }`$. We do this variationally using the gap angle, $`\varphi _k`$, (not the BCS angle) which leads to the gap equation, $`\delta \mathrm{\Omega }|H|\mathrm{\Omega }=0`$. After considerable mathematical reduction, the nonlinear integral gap equation follows

$$k\mathrm{sin}\varphi _km\mathrm{cos}\varphi _k=\frac{2}{3}\frac{d\stackrel{}{q}}{(2\pi )^3}\widehat{V}(|\stackrel{}{k}\stackrel{}{q}|)[\mathrm{sin}\varphi _k\mathrm{cos}\varphi _q\widehat{k}\widehat{q}\mathrm{sin}\varphi _q\mathrm{cos}\varphi _k].$$
(9)
The angular integrals can be analytically evaluated (see Appendix B) to give

$$k\mathrm{sin}\varphi _km\mathrm{cos}\varphi _k=\frac{2}{3}\frac{1}{(2\pi )^2}_0^{\mathrm{}}q^2𝑑q[\mathrm{cos}\varphi _q\mathrm{sin}\varphi _k\widehat{V}_1(k,q)\mathrm{cos}\varphi _k\mathrm{sin}\varphi _q\widehat{V}_0(k,q)]$$
(10)
where

$$\widehat{V}_0=\frac{16\pi \sigma }{(k^2q^2)^2}$$

$$\widehat{V}_1=\frac{2\pi \sigma }{k^2q^2}\left[ln\left(\frac{k+q}{kq}\right)^2+(k^2+q^2)\left(\frac{4qk}{(k^2q^2)^2}\right)\right]$$

corresponding to the linear potential above. Similar expressions for the Coulomb potential are given in Appendix B.

There are several alternative ways to derive this same gap equation. One is through the Ward identites. Another is by requiring cancellation of the anomalous Bogoliubov terms in the 2-body part of the newly normal ordered Hamiltonian. The latter is necessary to stabilize the vacuum and is also equivalent to minimizing the 0-body constant energy splitting the BCS and trivial vacua (see Ref. ).

Numerically we actually solve a different form of the gap equation, originally obtained by Adler and Davis , that is more familiar to the solid state community. They use the function $`\psi _k=\psi (k)`$ related to our gap angle by

$$\mathrm{sin}\varphi _k=\frac{2\psi _k}{1+\psi _k^2}$$

$$\mathrm{cos}\varphi _k=\frac{1\psi _k^2}{1+\psi _k^2}$$

with corresponding gap equation

$$k\psi _k\frac{m}{2}(1\psi _k^2)=\frac{2}{3(2\pi )^2}_0^{\mathrm{}}q^2𝑑q\frac{\widehat{V}_1\psi _k(1\psi _q^2)\widehat{V}_0\psi _q(1\psi _k^2)}{1+\psi _q^2}.$$
(11)
Examination of Eqs. (10,11) reveals that the divergence at $`k=q`$ is an integrable singularity for the linear potential since for all $`k`$ the integrands vanish at $`k=q`$. This is not the case for the Coulomb potential UV singularity. However, since it naturally emerges from the canonical QCD Hamiltonian (one gluon exchange), we retain the option of including this potential for selected calculations and use a cut-off to regulate its ultraviolet (UV) divergence.

The solution of the gap equation (see Sec. III) leads to a vacuum quark-antiquark condensate given by

$$\overline{q}q\mathrm{\Omega }|\overline{\mathrm{\Psi }}(0)\mathrm{\Psi }(0)|\mathrm{\Omega }=\frac{3}{\pi ^2}k^2\mathrm{sin}\varphi _kdk$$
(12)
which is quadratically divergent for non-zero current quark mass $`m0`$. We regulate this by subtracting the trivial condensate contribution giving

$$\overline{q}q_{reg}=\frac{3}{\pi ^2}k^2\left(\mathrm{sin}\varphi _k\frac{m}{E_k}\right)𝑑k.$$
(13)
Our model is color confining and does not permit free solitary particles since the self-energy or dispersion relation

$$ϵ_k=m\mathrm{sin}\varphi _k+k\mathrm{cos}\varphi _k\frac{2}{3}\frac{d\stackrel{}{q}}{(2\pi )^3}\widehat{V}(|\stackrel{}{k}\stackrel{}{q}|)(\mathrm{sin}\varphi _k\mathrm{sin}\varphi _q+\widehat{k}\widehat{q}\mathrm{cos}\varphi _k\mathrm{cos}\varphi _q)$$
(14)
obtained from the 1-body part of $`H`$, Eq. (1), is divergent (now there is no cancellation at the singular point $`\stackrel{}{k}=\stackrel{}{q}`$). Further, this divergence is also cancelled in the bound state equation but only for color singlet states (see below). Even though the self-energy diverges it is still useful to introduce the concept of an effective quasiparticle (constituent) mass, $`\mathrm{m}_q`$, which can be extracted from the low momentum behavior of the gap angle. We introduce a running, dynamical mass, $`\mathrm{m}(k)`$, by an effective Dirac spinor in canonical form

$$U_\lambda ^{eff}(\stackrel{}{k})=N\left[\begin{array}{c}\chi _\lambda \\ \frac{\stackrel{}{\sigma }\stackrel{}{k}}{E+\mathrm{m}(k)}\chi _\lambda \end{array}\right]$$
(15)
with normalization $`N`$ and $`E=\sqrt{\mathrm{m}^2(k)+k^2}`$. Then using this equation and Eq. (7) we equate the two relative normalizations between upper and lower spinor components yielding a relation between the running dynamical mass and gap angle

$$\frac{\sqrt{1+\mathrm{sin}\varphi (k)}}{\sqrt{1\mathrm{sin}\varphi (k)}}=\frac{E+\mathrm{m}(k)}{k},$$
(16)
or

$$\mathrm{sin}\varphi (k)=\frac{\mathrm{m}(k)}{E}=1\frac{k^2}{2\mathrm{m}^2(k)}+\vartheta (k^4).$$
(17)
We identify the dressed quark or quasiparticle mass as $`\mathrm{m}_q=maxima\mathrm{m}(k)`$ and extract it from the low momentum behavior of the gap angle (essentially inverse of the slope near zero momentum). The value of $`\mathrm{m}_q`$ characterizes the degree of chiral symmetry breaking and can be loosely regarded as the constituent quark mass associated with phenomenological quark models.

Note that our expression for the running mass is functionally identical to the perturbative expression $`\mathrm{tan}\alpha =\frac{m}{k}`$. Related, since the rotated quasiparticle spinors have a running momentum dependence, they no longer rigorously provide a representation of the Lorentz group. Our form of dynamical chiral symmetry breaking violates Lorentz invariance which implies a preferred reference frame, namely the condensate rest frame. For most static observables such as masses, condensates and decay constants, Lorentz symmetry is not important. However, for some observables, such as electromagnetic form factors, care is necessary and boost corrections may be important. This issue is under investigation and will be reported in a future communication.

### C Comparison to the Nambu-Jona-Lasinio model

The classical effective model of Nambu-Jona-Lasinio (NJL) (see Ref.  for review) entails various Lagrangian formulations, a common one being

$$=i\overline{\mathrm{\Psi }}\overline{)}\mathrm{\Psi }+G[(\overline{\mathrm{\Psi }}\mathrm{\Psi })(\overline{\mathrm{\Psi }}\mathrm{\Psi })(\overline{\mathrm{\Psi }}\gamma _5\mathrm{\Psi })(\overline{\mathrm{\Psi }}\gamma _5\mathrm{\Psi })]$$
(18)
where $`G`$ is a constant. It is customary to introduce the approximations $`\overline{\mathrm{\Psi }}\gamma _5\mathrm{\Psi }0`$ and $`\overline{\mathrm{\Psi }}\mathrm{\Psi }\overline{\mathrm{\Psi }}\mathrm{\Psi }_{vacuum}`$ to linearize the equations of motion and then extract a constituent quark mass from the NJL mass gap equation. In this fashion chiral symmetry breaking is achieved.

Our formulation extends beyond the NJL model in several important ways.

1. Our approach is more general and permits explicit gluonic degrees of freedom (see Refs. ). The unification of the quark and glue sectors is crucial for a comprehensive treatment of hadron structure, especially for glueball and hybrid meson systems.
2. Our formulation includes confinement and is renormalizable while the NJL model has neither. The NJL pointlike interaction would be recovered in the limit $`V(\stackrel{}{x}\stackrel{}{y})G\delta (\stackrel{}{x}\stackrel{}{y})`$ which removes all important nonlocalities.
3. Our model has a density-density interaction kernel with a different Lorentz structure, $`\gamma _0\gamma _0`$, which is the product of four-vector time components. As discussed in Ref. , a density-density (vector-vector) interaction is superior to the scalar-pseudoscalar displayed by the NJL model.
4. The chiral symmetry breaking mode of the NJL is extremely restrictive yielding a constant quasiparticle mass, $`m_{dyn}`$, and simple dispersion $`E=\sqrt{m_{dyn}^2+k^2}`$. Related, the NJL limit of our model also yields a more restricted gap angle since $`\mathrm{sin}\varphi =\frac{m_{dyn}}{\sqrt{m_{dyn}^2+k^2}}`$. Our method has a running mass and different quasiparticle dispersion which yields more realistic TDA and RPA hadron masses.

## III Numerical solution of the gap equation

The gap equation (10) has been previously solved for the harmonic oscillator potential, where it takes its simplest form as a differential equation (see Ref. ), and also for the linear potential (see Ref. ). Here we summarize our analysis which confirms and extends the latter results.

To numerically treat the integrable IR singularity a regularization must still be implemented even though the final results are independent of this procedure. We considered several different regularizations. We first tried an analytical regulator (equivalent to a deconfining correction to the potential). This was an unstable algorithm and and convergence could not be achieved. We next examined the method of Ref.  which off-sets the q-discretization by half a step in the kernel with respect to the k-discretization. This procedure was also rejected as it was less amenable for documenting the regulator sensitivity. We finally adopted the simplest method of omitting the singular point $`k=q`$. This also facilitated a controlled senstivity study by just increasing the number of mesh points. Related, we adopted a variable mesh size to integrate more efficiently and mapped the integration variable $`q`$ to $`v`$

$$q=\frac{v^2q_{max}}{1+q_{max}(1v)}$$

for N points uniformly distributed in the interval $`v(0,1)`$.

Following Ref.  we elected to solve the gap equation in form specified by Eq. (11) and also utilized the Gauss algorithm as described there. The Gauss method assures convergence but is rather inefficient for extensive sensivity studies in parameter space. We therefore modified our numerical approach by first finding a good approximate solution, $`\psi _k^0`$, to the non-linear gap equation and then obtained a linear equation for the desired correction, $`\delta _k`$, giving the final solution

$$\psi _k=\psi _k^0+\delta _k$$
(19)
to arbitrary accuracy. Substituting Eq. (19) in the gap equation, Eq. (11), and dropping higher powers of $`\delta _k`$ yields the approximate linear equation

$`\delta _k\left[(k+m)\psi _k^0{\displaystyle \frac{2}{3(2\pi )^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{q^2dq}{1+\psi _q^{0^2}}}\left[\widehat{V}_1(1\psi _q^{0^2})+2\widehat{V}_0\psi _q^0\psi _k^0\right]\right]`$
$`+{\displaystyle \frac{2}{3(2\pi )^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{q^2dq}{1+\psi _q^{0^2}}}\left[\widehat{V}_0\left(1{\displaystyle \frac{2\psi _q^{0^2}}{1+\psi _q^{0^2}}}\right)(1\psi _k^{0^2})+4\widehat{V}_1{\displaystyle \frac{\psi _k^0\psi _q^0}{1+\psi _q^{0^2}}}\right]\delta _q=`$
$`=k\psi _k^0+{\displaystyle \frac{m}{2}}(1\psi _k^{0^2}){\displaystyle \frac{2}{3(2\pi )^2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{q^2dq}{1+\psi _q^{0^2}}}\left[\widehat{V}_0\psi _q^0(1\psi _k^{0^2})\widehat{V}_1\psi _k^0(1\psi _q^{0^2})\right].`$
This equation is of the form

$$𝑑q𝒜(k,q)\delta (q)=(k)$$
(20)
which can be solved for $`\delta _k`$ by matrix inversion. We found the Gauss algorithm sufficient for obtaining the initial approximate solution $`\psi _k^0`$. To achieve full convergence required up to 12,000 mesh points, a factor of 60 more than the early calculations of Ref. .

We checked our computer codes by calculating two different toy kernels $`\widehat{V}_0=\frac{k}{q^2(1+q^2)}`$, $`\widehat{V}_1=0`$ and $`\widehat{V}_0=0`$, $`\widehat{V}_1=\frac{k}{q^2(1+q^2)}`$, each designed to yield a known constant value for $`\psi _k`$. We then performed a series of cut-off sensitivity runs and mapped out the convergence rate as a function of mesh point number $`N`$ which ranged from 100 to 12,000. We used the quark condensate as a test observable and also performed calculations for zero and non-zero current quark mass, $`m`$, with and without the Coulomb potential using $`\alpha _s=.4`$. For $`m=0`$, $`\alpha _s=0`$, we determined the sensitivity to $`N`$ (the effective cut-off parameter) was slightly higher than previously reported  and given by

$$\overline{\mathrm{\Psi }}\mathrm{\Psi }\left[\left(113\frac{1400}{N}\right)MeV\right]^3.$$
(21)
Note that this number is somewhat smaller than the commonly accepted lattice value of about $`(250MeV)^3`$. Including the Coulomb potential only increases the condensate to $`(119MeV)^3`$. We therefore conclude an improved model ground state is needed which can be provided by including additional terms in the Hamiltonian, such as the quark-gluon minimal coupling (hyperfine) interaction. This point is also affirmed below in our RPA treatment which does yield a more realistic condensate value.

Our other key result, which will be of interest in connection with chiral perturbation theory , is for the constituent quark mass and the BCS condensate as a function of the $`u`$, $`d`$ quark mass. Now it is necessary to use Eq. (13) and also impose an additional integration cut-off limit ($`q_{max}`$ around 10 $`GeV`$). This yields

$$\overline{\mathrm{\Psi }}\mathrm{\Psi }^{\frac{1}{3}}=2.03m+113.1$$

$$\mathrm{m}_q=1.6m+77.9$$

where the units are $`MeV`$. A precision calculation for $`\mathrm{m}_q`$ with $`m=0`$ yields the slightly higher value $`\mathrm{m}_q=`$ $`80.5MeV`$ (all values for the linear potential only).

Finally, we note that our Hamiltonian is $`SU_f(3)`$ flavor symmetric, broken only by the small current flavored quark mass term. However, and quite significant, the vacuum properties and gap angle exhibit substantial $`SU_f(3)`$ violations as evidenced by our strange quark calculations using a current mass of 150 $`MeV`$. Important violations occur even for a strange quark mass as low as 50 $`MeV`$. While this result may be model dependent it does suggest that certain chiral symmetry arguments in the literature regarding the strange quark sector should be taken with care.

## IV Many-Body Techniques

We now formulate mesons as excited states consisting of quasiparticles and seek approximate eigensolutions of our Hamiltonian. We first develop the TDA and then treat the RPA in the next subsection. Of the two, only the RPA preserves chiral symmetry, as we detail below. It is therefore more closely related to the Bethe-Salpeter formalism  incorporating the Schwinger-Dyson quark propagator using an instantaneous interaction in the rainbow approximation (equivalent to our gap equation). We will document this connection more formally in a future publication.

### A TDA equation of motion

The principle advantage of the TDA is that it is a controllable approximation which truncates the Fock-space expansion for a chosen level of calculational effort and resources. In terms of the quasiparticle operators introduced in Sec. II, we introduce the TDA meson creation operator

$$Q_{nJP}^{}(TDA)=\underset{\gamma \delta }{}\frac{d\stackrel{}{k}}{(2\pi )^3}\mathrm{\Psi }_{\gamma \delta }^{nJP}(\stackrel{}{k})B_\gamma ^{}(\stackrel{}{k})D_\delta ^{}(\stackrel{}{k}).$$
(22)
A meson with quantum numbers $`nJP`$ (radial-node number, $`n`$, total angular momentum, $`J`$, and parity, $`P`$) is then represented by the Fock space expansion

$$|\mathrm{\Psi }_{TDA}^{nJP}=Q_{nJP}^{}(TDA)|\mathrm{\Omega }$$
(23)
containing a quasiparticle and quasihole excited from the BCS vacuum. The Hamiltonian equation is then projected onto this 1p-1h truncated Fock sector giving the TDA equation

$$\mathrm{\Psi }_{TDA}^{nJP}|[H,B_\alpha ^{}D_\beta ^{}]|\mathrm{\Omega }=(E_{nJP}E_0)\mathrm{\Psi }_{\alpha \beta }^{nJP}.$$
(24)
In evaluating the commutator we note

$$[H_0,B_\alpha ^{}D_\beta ^{}]=0$$

$$\mathrm{\Psi }_{TDA}^{nJP}|[H_2,B_\alpha ^{}D_\beta ^{}]|\mathrm{\Omega }=(ϵ_k+ϵ_k)\mathrm{\Psi }_{\alpha \beta }^{nJP}$$

and for the two body potential

$$\mathrm{\Psi }_{TDA}^{nJP}|[H_4,B_\alpha ^{}D_\beta ^{}]|\mathrm{\Omega }=\frac{4}{3}\underset{\gamma \delta }{}\mathrm{\Psi }_{\delta \gamma }^{nJP}\widehat{V}(\stackrel{}{k}_\delta \stackrel{}{k}_\gamma )U_\alpha ^{}U_\delta V_\gamma ^{}V_\beta $$

where $`ϵ_k`$ is the BCS gap energy given by (14) and $`H_N`$ is the Hamiltonian component containing N field operators (after normal ordering with respect to the BCS vacuum).

We can exploit the rotational invariance of our Hamiltonian and reduce the linear TDA equation to a one dimensional, nonlocal equation by an angular momentum decomposition. Introducing the orbital and spin angular momenta $`\stackrel{}{L}`$ and $`\stackrel{}{S}`$, respectively, the meson state vector can be expanded in partial-waves involving a one-dimensional (radial) wavefunction $`\mathrm{\Psi }_{LS}^{nJP}`$

$$\mathrm{\Psi }_{\delta \gamma }^{nJP}(\stackrel{}{k})=\underset{LSm_Lm_S}{}Lm_LSm_S|Jm_J(1)^{\frac{1}{2}+\gamma }\frac{1}{2}\delta \frac{1}{2}\gamma |Sm_SY_L^{m_L}(\widehat{k})\mathrm{\Psi }_{LS}^{nJP}(k)$$
(25)
where again the color index is omitted. Note the phase factor and negative magnetic substate sign in the Clebsch-Gordan coefficient due to the transformation properties of antiparticles under the $`SU(2)`$ rotation group. A thorough discussion is given in Ref. .

Inserting Eq. (25) into Eq. (24) yields the TDA partial-wave equation of motion

$`(E_{nJP}E_02ϵ_k)\mathrm{\Psi }_{LS}^{nJP}(k)(2J+1)={\displaystyle \underset{\mathrm{\Lambda }\mathrm{\Sigma }m_\mathrm{\Lambda }m_\mathrm{\Sigma }m_Jm_Lm_S}{}}Jm_J|Lm_LSm_S\mathrm{\Lambda }m_\mathrm{\Lambda }\mathrm{\Sigma }m_\mathrm{\Sigma }|Jm_J`$
$`{\displaystyle 𝑑\mathrm{\Omega }_k𝑑\mathrm{\Omega }_qY_L^{m_L}(\widehat{k})Y_\mathrm{\Lambda }^{m_\mathrm{\Lambda }}(\widehat{q})\frac{4}{3}_0^{\mathrm{}}\frac{q^2dq}{(2\pi )^3}\mathrm{\Psi }_{\mathrm{\Lambda }\mathrm{\Sigma }}^{nJP}(q)\widehat{V}(|\stackrel{}{k}\stackrel{}{q}|)}`$
$`{\displaystyle \underset{\alpha \beta \gamma \delta }{}}(1)^{1+\beta +\gamma }h_{\gamma \delta }^{\alpha \beta }(k,q)Sm_S|{\displaystyle \frac{1}{2}}\alpha {\displaystyle \frac{1}{2}}\beta {\displaystyle \frac{1}{2}}\delta {\displaystyle \frac{1}{2}}\gamma |\mathrm{\Sigma }m_\mathrm{\Sigma }`$
where the function $`h_{\gamma \delta }^{\alpha \beta }(k,q)`$ contains the gap angle from contractions involving rotated spinors

$`h_{\gamma \delta }^{\alpha \beta }(k,q)={\displaystyle \frac{1}{4}}[c_kc_q(\delta _{\alpha \delta }g_{\gamma \beta }+\delta _{\gamma \beta }g_{\alpha \delta })+`$ $`(1+s_k)(1+s_q)\delta _{\alpha \delta }\delta _{\beta \gamma }+(1s_q)(1s_k)g_{\alpha \delta }g_{\gamma \beta }]`$
with

$$g_{\alpha \beta }=\chi _\alpha ^{}\stackrel{}{\sigma }\widehat{q}\stackrel{}{\sigma }\widehat{k}\chi _\beta .$$

Denoting the meson mass for state $`nJP`$ by $`M_{nJP}=E_{nJP}E_0`$ and using the multipole expansion formulas for the interaction yields the final TDA equation appropriate for numerical calculation

$$(M_{nJP}2ϵ_k)\mathrm{\Psi }_{LS}^{nJP}(k)=\underset{\mathrm{\Lambda }\mathrm{\Sigma }}{}_0^{\mathrm{}}K_{L\mathrm{\Lambda }S\mathrm{\Sigma }}^{JP}(k,q)\mathrm{\Psi }_{\mathrm{\Lambda }\mathrm{\Sigma }}^{nJP}(q)\frac{q^2dq}{12\pi ^2}.$$
(26)
Note the Hamiltonian spin dependence generates a kernel that couples different orbital and spin states.

We now apply these equations to the low lying meson spectrum with quantum states specified by $`I^G(J^{PC})`$ having $`C`$ parity, $`C=(1)^{L+S}`$, and $`G`$ parity, $`G=(1)^{L+S+I}`$. In our model we neglect the small electromagnetic (isospin violating) effects as well as coupling to the gluon sector so that $`I=0`$ and $`1`$ states are degenerate for the same $`J^{PC}`$. For pseudoscalar states, $`J^{PC}=0^+`$, $`S=L=J=0`$ giving only one wavefunction component (no coupling). This is also the case for scalar mesons ($`L=1,S=1`$) having $`J^{PC}=0^{++}`$. However for the vector meson sector $`J^{PC}=1^{}`$ both $`(L=0,S=1)`$ and $`(L=2,S=1)`$ waves are allowed and, in general, will be coupled. Similarly for low lying pseudovector mesons having $`(L=1,S=0`$ or $`S=1)`$ and tensor mesons, with $`(L=1,S=2)`$ and $`(L=3,S=2)`$, there will be coupled equations. Although these equations are not difficult to solve, we only include the lowest orbital partial-wave component and neglect all coupling since it has been computed small for the harmonic oscillator potential . There is then only one kernel $`K(k,q)`$ for each meson state with structure:

* pseudoscalar, $`L=S=J=0`$,

$$K(k,q)=2(c_kc_q\widehat{V}_1+(1+s_ks_q)\widehat{V}_0)$$
* scalar, $`L=S=1`$, $`J=0`$,

$$K(k,q)=2(c_kc_q\widehat{V}_0+(1+s_ks_q)\widehat{V}_1)$$
* vector, $`L=0`$, $`S=J=1`$ (neglecting the tensor $`L=2,S=J=1`$ coupling),

$$K(k,q)=2c_kc_q\widehat{V}_1+(1+s_k)(1+s_q)\widehat{V}_0+(1s_q)(1s_k)(\frac{4\widehat{V}_2\widehat{V}_0}{3})$$
* pseudovector, $`L=J=1`$, $`S=0`$ (degenerate with $`L=S=J=1`$),

$$K(k,q)=c_kc_q(\widehat{V}_0+\widehat{V}_2)+2(1+s_ks_q)\widehat{V}_1$$
* tensor, $`L=S=1`$, $`J=2`$ (neglecting $`L=3,S=J=2`$ coupling),

$$K(k,q)=c_kc_q(3\widehat{V}_2\widehat{V}_0)+(1+s_k)(1+s_q)\widehat{V}_1+(1s_k)(1s_q)\frac{12\widehat{V}_37\widehat{V}_1}{5}.$$

Note for the pseudovector mesons the kernels for $`S=0`$ and $`S=1`$ are identical which differs from Ref. .

We have also applied our approach to other flavored ($`s`$ and $`c`$) meson systems and have obtained similar, but more complicated TDA equations. As a representative result, consider the pseudoscalar $`D`$ meson with a $`u`$ (or $`d`$) and $`c`$ quark. The gap equation for the $`c`$ quark remains the same, except for current mass, now $`1.2GeV`$, which gives a different gap energy, $`ϵ_k^c`$, and angle, $`s_k^c`$. The TDA equation, however, has a different form and generalizes to

$`(M_Dϵ_k^uϵ_k^c)\mathrm{\Psi }_D(k)={\displaystyle \frac{1}{3}}{\displaystyle }{\displaystyle \frac{q^2dq}{(2\pi )^2}}\mathrm{\Psi }_D(q)`$
$`[(\sqrt{1+s_k^u}\sqrt{1+s_q^u}\sqrt{1+s_k^c}\sqrt{1+s_q^c}+\sqrt{1s_q^u}\sqrt{1s_k^u}\sqrt{1s_k^c}\sqrt{1s_q^c})\widehat{V}_0(k,q)`$
$`(\sqrt{1s_q^u}\sqrt{1s_k^u}\sqrt{1+s_k^c}\sqrt{1+s_q^c}+\sqrt{1+s_k^u}\sqrt{1+s_q^u}\sqrt{1s_k^c}\sqrt{1s_q^c})\widehat{V}_1(k,q)]`$
with obvious form for other mixed flavors. All equations are finite for $`k=q`$ as the IR divergence terms from the confining potential again cancel. We have also derived and solved the TDA equations for other spin parity states which is further detailed in Sec. V.

### B RPA and the quasiboson approximation

The TDA can be improved by utilizing a better vacuum with additional quasiparticle correlations beyond the BCS. Consistent with many-body applications in other disciplines we now formulate the RPA  and introduce a new vacuum, $`|\mathrm{\Omega }_{RPA}`$, having both fermion (two quasiparticles or Cooper pairs) and boson (four quasiparticles or meson pairs) correlations. The RPA meson state

$$|\mathrm{\Psi }_{RPA}^{nJP}=Q_{nJP}^{}(RPA)|\mathrm{\Omega }_{RPA}$$
(27)
involves a meson creation operator which is a generalization of Eq. (22)

$$Q_{nJP}^{}(RPA)=\underset{\lambda \mu }{}\frac{d\stackrel{}{k}}{(2\pi )^3}[X_{\lambda \mu }^{nJP}B_\lambda ^{}(\stackrel{}{k})D_\mu ^{}(\stackrel{}{k})Y_{\lambda \mu }^{nJP}B_\lambda (\stackrel{}{k})D_\mu (\stackrel{}{k})].$$
(28)
The RPA vacuum then satisfies

$$Q_{nJP}(RPA)|\mathrm{\Omega }_{RPA}=0$$

which, because of additional correlations from admixtures of particle-hole excitation states, is not true for the BCS vacuum.

To derive the RPA equations of motion we use Eq. (24) and replace the BCS vacuum with $`|\mathrm{\Omega }_{RPA}`$ and also substitute $`\mathrm{\Psi }_{RPA}^{nJP}`$ for $`\mathrm{\Psi }_{TDA}^{nJP}`$ to generate one equation for the $`X`$ component. We then repeat, changing the commutator to $`[H,B_\alpha D_\beta ]`$ to obtain the $`Y`$ equation. Following standard treatments in other fields of physics, we also invoke the quasiboson approximation and treat the fermion pair operator $`BD`$ as a pure boson operator. This significantly reduces the commutator algebra complexity and generates one of the two coupled equations for the RPA wavefunctions $`X`$ and $`Y`$.

For the important pseudoscalar meson channel we obtain for the excited state $`n`$

$`2ϵ_kX^n(k)+{\displaystyle \frac{1}{3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{q^2dq}{(2\pi )^2}}\left[X^n(q)F(k,q)+Y^n(q)G(k,q)\right]=M_nX^n(k)`$ (29)
$`2ϵ_kY^n(k)+{\displaystyle \frac{1}{3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{q^2dq}{(2\pi )^2}}\left[Y^n(q)F(k,q)+X^n(q)G(k,q)\right]=M_nY^n(k)`$ (30)
where

$`F(k,q)=2c_qc_k\widehat{V}_1+2(1+s_qs_k)\widehat{V}_0`$ (31)
$`G(k,q)=2c_qc_k\widehat{V}_12(1s_qs_k)\widehat{V}_0.`$ (32)
Similar expressions directly follow for the other spin-parity states.

We adopt the standard normalization for the RPA wavefunctions

$$\nu ^{}|\nu =\mathrm{\Omega }_{RPA}|Q_\nu ^{}Q_\nu ^{}|\mathrm{\Omega }_{RPA}=\delta _{\nu \nu ^{}}$$
(33)
yielding

$$_0^{\mathrm{}}k^2𝑑k(X^\nu ^{}(k)^{}X^\nu (k)Y^\nu ^{}(k)^{}Y^\nu (k))=(2\pi )^3\delta _{\nu \nu ^{}}.$$

The RPA equations, which reduce to the TDA equations in the limit $`Y`$ or $`G`$ $`0`$, are again an eigenvalue problem for $`M_{nJP}`$ which can be easily diagonalized. Related, the matrix size can be reduced by a factor of 2 using the variables $`X+Y`$ and $`XY`$. Finally, the equations are also IR finite for the sigular point $`k=q`$.

### C Weak decay constants

A crucial test of any approach is the ability to describe hadronic decays. In this paper we compute weak decays and defer our analysis of hadronic decays to a subsequent publication. For a pseudoscalar meson $`P`$ with momentum $`p_\mu `$, energy $`E_P`$, mass $`M_P`$, the weak decay constant, $`f_P`$, is defined for our normalization by

$$\mathrm{\Omega }|A_\mu (0)|P(\stackrel{}{p})=\frac{1}{\sqrt{E_P}}f_Pp_\mu .$$
(34)
Here $`A_\mu (\stackrel{}{x})=\overline{\mathrm{\Psi }}(\stackrel{}{x})\gamma _\mu \gamma _5\mathrm{\Psi }(\stackrel{}{x})`$ is the axial current which specifies the chiral charge operator

$$Q_5=𝑑\stackrel{}{x}A_0(\stackrel{}{x})=𝑑\stackrel{}{x}\mathrm{\Psi }^{}(\stackrel{}{x})\gamma _5\mathrm{\Psi }(\stackrel{}{x}).$$
(35)
Simplifying Eq. (34) for a meson at rest yields

$$f_P=\frac{1}{\sqrt{M_P}}\mathrm{\Omega }|\mathrm{\Psi }^{}(0)\gamma _5\mathrm{\Psi }(0)|P(0).$$
(36)
Applying this result for the TDA pion wavefunction gives the TDA pion decay constant

$$f_\pi ^{TDA}=\frac{1}{\pi \sqrt{(2\pi )^3M_\pi }}_0^{\mathrm{}}\mathrm{\Psi }_{TDA}^\pi (q)s_qq^2𝑑q.$$
(37)
Similarly, for the RPA pion wavefunction we obtain

$$f_\pi ^{RPA}=\frac{1}{\pi \sqrt{(2\pi )^3M_\pi }}_0^{\mathrm{}}s_q(X^\pi (q)+Y^\pi (q))q^2𝑑q.$$
(38)
Our results easily generalize to the $`SU_f(3)`$ flavor nonet. Now there are nine axial charges given by

$$Q_5^a=𝑑\stackrel{}{x}A_0^a(\stackrel{}{x})=𝑑\stackrel{}{x}\mathrm{\Psi }^{}(\stackrel{}{x})\gamma _5\frac{\lambda ^a}{2}\mathrm{\Psi }(\stackrel{}{x})$$

where the eight Gell-Mann $`\lambda ^a`$ matrices are supplemented by $`\lambda ^0=\sqrt{\frac{2}{3}}I`$ to obtain both the octet and the singlet under $`SU_f(3)`$ transformations. The appropriate generalizations of Eq. (37) are then:

$$f_K^{TDA}=\frac{1}{2\pi \sqrt{(2\pi )^3M_K}}_0^{\mathrm{}}q^2𝑑q\mathrm{\Psi }_{TDA}^K(q)(\sqrt{1+s_q^s}\sqrt{1+s_q^u}\sqrt{1s_q^s}\sqrt{1s_q^u})$$

$$f_{\eta _8}^{TDA}=\frac{1}{3\pi \sqrt{(2\pi )^3M_{\eta _8}}}_0^{\mathrm{}}\mathrm{\Psi }_{TDA}^{\eta _8}(s_q^u+2s_q^s)q^2𝑑q$$

$$f_{\eta _0}^{TDA}=\frac{1}{3\pi \sqrt{(2\pi )^3M_{\eta _0}}}_0^{\mathrm{}}\mathrm{\Psi }_{TDA}^{\eta _0}(2s_q^u+s_q^s)q^2𝑑q.$$

We will use the above results in the next subsection to derive a generalized Gell-Mann-Oakes-Renner relation and also in Sec. V were we report numerical results.

### D Chiral symmetry and the Gell-Mann-Oakes-Renner relation

Chiral symmetry is a significant element of hadronic QCD and should be present in all realistic models. Even though our vacuum properly exhibits dynamic chiral symmetry breaking, our model Hamiltonian does indeed respect this symmetry since the commutator

$$[H,Q_5]=m𝑑\stackrel{}{x}\mathrm{\Psi }^{}(\stackrel{}{x})\mathrm{\Psi }(\stackrel{}{x})0$$

essentially vanishes, consistent with the small $`u`$, $`d`$ quark mass. Related, our RPA states also preserve chiral symmetry as the RPA meson creation operator commutes with the chiral charge in the chiral limit ($`m0`$)

$$[Q^{}(RPA),Q_5]=0.$$

However, the TDA operator, Eq. (22) above, does not commute with $`Q_5`$ and violates chiral symmetry since it is not fully symmetric in operator structure (only contains $`B^{}D^{}`$ and not $`BD`$). This can also be documented by chiral transforming the TDA meson state verifing that $`B^{}D^{}`$ rotates to combinations of $`B^{}B`$, $`DD^{}`$ and $`DB`$. Hence, the TDA ansatz is not closed under a chiral rotation and Goldstone bosons will not appear in the TDA spectrum. We therefore expect significant, but unphysical, chiral symmetry violations in the TDA calculations and anticipate the TDA pion mass to be much larger than in the RPA which is confirmed in Sec. V as only the RPA calculations yield a Goldstone pion in the chiral limit.

On the other hand, the RPA ansatz, Eq. (28), is chirally invariant since it is closed under this rotation in the chiral limit ($`|X|=|Y|`$). Hence the mechanism of spontaneous symmetry breaking is present and the pion mass is zero according to its Goldstone boson nature as we numerically verify in the next section.

With these results we now derive two different Gell-Mann-Oakes-Renner (GMOR) relations, one based upon exact model eigenstates while the other relates RPA states and the BCS vacuum. Assume we have the complete set of exact eigenstates, $`|n`$, to our QCD model Hamiltonian, including the vacuum ground state $`|\mathrm{\Omega }_{exact}`$. Evaluating the double commutator

$$\mathrm{\Omega }_{exact}|[Q_5,[Q_5,H]]|\mathrm{\Omega }_{exact}=4m\overline{q}q_{exact}$$
(39)
then generates the exact quark condensate. Evaluating the double commutator again, but now invoking twice the completeness relation, $`1=_n|nn|`$, and identifying the decay constant relation, Eq. (36), leads to the generalized GMOR relation

$$2m\overline{q}q=\underset{n}{}M_n^2f_n^2$$
(40)
summed over all (ground and excited) pseudoscalar meson states with mass $`M_n`$ and decay constants $`f_n`$.

This can be extended to flavor using two of the nine $`SU_f(3)`$ axial charge operators, $`Q_5^a`$ and $`Q_5^b`$, to derive the following GMOR relations:

* a= b = 1, 2, or 3

$$2m_u\overline{u}u=\underset{\pi ^n}{}M_{\pi ^n}^2f_{\pi ^n}^2$$
* a = b = 4, 5, 6, or 7

$$\left(\frac{m_u+m_s}{2}\right)\overline{u}u+\overline{s}s=\underset{K^n}{}M_{K^n}^2f_{K^n}^2$$
* a = b = 8,

$$\left(\frac{2}{3}m_u\overline{u}u+\frac{4}{3}m_s\overline{s}s\right)=\underset{\eta _8^n}{}f_{\eta _8^n}^2M_{\eta _8^n}^2$$
* a = b = 0

$$\left(\frac{4}{3}m_u\overline{u}u+\frac{2}{3}m_s\overline{s}s\right)=\underset{\eta _0^n}{}f_{\eta _0^n}^2M_{\eta _0^n}^2.$$

Because the exact eigenstates will generally not be available, these equations are of limited value. However, they do provide testing criterion for approximation solutions. Further, since decay constants are suppressed for excited states $`n`$ we can drop higher terms to obtain more useful relations involving ratios such as

$$\left(\frac{M_Kf_K}{M_\pi f_\pi }\right)^2=\frac{m_u+m_s}{2m_u}\left(\frac{\overline{u}u+\overline{s}s}{2\overline{u}u}\right)$$

which we will use in Sec. V in connection with our discussion of the kaon mass calculation.

Next we utilize Thouless’ theorem  applied to the chiral charge operator

$$\mathrm{\Omega }|[Q_5,[Q_5,H]]|\mathrm{\Omega }=2\underset{n}{}|\mathrm{\Psi }_{RPA}^{nJP}|Q_5|\mathrm{\Omega }_{RPA}|^2(E_nE_0)_{RPA}$$
(41)
to immediately derive the RPA GMOR relation

$$2m\overline{q}q_{BCS}=(M_\pi ^2f_\pi ^2)_{RPA}$$
(42)
where we have again dropped the excited state decay constants. Note that the left hand side entails the BCS vacuum while the right hand side involves the RPA states and energies. This relation clearly predicts the RPA pion is a Golstone boson in the chiral limit (i.e $`M_\pi 0`$ for $`m0`$) which we numerically confirm in the next section.

Finally, calculating the RPA meson mass spectrum does not require obtaining the RPA vacuum (see Eq.(29)), but computing the RPA decay constants does. Determining $`|\mathrm{\Omega }_{RPA}`$ is actually quite difficult, however, the leading correction can be approximately calculated using another theorem by Thouless

$$|\mathrm{\Omega }_{RPA}|\mathrm{\Omega }+(F^{})^2|\mathrm{\Omega }$$
(43)
where $`F^{}`$ is a TDA type operator given by

$$F^{}=\underset{\alpha \beta }{}𝑑\stackrel{}{q}f_{\alpha \beta }(\stackrel{}{q})^{}B_\alpha ^{}(\stackrel{}{q})D_\beta ^{}(\stackrel{}{q})$$

with the fermion pair operator $`B^{}D^{}`$ now obeying bosonic commutation relations. Here $`f^{}`$ is an unknown amplitude assumed to be small (this is not true in the chiral limit as discussed in the next section). The approximate RPA vacuum is thus described as a mixture of the BCS vacuum having Cooper pairs (Eq. 8) and two quasibosons (mesons) coupled to vacuum quantum numbers ($`0^{++}`$). Since, as shown in the next section, the RPA is of importance for mainly the low-lying pseudoscalar states, we assume that only these states contribute to $`|\mathrm{\Omega }_{RPA}`$ and neglect all others. Our result easily generalizes to scalar or other mesons (appropriately coupled to $`J^P=0^{++}`$). Within this approximation we obtain an improved quark condensate to be compared with Eq. (12)

$$\overline{q}q_{RPA}\overline{q}q_{BCS}+\frac{8_{\alpha \beta }𝑑\stackrel{}{q}s_q|f_{\alpha \beta }(\stackrel{}{q})|^2}{1+2_{\alpha \beta }𝑑\stackrel{}{q}|f_{\alpha \beta }(\stackrel{}{q})|^2}.$$

To obtain $`f_{\alpha \beta }`$, we impose $`Q_\pi (RPA)|\mathrm{\Omega }_{RPA}=0`$ and use Eq. (43) neglecting higher order terms corresponding to Fock states with more than two pions. We also only retain ground state meson (pion) contributions from $`X`$ and $`Y`$ yielding

$$f_{\alpha \beta }(q)=N^1Y(q)(1)^{\frac{1}{2}+\alpha +\beta }\frac{1}{2}\alpha \frac{1}{2}\beta |00$$
(44)
with a normalization constant depending on both RPA wavefunction components

$$N^2=2\frac{d\stackrel{}{q}}{(2\pi )^3}Y(q)X(q).$$

The improved RPA quark condensate is

$$\overline{q}q_{RPA}\overline{q}q_{BCS}+c_0^{\mathrm{}}s_qY(q)^2q^2𝑑q$$
(45)
where the constant c is given by

$$c=\frac{2Y(q)^2q^2𝑑q}{(Y(q)X(q)q^2𝑑q)^2+\frac{1}{2}(Y(q)^2q^2𝑑q)^2}.$$

Also the improved pion decay constant to this vacuum is

$$f_\pi \frac{1}{\pi \sqrt{(2\pi )^3M_\pi }}s_q(X(q)+c\pi ^2N^2Y(q))q^2𝑑q$$
(46)
in contrast to Eq. (38). We have found that for meson masses above 800 $`MeV`$ $`Y`$ is very small and there is no essential difference between RPA and TDA.

We now comprehensively apply the above formulas and conduct a comparative analysis of the TDA and RPA approaches.

## V Applications and numerical results

In this section we first present and discuss our TDA spectra for the pseudoscalar and vector mesons and then in subsection B we compare to our RPA results for both meson masses and the decay constants. Subsection C treats mixing for the $`\eta `$ and $`\eta ^{}`$ mesons while subsection D details applications to charmed mesons, especially $`J/\mathrm{\Psi }`$ states.

### A Tamm-Dancoff spectrum

Solving the TDA eigenvalue problem is straight forward. The diagonal part of the matrix contains the IR singularity that is rigorously controlled by cancellation, again permitting numerical regulation by simply skipping the point $`q=k`$ as in the gap equation solution. Due to the linear nature of this system, results are convergent for a mesh as sparce as 700 points.

Tables I and II summarize the resulting TDA meson spectrum corresponding to the five $`JLS`$ kernels specified in Sec. IV A. The energy difference between pseudoscalar and vector states is about 200 $`MeV`$ for the $`u`$/$`d`$ quark mesons, 80 $`MeV`$ for the open flavored ($`K`$, $`K^{}`$) and only 50 $`MeV`$ for the pure strange composites ($`\varphi `$). Our one parameter model can not accurately describe the entire observed splitting indicating additional dynamics beyond simple spin interactions from a Dirac spinor field is needed. In general, the masses are in good agreement with the PDG  accepted values for the vector mesons, but the TDA low-lying pseudoscalar sector is deficient. This is expected since in this channel vacuum (chiral symmetry) effects are most prominent. In the scalar channel the situation is more confusing, since other hadron states, some with explicit gluonic structure, can more easily mix. Since our unified model allows us to treat glueballs, mesons and hybrids comprehensively, future work will further address understanding this channel. Interestingly, our lightest $`f_0`$ mass, which has a P wave oribital excitation, is below 1 $`GeV`$. We also find that the computed TDA masses are not very sensitive to details of the vacuum that enters via the gap function characterizing the BCS ground state.

As described above, we may also use a simplified TDA equation to extract a constituent quark mass. In the chiral limit, the mass obtained is 51 $`MeV`$ and this dressing is roughly constant, consistent with $`SU_f(3)`$ symmetry, up to current masses of 150 $`MeV`$, where the generated constituent mass is 203 $`MeV`$. We are therefore dealing with light dressed quarks, even for strange quarks which may prove a deficiency when calculating electromagnetic form factors. We shall also address this issue in the near future.

Another important feature of our relativistic effective Hamiltonian is that it naturally includes the kinematical and spin dependent interactions (e.g. spin-spin, spin-orbit, tensor). These effects are very important in the light quark sector even after chiral symmetry breaking, because the light quark constituent mass generated in our scheme, around 80 $`MeV`$, is still small as compared to our interaction scale (424 $`MeV`$). In particular, notice in Table II the spin splitting between the $`0^{++}`$, $`1^{++}`$ and $`2^{++}`$ mesons, all having the same $`L`$ and $`S`$ quantum numbers (also observe the large radial excitation in each channel). The level spacing is consistent with Ref.  but very different from the naive expectations from the constituent quark model. The is due to the confining (nonperturbative) $`\stackrel{}{L}\stackrel{}{S}`$ coupling, which is only of order $`\alpha _s`$ in the constituent quark model. The spin spacing is governed by the matrix element

$$LSJ|\stackrel{}{L}\stackrel{}{S}|LSJ=\frac{1}{2}[J(J+1)L(L+1)S(S+1)]$$

which, for $`L=S=1`$, reduces to $`2+\frac{J(J+1)}{2}`$. The $`2`$ in this expression describes the light scalar meson, whereas the $`J(J+1)`$ contribution explains the splittings in Table II. The calculated $`2^{++}`$ mass is much heavier than the lightest observed $`f_2`$ at 1270 $`MeV`$. Clearly, including coupling to the $`L=3`$, $`S=1`$ channel as well as multi-quark Fock states will alter, but improve, our prediction. Further, and as mentioned previously, additional quark-gluon interaction terms should also be included in our Hamiltonian which have a different Lorentz spin structure. In our model this would generate a weaker hyperfine interaction of order $`\alpha _s`$. This would also provide additional splitting between $`S=0`$ and $`S=1`$ levels that would improve our TDA, and especially RPA, $`\pi `$-$`\rho `$ mass underprediction. Such an analysis would fully clarify the relative importance of chiral symmetry versus spin interactions as the latter is generally attributed the dominant effect in conventional constituent quark models having color magnetic, effective one gluon exchange potentials (see Ref. ). We are currently examining this issue and will report results in a future paper.

Finally and also related is the spin-orbit splitting for other flavors which is summarized in Tables I and IV. In Table IV we illustrate the familiar $`J(J+1)`$ dependence by fitting the TDA spectrum for different flavors.

### B RPA spectrum and decay constants

We now present our RPA results and compare with both TDA and observation. Table III and Figs. 2 through 7 highlight our key results. In general the RPA and TDA masses agree except for the light pseudoscalar mesons where the RPA provides a better description. This is because, as discussed above, only the RPA correctly implements chiral symmetry as illustrated in Fig. 3 where the lightest scalar and pseudoscalar masses are plotted as a continuous function of the current quark mass. Note that only the pseudoscalar mass approaches zero in the chiral limit consistent with Goldstone’s theorem (the observed pion mass is reproduced for a $`u`$ quark mass of about 2 $`MeV`$).

The kaon system reveals the largest model deficiency (see Table I). Even the RPA in the chiral limit produces a too massive kaon, about 850 $`MeV`$. As indicated by Fig. 4, to reproduce the observed kaon requires a strange quark mass of about 50 $`MeV`$. A detailed analysis reveals that the explicit contribution of the current quark mass to the RPA equation -through its appearance in $`ϵ(k)`$\- is only additive. It is the gap angle, introducing an implicit flavor dependence, which inhibits a lighter kaon mass. We could fit both pion and kaon masses by adjusting the current quark masses, but we prefer awaiting improvements from renormalizing the quark gap equation (work in progress and will be subsequently publish). The gap angle for a non-zero current quark mass is very sensitive at high momentum to the dominanting Coulomb potential and sizeable corrections are expected. Also for higher lying excited meson states there is also the issue of two particle, two hole Fock state contributions.

In Fig. 5 we compare our scalar, pseudovector and tensor meson TDA and RPA predictions to data. In general there is qualitative agreement. Note that for the higher excited states above 1 $`GeV`$ the TDA and RPA results are identical and it is clear that these systems are not governed by chiral symmetry.

We also illustrate the behavior of various wavefunctions. In Fig. 6 we detail the difference between TDA wavefunctions for the $`f_0`$ and $`\pi `$ mesons. Figure 7 compares the TDA and RPA wavefunctions for the pion.

Finally we discuss decay constants. Both the TDA and RPA pion decay constants are too small, about 17 $`MeV`$, in contrast to the observed value of 93 $`MeV`$. This is consistent with the Orsay  results. We again attribute this to the model Hamiltonian and vacuum as the BCS angle does not have sufficiently high momentum components. However, if we use the approximate RPA vacuum in the quasiboson approximation (see Eq. (4.25)) the improved RPA decay constant increases to 57 $`MeV`$. Appropriately, the condensate also significantly increases to $`(320MeV)^3`$, in much better agreement with lattice results $`(250MeV)^3`$. The latter result is consistent with applications in nuclear physics where the RPA tends to overcorrelate the ground state.

Since we have only approximately evaluated the RPA vacuum by truncating at the two pion level it is not surprising that $`f_\pi `$ does not agree with measurement (or chiral perturbation theory results) and therefore needs further refining. We expect the truncation to the two pion level, Eq. (43), to be reasonable provided we are not in the chiral limit, since then $`YX`$. However, we would like to be able to calculate in that regime. Further, the calculation in Ref.  points to a necessary decrease of the BCS condensate when including coupled channels. They argue that for a pure chiral pion, coupling additional channels would decrease its mass which might even become negative, destabilizing the vacuum. This argument seems sound and since we are above the chiral quark mass limit, at a model value $`m=5`$ $`MeV`$ where our pion is too massive (277 $`MeV`$), this decrease is a welcome improvement. We defer further discussion until publication by our collaborative effort with the Lisbon group which will also clarify pionic correlations in the ground state of our approach.

Recalling that our computed condensates require renormalization except in the chiral limit, we can only test the GMOR relation for $`m=0`$ which is trivially satisfied in the RPA. However, we note that the decay constants for excited pion states are much smaller and rapidly approach zero in the chiral limit. Hence the GMOR relation is satisfied with predominatly the first state. This is not true for heavier quarks and to numerically satisfy the GMOR requires including several eigenstates.

### C $`\eta \eta ^{}`$ mixing

An extremely challenging but still not understood problem is the $`\eta `$, $`\eta ^{}`$ system and attending flavor mixing of the light quarks. Although our effective Hamiltonian has an explicit flavor dependence through the current quark masses, it still conserves flavor. However, if the gluon sector is included, such as through the hyperfine, minimal coupling interaction, an effective flavor dependence naturally emerges though higher order quark-glue-quark effects and dynamic mixing of flavor states is possible. We are currently deriving such a term which is similar, but more rigorous than the t’Hooft interaction based upon instantons (classical glue). This will be reported in a future communication, however, it is still of interest to perform a simple $`\eta `$, $`\eta ^{}`$ mixing analysis by introducing a flavor off-diagonal interaction as we now detail.

With no dynamic flavor mixing, the $`\eta `$ and $`\eta ^{}`$ are (poorly) modeled as $`SU_f(3)`$ octet, $`\eta _8`$, and singlet, $`\eta _0`$, states respectively given by (we adopt the convention used in Refs. )

$`\eta `$ $`\eta _8=`$ $`\sqrt{{\displaystyle \frac{1}{3}}}n\overline{n}\sqrt{{\displaystyle \frac{2}{3}}}s\overline{s}=\mathrm{cos}\theta _{SU_f(3)}n\overline{n}\mathrm{sin}\theta _{SU_f(3)}s\overline{s}`$ (47)
$`\eta ^{}`$ $`\eta _0=`$ $`\sqrt{{\displaystyle \frac{2}{3}}}n\overline{n}+\sqrt{{\displaystyle \frac{1}{3}}}s\overline{s}=\mathrm{sin}\theta _{SU_f(3)}n\overline{n}+\mathrm{cos}\theta _{SU_f(3)}s\overline{s}`$ (48)
involving the isoscalar $`n\overline{n}=(u\overline{u}+d\overline{d})/\sqrt{2}`$ and $`s\overline{s}`$ states and pure $`SU_f(3)`$ mixing angle $`\theta _{SU_f(3)}=54.74^{}`$. We have computed the TDA masses of the pure $`n\overline{n}`$ and $`s\overline{s}`$ meson states to be $`M_n=`$ 612 $`MeV`$ and $`M_s=`$ 1002 $`MeV`$, respectively. Hence the predicted, pure $`SU_f(3)`$, $`\eta `$, $`\eta ^{}`$ masses are

$`M_\eta `$ $`M_{\eta _8}=`$ $`\mathrm{cos}^2\theta _{SU_f(3)}M_n+\mathrm{sin}^2\theta _{SU_f(3)}M_s=872MeV`$ (49)
$`M_\eta ^{}`$ $`M_{\eta _0}=`$ $`\mathrm{sin}^2\theta _{SU_f(3)}M_n+\mathrm{cos}^2\theta _{SU_f(3)}M_s=742MeV.`$ (50)
For the RPA, $`M_n=`$ 290 $`MeV`$ and $`M_s=`$ 978 $`MeV`$, yielding $`M_\eta =749`$ $`MeV`$ and $`M_\eta ^{}=519`$ $`MeV`$. Similarly the $`SU_f(3)`$ $`\eta `$, $`\eta ^{}`$ decay constants are given in terms of the $`u/d`$ isoscalar, $`f_n`$, and strange quark, $`f_s`$, decay constants

$`f_\eta `$ $``$ $`f_{\eta _8}=\mathrm{cos}^2\theta _{SU_f(3)}f_n+\mathrm{sin}^2\theta _{SU_f(3)}f_s`$ (51)
$`f_\eta ^{}`$ $``$ $`f_{\eta _0}=\mathrm{sin}^2\theta _{SU_f(3)}f_n+\mathrm{cos}^2\theta _{SU_f(3)}f_s.`$ (52)
Note the quadratic dependence on angles due to expressing both the axial current and meson states in their $`SU_f(3)`$ representations. Using the TDA computed values of $`f_n=17`$ $`MeV`$ and $`f_s=75`$ $`MeV`$, yields $`f_{\eta _8}=56`$ $`MeV`$ and $`f_{\eta _0}=`$ 40 $`MeV`$. The RPA values are $`f_n=57`$ $`MeV`$ and $`f_s=75`$ $`MeV`$, giving $`f_{\eta _8}=68`$ $`MeV`$ and $`f_{\eta _0}=`$ 63 $`MeV`$.

We can now improve these results by generalizing our Hamiltonian, which is diagonal in flavor space, to have off-diagonal matrix elements. The simplest prescription is to just add a constant, $`q\overline{q}|H^{}|q\overline{q}=\lambda `$, to both diagonal, $`n\overline{n}|H^{}|n\overline{n}=2\lambda `$, $`s\overline{s}|H^{}|s\overline{s}=\lambda `$, and off-diagonal $`n\overline{n}|H^{}|s\overline{s}=\sqrt{2}\lambda `$, terms giving

$$H=\left(\begin{array}{cc}M_n+2\lambda & \sqrt{2}\lambda \\ \sqrt{2}\lambda & M_s+\lambda \end{array}\right).$$
(53)
Diagonalizing $`H`$ leads to the new, mixed mass eigenvalues

$`M_\eta `$ $`=`$ $`{\displaystyle \frac{M_n+M_s+3\lambda }{2}}{\displaystyle \frac{1}{2}}\sqrt{M_n^2+M_s^2+9\lambda ^2+2(\lambda M_nM_sM_nM_s\lambda )}`$ (54)
$`M_\eta ^{}`$ $`=`$ $`{\displaystyle \frac{M_n+M_s+3\lambda }{2}}+{\displaystyle \frac{1}{2}}\sqrt{M_n^2+M_s^2+9\lambda ^2+2(\lambda M_nM_sM_nM_s\lambda )}`$ (55)
and eigenstates

$`\eta `$ $`=`$ $`\mathrm{cos}(\theta _{SU_f(3)}+\theta _P)n\overline{n}\mathrm{sin}(\theta _{SU_f(3)}+\theta _P)s\overline{s}=\mathrm{cos}\theta _P\eta _8\mathrm{sin}\theta _P\eta _0`$ (56)
$`\eta ^{}`$ $`=`$ $`\mathrm{sin}(\theta _{SU_f(3)}+\theta _P)n\overline{n}+\mathrm{cos}(\theta _{SU_f(3)}+\theta _P)s\overline{s}=\mathrm{cos}\theta _P\eta _0+\mathrm{sin}\theta _P\eta _8`$ (57)
involving rotation by an additional angle $`\theta _P`$ that is a function of $`\lambda `$. The mixed, presumably more physical, decay constants are then

$`f_\eta `$ $`=`$ $`\mathrm{cos}\theta _Pf_{\eta _8}\mathrm{sin}\theta _Pf_{\eta _0}`$ (58)
$`f_\eta ^{}`$ $`=`$ $`\mathrm{cos}\theta _Pf_{\eta _0}+\mathrm{sin}\theta _Pf_{\eta _8}.`$ (59)
Performing a least squares fit to the observed masses ($`M_\eta =`$ 547 $`MeV`$, $`M_\eta ^{}=`$ 958 $`MeV`$), yields $`\lambda =33MeV`$ ($`\theta _P=`$ $`61^{}`$) for the TDA, which in turn produces $`M_\eta =`$ 541 $`MeV`$, $`M_\eta ^{}=`$ 974 $`MeV`$, $`f_\eta =`$ 62 $`MeV`$ and $`f_\eta ^{}=`$ -30 $`MeV`$. For the RPA, $`\lambda =`$ 82 $`MeV`$ ($`\theta _P=`$ $`44^{}`$) generating $`M_\eta =`$ 433 $`MeV`$, $`M_\eta ^{}=`$ 1081 $`MeV`$, $`f_\eta =`$ 93 $`MeV`$ and $`f_\eta ^{}=`$ -3 $`MeV`$. It is interesting that while the simple mixing provides improvement, the TDA masses are in better agreement than the RPA. We could also improve the decay constants utilizing a two angle mixing formalism  but refrain since clearly a more sophisticated treatment is necessary which will be provided by our quark-gluon coupling formulation in the near future.

### D Heavy mesons

The constituent quark models and non-relativistic expansions of QCD offer more reliable results for heavy quark systems where physical intuition from Quantum Mechanical bound states is more appropriate. Hence the charmed mesons afford a good limiting testing for our relativistic approach. We again calculate the spectrum of the charmed mesons (see Figures ) and of charmonium (Fig. ) using our many-body model. The TDA is now sufficient since chiral symmetry, crucial for the light mesons, is not a constraint and the RPA will produce the same results. Using a charmed quark mass of 1200 $`MeV`$, the general features of the spectra are well reproduced and the radial excitations, $`\mathrm{\Psi }(2S)`$ and $`\eta _c`$, are adequately described. We therefore expect our predictions for the remaining unconfirmed states to be reasonable.

The angular momentum splittings of these systems are known to be dominated by the one gluon exchange potential (OGE) which we have not included in this calculation. Hence there will be improvement from future calculations based upon our renormalized project. A general feature reflected by our charmed spectra is the near vanishing of the spin-spin interaction, leading to degenerate $`0^+`$ and $`1^{}`$ states. Since the 100-150 $`MeV`$ hyperfine splittings in these systems will presumably be recovered when we include the perturbative OGE, we do not comment further. As for the spin-orbit splittings, our results are too large for the $`D`$ mesons and too small for the $`\chi _c`$ mesons. These splittings are adequately explained in non-relativistic quark models (see Ref. ) where the absence of a large spin-orbit effect for light quark masses is attributed to the cancellation between the Thomas precession in the confining potential and the one gluon exchange effective potential, although it is a bit concerning that the actual splitting between the $`\chi _{c0}`$ and the $`\chi _{c1}`$ is twice the size of the $`\chi _{c1}`$-$`\chi _{c2}`$ splitting, when according to the spin-orbit $`J(J+1)`$ rule, it should be a half.

Finally, we note that by including D waves in the charmonium spectrum we are able to resolve the long standing ”overpopulation” problem of $`J/\mathrm{\Psi }`$ states relative to observation. This is clearly illustrated in Fig. 10. This (previous) deficiency in number of $`c\overline{c}`$ states has been characterized as evidence for glueballs and/or hybrid mesons because the $`J/\mathrm{\Psi }`$ system is believed to be gluon rich. Our result suggests, however, that simple level counting may not be effective in identifying hadrons with explicit gluonic degrees of freedom.

## VI Outlook and summary

Before summarizing our results we comment on the strengths and weaknesses of our many-body approach as well as some attending, open hadronic physics issues. Beginning with our Hamiltonian, $`H`$, the current-current (density-density) color interaction forbids free, isolated colored objects in the theory. For the vacuum this is realized in the BCS by an infinite shift in the free quark self-energy due to the integral of $`\widehat{V}(|\stackrel{}{k}\stackrel{}{q}|)`$. Similarly for hadrons in both TDA and RPA, colored composite objects (e.g. diquarks) are precluded by the appearance of $`\widehat{V}(0)`$ which is divergent, whereas in the singlet channel this divergence is removed by vanishing color factors. Next we note that $`H`$ conserves chiral symmetry yet our BCS vacuum properly exhibits dynamic chiral symmetry breaking. Further, our RPA pion emerges as a Goldstone boson in the chiral limit. This is not true for the TDA pion since only the RPA excitation operator commutes with the chiral charge. Another, significant model feature is that Fock state truncation is a controlable approximation amendable to systematic improvement. Thus our Hamiltonian many-body approach is an attractive, promising method for comprehensively investigating hadronic structure as it embodies confinement, chiral symmetry breaking and orderly construction of multiparticle Fock states. It also provides an excellent vehicle for testing more fundamental effective Hamiltonians as well as affording a powerful phenomenological framework for hadron structure. Considering the form of $`H`$, with only a single predetermined dynamical parameter, it is encouraging that the chiral limit is adequately reached in the RPA and that the meson spectrum is in qualitative agreement with experiment. To achieve detailed, quantitative descriptions will require further improvements in both the Hamiltonian and effects from including higher Fock space components. In particular, both the high energy behavior (Coulomb potential) and quark-gluon coupling effects (efectively instantons) will be incorporated and reported in a future publication. Finally, our current study is similar, but more extensive than the Orsay analysis  due to our application to multiflavor systems. Our results are also more realistic (and numerically more difficult) then that work since we have utilized a linear confining interaction, determined by lattice and Regge phenomenology, which generates complicated nonlocal integral equations, rather than solving a simpler differential equation for a harmonic oscillator potential.

Summarizing, we have performed approximate, but large-scale diagonalizations of an effective Coulomb gauge Hamiltonian utilizing standard many-body techniques. Using the BCS, a non-linear gap equation has been derived and accurately solved to provide vacuum properties (quasiparticles and condensates). Incorporating only predetermined parameters (string tension, $`\sigma `$, and reasonable current quark masses), we have qualitatively reproduced the low energy $`u,d,s`$ and $`c`$ meson spectra. Most importantly, we have obtained a chiral pion, detailed that chiral symmetry is responsible for the large $`\pi \rho `$ mass splitting and resolved the problem of overpopulation of theoretical $`J/\mathrm{\Psi }`$ states.

Future work will address the full Hamiltonian in the combined quark and gluon sectors. In particular, we will obtain improved spin (hyperfine) and flavor (t’Hooft) interactions from quark-glue coupling. This should provide a better description of vacuum properties and the scalar/pseudoscalar masses, especially the $`\pi ,K,\eta `$ and $`\eta ^{}`$. We will also include more complex 2 quasiparticle-2 quasihole Fock states for heavier $`u/d`$ mesons as well as 3 quasiparticles for baryons and hybrids. Much of this work is in progress and will soon be reported.

###### Acknowledgements.

The authors are grateful for comments and discussions with P. Bicudo, J. E. Ribeiro, A. Szczepaniak and members (present and former) of the NCSU theory group. This work is supported in part by grants DOE DE-FG02-97ER41048 and NSF INT-9807009. F. J. Llanes-Estrada was a SURA-Jefferson Laboratory graduate fellowship recipient. Supercomputer time from NERSC is also acknowledged.

## A BCS Vacuum State

We further discuss the relation between the BCS rotated, $`|\mathrm{\Omega }`$, and the trivial or perturbative, $`|0`$, vacua. We first note that the BCS vacuum state given by Eq. (8) is not a unitary transformation and does not have a finite normalization. This is because the operator in the exponential is not antihermitian. It is therefore necessary to normalize matrix elements by dividing with $`\mathrm{\Omega }|\mathrm{\Omega }`$ and this is implicit in our presentation. Alternatively, and equivalently, $`|\mathrm{\Omega }`$ can be represented by a norm preserving unitary transformation of the form

$$|\mathrm{\Omega }=e^{A^{}A}|0$$

where

$$A^{}=\underset{\lambda _1\lambda _2}{}𝑑\stackrel{}{k}\mathrm{tan}\theta _kM_{\lambda _1\lambda _2}b_{\lambda _1}^{}(\stackrel{}{k})d_{\lambda _2}^{}(\stackrel{}{k})$$

and all flavor and color indices are suppressed. Here $`M_{\lambda _1\lambda _2}`$ are matrix elements of the Pauli matrices $`M=\stackrel{}{\sigma }\widehat{k}`$

$`M_{\lambda _1\lambda _2}={\displaystyle \underset{\mu }{}}(1)^\mu \widehat{k}^\mu \sigma _{\lambda _1\lambda _2}^\mu `$ (A1)
$`=\sqrt{8\pi }(1)^{\frac{1}{2}\lambda _2}{\displaystyle \underset{\mu \nu }{}}{\displaystyle \frac{1}{2}}\lambda _1{\displaystyle \frac{1}{2}}\lambda _2|1\nu 1\mu 1\nu |00Y_1^\mu (\widehat{k}).`$ (A2)
It is interesting to note that the BCS vacuum state is orthogonal to the trivial vacuum, $`0|\mathrm{\Omega }=0`$, in the infinite volume limit. Further, the Hilbert space vectors constructed from the two different vacua are also orthogonal provided the BCS angle is nonzero (vacuum condensates are present). Because of this property the BCS rotation has been called a pseudounitary transformation . Consult this reference for further details (note they have a different phase convention and use $`M=\stackrel{}{\sigma }\widehat{k}(i\sigma _2)`$).

## B General TDA Equation

Using the phase convention of Ref. , the general TDA meson equation for arbitrary angular momentum is

$`(M_{nJP}2ϵ_k)\mathrm{\Psi }_{\mathrm{\Lambda }\mathrm{\Sigma }}^{nJP}(k)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{q^2dq}{\pi ^2}}[{\displaystyle \frac{(1+s_k)(1+s_q)}{12}}\widehat{V}_\mathrm{\Lambda }(k,q)\mathrm{\Psi }_{\mathrm{\Lambda }\mathrm{\Sigma }}^{nJP}(q)+`$
$`+{\displaystyle \underset{lLSf}{}}{\displaystyle \frac{c_kc_q}{2}}\mathrm{\Psi }_{LS}^{nJP}(q)\widehat{V}_l(k,q){\displaystyle \underset{1}{}}`$ $`+`$ $`{\displaystyle \underset{lLSfghL_1L_2}{}}3(1s_q)(1s_k)\widehat{V}_l(k,q)\mathrm{\Psi }_{LS}^{nJP}(q){\displaystyle \underset{2}{}}]`$
where the angular momentum products are

$`{\displaystyle \underset{1}{}}=10L0|l010\mathrm{\Lambda }0|l0(2f+1)\sqrt{(2\mathrm{\Sigma }+1)(2S+1)(2\mathrm{\Lambda }+1)(2L+1)}`$
$`W(\mathrm{\Sigma }{\displaystyle \frac{1}{2}}S{\displaystyle \frac{1}{2}};{\displaystyle \frac{1}{2}}f)W({\displaystyle \frac{1}{2}}1{\displaystyle \frac{1}{2}}1;{\displaystyle \frac{1}{2}}f)W(L1\mathrm{\Lambda }1;lf)W(LS\mathrm{\Lambda }\mathrm{\Sigma };Jf)(1)^{J+L+1}(1+(1)^{S+\mathrm{\Sigma }})`$
and

$`{\displaystyle \underset{2}{}}=(2f+1)(2g+1)(2h+1)\sqrt{(2L_1+1)(2L_2+1)(2L+1)(2\mathrm{\Lambda }+1)(2\mathrm{\Sigma }+1)(2S+1)}`$
$`(1)^{(1+2(fh)+L\mathrm{\Lambda }+lJ)}1010|L_101010|L_20L_10\mathrm{\Lambda }0|0L_20L0|l0`$
$`W({\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}1L_1;1f)W(f{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}};1\mathrm{\Sigma })W(L_21{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}};fg)W({\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{2}}g{\displaystyle \frac{1}{2}};1S)`$
$`W(L_2LL_1\mathrm{\Lambda };lh)W(SL\mathrm{\Sigma }\mathrm{\Lambda };Jh)W(f\mathrm{\Sigma }gS;{\displaystyle \frac{1}{2}}h)W(L_2gL_1f;{\displaystyle \frac{1}{2}}h).`$
These formulas have been applied to several different meson spin, parity states as presented in Sec. IV.

The moments of the angular integrations for the linear potential are obtained from

$$\widehat{V}_n^L=8\pi \sigma _1^1\frac{1}{|\stackrel{}{k}\stackrel{}{q}|^4}x^n𝑑x$$

where $`x=\widehat{k}\widehat{q}`$. These can be calculated using the recurrence relation

$$\widehat{V}_n^L=\frac{k}{q}\widehat{V}_{n1}^L\frac{4\pi \sigma }{q}\frac{d}{dk}\left(\frac{1}{2qk}_1^1\frac{x^{n1}dx}{x\frac{k^2+q^2}{2qk}}\right)$$

or by explicit integration

$$\widehat{V}_n^L=\frac{k}{q}\widehat{V}_{n1}^L\frac{4\pi \sigma }{2q^2}\frac{d}{dk}\left[\frac{1}{k}\left(\frac{x^{n1}}{n1}+\frac{\omega x^{n2}}{n2}+\mathrm{}+x\omega ^{n2}\right)_1^1+\omega ^{n1}\mathrm{log}\left(\frac{\omega 1}{\omega +1}\right)\right]$$

where $`\omega \frac{k^2+q^2}{2qk}`$. Evaluating the first three moments yields

$$\widehat{V_0}^L=\frac{16\pi \sigma }{(k^2q^2)^2}$$

$$\widehat{V_1}^L=\frac{2\pi \sigma }{k^2q^2}\left[ln\left(\frac{k+q}{kq}\right)^2+(k^2+q^2)\left(\frac{4qk}{(k^2q^2)^2}\right)\right]$$

$$\widehat{V_2}^L=3\pi \sigma \frac{k^2+q^2}{k^3q^3}ln\left(\frac{k+q}{kq}\right)^2\frac{8\pi \sigma }{k^2q^2}\frac{k^4+q^4}{(k^2q^2)^2}.$$

For the Coulomb potential, $`V_C=\frac{\alpha _s}{|\stackrel{}{x}\stackrel{}{y}|}`$, $`\widehat{V}_C=\frac{4\pi \alpha _s}{|\stackrel{}{k}\stackrel{}{q}|^2}`$, and the angular integrals are

$$\widehat{V_0}^C(k,q)=\frac{2\pi \alpha _s}{qk}\mathrm{ln}\left(\left(\frac{k+q}{kq}\right)^2\right)$$

$$\widehat{V_1}^C=\frac{4\pi \alpha _s}{qk}+\frac{k^2+q^2}{2qk}\widehat{V}_0^C$$

$$\widehat{V_2}^C=\frac{k^2+q^2}{2qk}\widehat{V}_1^C.$$


format: text

text:
# Hadronic Loop Corrections to the Muon Anomalous Magnetic Moment

## Abstract

The dominant theoretical uncertainties in both, the anomalous magnetic moment of the muon and the value of the electromagnetic coupling at the $`Z`$ scale, $`M_Z`$, arise from their hadronic contributions. Since these will ultimately dominate the experimental errors, we study the correlation between them, as well as with other fundamental parameters. To this end we present analytical formulas for the QCD contribution from higher energies and from heavy quarks. Including these correlations affects the Higgs boson mass extracted from precision data.

preprint:

The magnetic moment of the electron, $`g_e`$, provides the best determination of the fine structure constant, but is currently not measured precise enough to give a sensitive probe of electroweak physics which is suppressed by a factor $`(\alpha /\pi )m_e^2/M_W^2𝒪(10^{13})`$. Meaningful (but still weak) bounds on new physics contributions can be set only if they enter without loop suppression. On the other hand, electroweak contributions to the anomalous magnetic moment of the muon , $`a_\mu =(g_\mu 2)/2`$, are enhanced by a factor $`m_\mu ^2/m_e^24\times 10^4`$, which renders them sizable enough to be detectable at the ongoing E821 experiment at the AGS at BNL. ($`a_\tau `$ has not yet been observed experimentally). E821 already reduced the experimental error to $`\pm 1.6\times 10^9`$ . The anticipated final error of about $`\pm 0.4\times 10^9`$ will mean a factor of 20 improvement relative to previous results . $`a_\mu `$ therefore provides a good laboratory to test the Standard Model (SM) and probe theories beyond it . For example, scenarios of low-energy supersymmetry with large $`\mathrm{tan}\beta `$ and moderately light superparticle masses can give large contributions to $`a_\mu `$ .

Unfortunately, the interpretation of $`a_\mu `$ is compromised by a large theoretical uncertainty introduced by hadronic effects. Two and three loop vacuum diagrams containing light quark loops cannot be calculated reliably in perturbative QCD (PQCD). Instead they are obtained by computing dispersion integrals over measured (at low energies) and theoretical (at higher energies) hadronic cross sections. At two loop ,

$$a_\mu (\text{had; 2-loop})=\left(\frac{\alpha m_\mu }{3\pi }\right)^2_{4m_\pi ^2}^{\mathrm{}}\frac{ds}{s^2}\widehat{K}\left(s\right)R(s),$$
(1)
where $`R(s)`$ is the cross section of $`e^+e^{}\text{ hadrons}`$, normalized to the tree level cross section of $`e^+e^{}\mu ^+\mu ^{}`$, and where

$$\widehat{K}(s)=_0^1𝑑x\frac{3x^2(1x)}{1x+x^2m_\mu ^2/s}.$$
(2)
The uncertainty introduced by this procedure is significantly larger than the anticipated experimental one. For example, Ref.  quotes an error of $`\pm 0.67\times 10^9`$, while other evaluations  give larger uncertainties. An analogous uncertainty occurs in the QED coupling constant, $`\widehat{\alpha }(\mu )`$, ($`\overline{\mathrm{MS}}`$ quantities will be marked by a caret), preventing its precise theoretical computation from the fine structure constant, $`\alpha `$, for $`\mu \stackrel{>}{_{}}2m_{\pi ^0}`$. Knowledge of $`\widehat{\alpha }(M_Z)`$ is indispensable for the extraction of the Higgs boson mass, $`M_H`$, from the mass of the $`W`$ boson, $`M_W`$, and the weak mixing angle, $`\widehat{s}_W^2`$. Again one must rely on a dispersion integral to be taken over $`R(s)`$, but using a different kernel function. As a result, these hadronic uncertainties are strongly correlated with each other, and also with other fundamental SM parameters, such as the strong coupling constant and the heavy quark masses.

To address these correlations we obtain compact analytical expressions for $`a_\mu (\text{had})`$ wherever possible, i.e. for the charm and bottom quarks, as well as for the high energy integration region (above $`1.8`$ GeV) for light quarks where PQCD appears to be applicable. A similar idea has been pursued in Ref. , for the calculation of $`\widehat{\alpha }(M_Z)`$.

To lowest order in QED, the anomalous magnetic moments are lepton universal and mass independent ,

$$a_{\mathrm{}}^{1\mathrm{loop}}=\frac{\alpha }{2\pi }.$$
(3)
To two-loop order one has,

$$a_\mu ^{2\mathrm{loop}}=a_\mu ^\mu +a_\mu ^e+a_\mu ^\tau +a_\mu ^b+a_\mu ^c+a_\mu ^{uds}+a_\mu ^{\mathrm{had}},$$
(4)
where the lepton contributions  amount to $`a_\mu ^\mu +a_\mu ^e+a_\mu ^\tau =4132.18\times 10^9`$. $`a_\mu ^b`$ and $`a_\mu ^c`$ are heavy quark contributions. To proceed, we used the expansion,

$$\widehat{K}(s)=1+k_1(\sqrt{s})\frac{m_\mu ^2}{s}+k_2(\sqrt{s})\frac{m_\mu ^4}{s^2}+𝒪(\frac{m_\mu ^6}{s^3}),$$
(5)
where

$$k_1(x)=\frac{25}{4}+3\mathrm{ln}\frac{m_\mu ^2}{x^2},k_2(x)=\frac{291}{10}+18\mathrm{ln}\frac{m_\mu ^2}{x^2}.$$
(6)
In Eq. (1) the integration is over the imaginary part of the photon polarization function which is related to $`R`$ by $`R(s)=12\pi \mathrm{Im}\mathrm{\Pi }(q^2+iϵ)`$ and is taken along the real axis from the quarkonium threshold to infinity. Analytic continuation allows us to integrate instead over the full function, $`\mathrm{\Pi }(s)`$, along a circle of radius $`s=\widehat{m}_q^2(\widehat{m}_q^2)`$ counterclockwise around the origin. This avoids a complicated integration over resonances, i.e. a region where perturbative QCD is not applicable, at the expense of introducing a new uncertainty from the imperfect knowledge of the $`\overline{\mathrm{MS}}`$ quark mass (see the discussion in the next paragraph). With $`\mathrm{\Pi }(s)`$ known up to $`𝒪(\alpha _s^2)`$ , we find for a quark of charge $`Q_q`$,

$$a_\mu ^q=\frac{Q_q^2}{4}\frac{\alpha ^2}{\pi ^2}\{\frac{m_\mu ^2}{4\widehat{m}_q^2}[\frac{16}{15}+\frac{3104}{1215}a_s+a_s^2\times $$

$$(0.5099+\frac{2414}{3645}n_l)]+\frac{m_\mu ^4}{16\widehat{m}_q^4}[\frac{108}{1225}0.1943a_s+$$
(7)
$$3\mathrm{ln}\frac{m_\mu ^2}{\widehat{m}_q^2}(\frac{16}{35}+\frac{15728}{14175}a_s+1.4123a_s^2+\frac{290179}{637875}n_la_s^2)]\},$$

where $`a_s=\alpha _s(\widehat{m}_q)/\pi `$, $`n_l=3`$ for charm, and $`n_l=4`$ for bottom. Higher orders in $`m_\mu ^2/\widehat{m}_q^2`$ can safely be neglected. Terms of $`𝒪(\alpha _s^2m_\mu ^4)`$ can also be dropped unless they are logarithmically enhanced. Non-perturbative effects in the operator product expansion were also computed and found to be negligible. Using (here and in the following) $`\alpha _s(M_Z)=0.120`$, $`\widehat{m}_c=1.31`$ GeV, and $`\widehat{m}_b=4.24`$ GeV, we find $`a_\mu ^c=1.39\times 10^9`$ and $`a_\mu ^b=0.03\times 10^9`$.

In contrast to the numerical integration over resonances, Eq. (7) is a simple and transparent representation of the heavy quark contribution. More importantly, the uncertainty implied by Eq. (7) is smaller. To see this, notice first the excellent behavior of its $`\alpha _s`$ expansion implying a small truncation error which we evaluated to be $`\pm 0.05\times 10^9`$. The dominant uncertainty is induced via the $`\overline{\mathrm{MS}}`$ charm mass, but it is necessarily smaller than in the conventional approach regardless of the available data in the resonance region: if one compares the two treatments — integration over the real axis heavily relying on experimental results and integration over the circle contour relying only on PQCD — one can determine the heavy quark masses by demanding consistency. This effectively amounts to deriving a specific QCD sum rule. The point we are making here, is that this is only one of a large number of possible sum rules, but not the one which uses the available information most efficiently. The charm mass extracted from the most efficient sum rule  (and other precise quark mass determinations) can then be used for Eq. (7). An unweighted average of various determinations yields  $`\mathrm{\Delta }\widehat{m}_c=\pm 0.07`$ GeV and $`\mathrm{\Delta }a_\mu ^c=\pm 0.16\times 10^9`$. In the near future these uncertainties are likely to reduce to about $`\pm 0.04`$ GeV and $`\pm 0.09\times 10^9`$, respectively. Similarly, $`\mathrm{\Delta }\widehat{m}_b=\pm 0.11`$ GeV, but the induced error in $`a_\mu `$ is negligible.

The remaining terms in Eq. (4) are due to $`u`$, $`d`$, and $`s`$ quark effects, which we separated into the contributions from $`\sqrt{s}\mu _0=1.8`$ GeV ($`a_\mu ^{uds}`$) and $`\sqrt{s}\mu _0`$ ($`a_\mu ^{\mathrm{had}}`$). $`a_\mu ^{uds}`$ can be written as an expansion in $`m_\mu ^2/\mu _0^2`$. For the leading contribution we find,

$`a_\mu ^{uds;\mathrm{LO}}={\displaystyle \frac{2}{9}}{\displaystyle \frac{\alpha ^2}{\pi ^2}}{\displaystyle \frac{m_\mu ^2}{\mu _0^2}}\left(1+B_1+{\displaystyle \underset{n=2}{\overset{\mathrm{}}{}}}d_nB_n\right),`$ (8)
where $`d_2=299/249\zeta (3)`$ and $`d_3=58057/288779\zeta (3)/4+75\zeta (5)/2`$ are coefficients of the Adler $`D`$-function, and

$`B_n={\displaystyle \frac{1}{2\pi i}}{\displaystyle _{|s|=\mu _0^2}}{\displaystyle \frac{ds}{s}}\left(1{\displaystyle \frac{\mu _0^2}{s}}\right)\left[{\displaystyle \frac{\alpha _s(s)}{\pi }}\right]^n,`$ (9)
which we compute with a 4-loop renormalization group improvement. For a representative value $`\alpha _s(\mu _0)/\pi =0.1`$,

$`B_1=7.069\times 10^2,B_2=4.514\times 10^3,`$ (10)
$`B_3=2.562\times 10^4,B_4=1.243\times 10^5.`$ (11)
Notice, that $`B_4`$ is small enough, that even with the fourth order coefficient, $`d_4`$, unknown this treatment keeps the truncation error at a negligible level. (For an estimate of the uncertainty of $`d_4`$, see Ref. .) Denoting $`k_1^{}=k_1(\mu _0)\frac{3}{2}`$ and $`k_2^{}=k_2(\mu _0)6`$, we find for the subleading contributions,

$`{\displaystyle \frac{\pi ^2}{\alpha ^2}}a_\mu ^{uds;\mathrm{rem}}={\displaystyle \frac{2}{9}}a_s{\displaystyle \frac{m_\mu ^2\widehat{m}_s^2(\mu _0^2)}{\mu _0^4}}+a_s^2{\displaystyle \frac{m_\mu ^2}{\mu _0^2}}G\left({\displaystyle \frac{\mu _0^2}{\widehat{m}_c^2(\widehat{m}_c^2)}}\right)`$ (12)
$`+\left\{(1+a_s){\displaystyle \frac{k_1^{}}{9}}+a_s^2\left[\left({\displaystyle \frac{34}{27}}\zeta (3)\right)k_1^{}+{\displaystyle \frac{3}{16}}\right]\right\}{\displaystyle \frac{m_\mu ^4}{\mu _0^4}}`$ (13)
$`+{\displaystyle \frac{2}{3}}\left\{(1+a_s){\displaystyle \frac{k_2^{}}{9}}+a_s^2\left[\left({\displaystyle \frac{281}{216}}\zeta (3)\right)k_2^{}+{\displaystyle \frac{1}{2}}\right]\right\}{\displaystyle \frac{m_\mu ^6}{\mu _0^6}},`$ (14)
where in $`𝒪(\widehat{\alpha }_s)`$ we kept the small $`s`$ quark mass effect ($`3\times 10^{12}`$). $`G(x)`$ arises from virtual charm quark effects inside a light quark loop (double bubble diagram). Even below threshold it can be well approximated as an expansion in $`x`$ despite the fact that $`\widehat{m}_c^2<\mu _0^2`$,

$$G(x)\frac{x}{1215}\left[\frac{3503}{75}\frac{2\pi ^2}{3}\frac{88}{5}\mathrm{ln}x+2\mathrm{ln}^2x\right]$$

$$+\frac{x^2}{11340}\left[\frac{1723}{420}\mathrm{ln}x\right].$$
(15)
$`G(x)`$ also applies to $`b`$ quarks, but this contribution can be safely neglected. We find $`a_\mu ^{uds}=4.38\times 10^9`$ and the leading order resummation in Eq. (8) renders the truncation error negligible.

Leading order non-perturbative contributions due to gluon and light quark condensates are given by,

$$a_\mu ^{uds;\mathrm{NP}}=\frac{2\alpha ^2}{81}a_s^2\frac{m_\mu ^2}{\mu _0^2}\left[\frac{11}{4}\frac{<a_sGG>}{\mu _0^4}\frac{<m_s\overline{s}s>}{\mu _0^4}\right]$$
(16)
where the condensates are of order $`\mathrm{\Lambda }_{\mathrm{QCD}}^4`$ and $`m_K^2f_\pi ^2`$, respectively. These terms are suppressed by two powers of $`\alpha _s`$ and they change $`a_\mu `$ by less than $`10^{12}`$, an effect completely negligible. Effects from up and down quark condensates are suppressed by a further factor of $`m_\pi ^2/m_K^2`$, and quartic mass terms are tiny, as well. Uncertainties from non-perturbative effects not accounted for by the operator product expansion, which are due to the transition from the data region to PQCD at $`\mu _0`$, have been estimated in Ref.  which quotes $`\pm 0.024\times 10^9`$.

We take the low energy contribution from Ref. ,

$$a_\mu ^{\mathrm{had}}=(63.43\pm 0.60)\times 10^9,$$
(17)
which includes a QCD sum rule improvement, PQCD down to a relatively low $`\mu _0=1.8`$ GeV, as well as additional information from $`\tau `$-decays. The quoted error is not uncontroversial  and needs to be confirmed. Note, however, that inclusion of the $`\tau `$-decay data decreases the difference between the SM prediction and the current experimental result . Our reason to use Ref.  is that it quantifies the correlation (69%) with the corresponding result on $`\mathrm{\Delta }\alpha _{\mathrm{had}}=(56.53\pm 0.83)\times 10^4`$. We also account for the almost perfect anti-correlation with higher order hadronic uncertainties, as will be discussed below. The total two-loop quark contribution is $`a_\mu ^b+a_\mu ^c+a_\mu ^{uds}+a_\mu ^{\mathrm{had}}=69.23\times 10^9`$.

This completes the discussion of the terms appearing in Eq. (4). The leading contribution at $`𝒪(\alpha ^3)`$ is from the light-by-light diagram containing an electron loop ,

$$a_\mu ^{\mathrm{e};\mathrm{lbl}}=\frac{\alpha ^3}{\pi ^3}\left[\frac{59}{3}\zeta (4)+4\zeta (2)\left(\mathrm{ln}\frac{m_\mu }{m_e}5\right)3\zeta (3)+\frac{2}{3}\right]$$

where we displayed only the leading term in explicit form, but the terms suppressed by one or two powers of $`m_e/m_\mu `$ are significant, as well. The contribution due to $`a_\mu ^{\mathrm{e};\mathrm{lbl}}=20.9479\alpha ^3/\pi ^3=262.54\times 10^9`$ is almost four times larger than the entire hadronic contribution at $`𝒪(\alpha ^2)`$. The corresponding contribution involving a muon loop is $`a_\mu ^{\mu ;\mathrm{lbl}}=0.3710\alpha ^3/\pi ^3=4.65\times 10^9`$, and the $`\tau `$ contributes ,

$$a_\mu ^{\tau ;\mathrm{lbl}}=\frac{\alpha ^3}{\pi ^3}\frac{m_\mu ^2}{m_\tau ^2}\left[\frac{3}{2}\zeta (3)\frac{19}{16}\right]=0.03\times 10^9.$$
(18)
An evaluation of the hadronic light-by-light contribution yields ,

$$a_\mu ^{\mathrm{had};\mathrm{lbl}}=(0.792\pm 0.154)\times 10^9,$$
(19)
which is consistent with the finding in Ref. . Notice, hadrons and leptons contribute with opposite signs. Since these results are based on model calculations, an independent confirmation would be desirable.

The purely leptonic 3-loop vacuum polarization contribution is again dominated by electron loops and given by $`a_\mu ^{\mathrm{};\mathrm{vpol}}=2.7294\alpha ^3/\pi ^3=34.21\times 10^9`$. On the other hand, the contribution from two hadronic loops is suppressed by a factor $`m_\mu ^4/16m_{\pi ^\pm }^4`$ and therefore small, $`a_\mu ^{\mathrm{had};\mathrm{vpol}}=(0.027\pm 0.001)\times 10^9`$ . Similar to our strategy at two loops, we separated the mixed leptonic-hadronic contribution into heavy quarks, light quarks ($`\sqrt{s}\mu _0`$), and light hadrons ($`\sqrt{s}\mu _0`$). The $`\tau `$-hadronic contribution is of the order of $`10^{12}`$ and negligible. As for the other leptons, we use the kernel functions from Ref.  and obtain for charm and bottom quarks,

$$a_\mu ^{\mathrm{}q;\mathrm{vpol}}=\frac{\alpha ^3}{\pi ^3}Q_q^2\frac{m_\mu ^2}{4\widehat{m}_q^2}[(\frac{8}{5}+\frac{1552}{405}a_s)$$

$$(\frac{110}{27}\frac{\pi ^2}{3}+\frac{23}{36}\mathrm{ln}\frac{m_\mu ^2}{\widehat{m}_q^2}+\frac{1}{9}\mathrm{ln}\frac{m_\mu ^2}{m_e^2})\frac{1771}{675}],$$
(20)
which amounts to $`0.05\times 10^9`$ and $`0.002\times 10^9`$, respectively. For the light quarks we obtain,

$$a_\mu ^{\mathrm{}uds;\mathrm{vpol}}=\frac{\alpha ^3}{\pi ^3}\frac{1+a_s}{9}[\frac{m_\mu ^2}{\mu _0^2}(\frac{371}{9}4\pi ^2+\frac{23}{3}\mathrm{ln}\frac{m_\mu ^2}{\mu _0^2}+$$

$$\frac{4}{3}\mathrm{ln}\frac{m_\mu ^2}{m_e^2})+\frac{m_\mu ^4}{\mu _0^4}(\frac{20359}{576}\frac{103}{24}\pi ^2+\frac{509}{72}\mathrm{ln}\frac{m_\mu ^2}{\mu _0^2}+$$
(21)
$$\frac{19}{6}\mathrm{ln}\frac{m_\mu ^2}{m_e^2}\frac{5}{24}\mathrm{ln}^2\frac{m_\mu ^2}{\mu _0^2}2\mathrm{ln}\frac{m_\mu ^2}{\mu _0^2}\mathrm{ln}\frac{m_\mu ^2}{m_e^2})]=0.15\times 10^9.$$

We found the vacuum polarization contribution arising from an electron in one loop and light hadrons in the other one, $`a_\mu ^{e\mathrm{had};\mathrm{vpol}}=0.97\times 10^9`$, by constructing a simplified function $`R(s)`$ which reproduces the results in Ref. . $`a_\mu ^{e\mathrm{had};\mathrm{vpol}}`$ is almost completely ($`>99.9\%`$) correlated with $`a_\mu ^{\mathrm{had}}`$ in Eq. (17). On the other hand, $`a_\mu ^{\mu \mathrm{had};\mathrm{vpol}}=1.80\times 10^9`$ is very strongly ($`97\%`$) anti-correlated with $`a_\mu ^{\mathrm{had}}`$. The small uncorrelated error contributions are clearly negligible, and the correlated ones can be added (subtracted) linearly, slightly reducing the error in Eq. (17) to $`\pm 0.59\times 10^9`$. Including the uncertainty from the light-by-light contribution in Eq. (19) we obtain $`\pm 0.61\times 10^9`$ as the total hadronic error excluding parametric uncertainties.

Taking the other theoretical uncertainties mentioned earlier to be 100% correlated with the corresponding uncertainties in $`\mathrm{\Delta }\alpha _{\mathrm{had}}`$ (which is a fit parameter), we obtain a residual correlation of 55% or $`\pm 0.34\times 10^9`$ (the uncorrelated error is $`\pm 0.51\times 10^9`$).

The four-loop contribution , $`a_\mu ^{4\mathrm{loop}}=`$ $`126.04\alpha ^4/\pi ^4`$ $`=3.67\times 10^9`$, and the five-loop estimate , $`a_\mu ^{5\mathrm{loop}}=930\alpha ^5/\pi ^5=0.06\times 10^9`$, are also included.

One-loop electroweak corrections due to $`W`$ and $`Z`$ boson loops are given by ,

$$a_\mu ^{\mathrm{EW};1\mathrm{loop}}=\frac{G_Fm_\mu ^2}{24\sqrt{2}\pi ^2}(5+(14s_W^2)^2).$$
(22)
Two-loop corrections to Eq. (22) are significant due to large logarithmic contributions ,

$$\frac{\alpha G_Fm_\mu ^2}{24\sqrt{2}\pi ^3}\left[41\frac{124}{3}s_W^2(12s_W^2)\right]\mathrm{ln}\frac{M_Z^2}{m_\mu ^2}.$$
(23)
The fermionic two-loop result, including some additional logarithms and sub-leading contributions, has been obtained in Ref.  for small and large values of $`M_H`$. We constructed an interpolation formula for other values of $`M_H`$ which reads in units of $`\alpha (m_\mu )G_Fm_\mu ^2/(3\sqrt{2}\pi ^3)`$,

$$\frac{13+6\mathrm{ln}x}{9}(1\omega )+x\left[3+\frac{\pi ^2}{3}+(\mathrm{ln}x+1)^2\right]\omega ,$$
(24)
where $`x=\widehat{m}_t^2/M_H^2`$ and $`\omega =e^{3.31849x}`$. Eq. (24) reproduces the exact result  for $`M_H=\widehat{m}_t`$. The bosonic two-loop corrections have been obtained only for large values of $`M_H`$ and as an expansion in $`s_W^2`$ . We take the leading non-logarithmic contribution in the $`M_H\mathrm{}`$ limit ($`a_2`$ in Ref. ) as the uncertainty induced by sub-leading bosonic two-loop effects. The leading logarithms to three-loop order have also been computed  and included in our analysis. We obtain $`a_\mu ^{\mathrm{EW}}=(1.52\pm 0.03)\times 10^9`$.

Summation of the various contributions gives $`(g_\mu 2\frac{\alpha }{\pi })/2=4506.28\times 10^9`$. Including our evaluation of $`a_\mu `$ into a global analysis of electroweak data yields,

$$\frac{1}{2}(g_\mu 2\frac{\alpha }{\pi })=(4506.35\pm 0.37\pm 0.51)\times 10^9,$$
(25)
where the first error includes $`\mathrm{\Delta }\alpha _{\mathrm{had}}`$ and all other parametric uncertainties, such as the $`\pm 0.028`$ error in $`\alpha _s(\mu _0)`$ obtained from current global fits. The extracted Higgs mass from this fit is $`M_H=88_{33}^{+49}`$ GeV. The current experimental world average, $`(4510.55\pm 1.51)\times 10^9`$ (which differs by $`2.6\sigma `$ from Eq. (25)) needs to be included into the analysis, as well. This is achieved by combining the experimental error with the second error in Eq. (25) in quadrature (for a total of $`1.59\times 10^9`$). The global fit then yields $`4506.52\pm 0.36`$ (parametric) with a pull of $`2.5`$. The crucial observation here is that now $`M_H=83_{31}^{+47}`$ GeV, i.e. both, the central value and the uncertainty of $`M_H`$ decrease. The increase in the precision in $`M_H`$ gained by properly correlating it to $`\mathrm{\Delta }\alpha _{\mathrm{had}}`$ (and assuming the SM) is almost identical to that provided by the $`M_W`$ measurement at the Tevatron Run I.

To summarize, we carefully examined the error correlations between $`a_\mu `$ and other quantities entering electroweak tests, especially $`\mathrm{\Delta }\alpha (M_Z)`$, $`\alpha _s`$, and the quark masses. We derived new analytical results for the hadronic contributions and showed that the proper treatment discussed in this article has a significant effects on the extraction of $`M_H`$, which could (depending on the future experimental findings) become even more dramatic.

Acknowledgements: It is a pleasure to thank P. Langacker and W. Marciano for valuable discussions. This work was supported in part by the US Department of Energy grant EY-76-02–3071, CNSF-10047005, and a fund for Trans-Century Talents.