text:
# 1 Introduction

## 1 Introduction

The concept of Operator Product Expansion (OPE) was applied to QCD sum rules in  to parametrize the nonperturbative effects. The operators of increasing dimension, constructed from quark and gluon fields, or condensates, constitute the OPE series, which is added to the perturbative ones. In case of heavy quark correlators the quark condensates are not essential and OPE series start from the dimension $`D=4`$ gluon condensate

$$\frac{\alpha _s}{\pi }G_{\mu \nu }^aG_{\mu \nu }^a$$

for which the authors of  have obtained the estimation $`0.012\mathrm{GeV}^4`$ from vector charmonium sum rules. In  the charmonium sum rules were studied in pseudoscalar channel and it was predicted the mass of the lowest $`\eta _c`$ state $`3.00\pm 0.03\mathrm{GeV}`$. This result was in contradiction with available to that time experimental information. Later measurements found the mass of $`\eta _c`$ close to $`3.0\mathrm{GeV}`$, which was considered as a triumph of QCD.

Since then various sum rules were analyzed in many publications111See  for the list of publications in order to obtain or specify the value of the gluon condensate. In the recent paper  the vector charmonium sum rules were reconsidered with $`\alpha _s^2`$-corrections of the perturbative series and $`\alpha _s`$-corrections to the condensate contribution, and up to date experimental data. The analysis of  resulted to the gluon condensate $`0.009\pm 0.007\mathrm{GeV}^4`$ and $`c`$-quark mass $`m_c(m_c)=1.275\pm 0.015\mathrm{GeV}`$. Despite high accuracy of experimental data and $`c`$-quark mass determination, the accuracy of the gluon condensate remains $`\pm 100\%`$, and zero value is not excluded.

It also seems reasonable to reanalyze the pseudoscalar sum rule, taking into account the information obtained in . Now the mass of $`\eta _c`$ is known experimentally with high accuracy $`2979.7\pm 1.5\mathrm{MeV}`$ , so we invert the problem and find a restriction on the charm quark mass and gluon condensate, imposed by this sum rule. A special attention should be paid to the correlation between these two values, since a variation of one parameter leads to the change of another.

The sum rule technique goes as follows. The correlator of the pseudoscalar charm currents is defined as:

$$\mathrm{\Pi }^p(q^2)=i𝑑xe^{iqx}TJ^p(x)J^p(0),J^p=2im\overline{c}\gamma _5c$$
(1)
We define the pseudoscalar current as $`J^p=_\mu J_\mu ^a`$, $`J_\mu ^a=\overline{c}\gamma _5\gamma _\mu c`$ is axial vector current. Within the narrow width approximation the imaginary part is:

$$\mathrm{Im}\mathrm{\Pi }^p(q^2+i0)=\pi \underset{\eta }{}\delta (q^2m_\eta ^2)\left|0|J^p(0)|\eta \right|^2$$
(2)
The sum goes over the pseudoscalar states222The next to lightest pseudoscalar state $`\eta _c(2S)`$ with mass $`3654\pm 6(stat)\pm 8(syst)\mathrm{MeV}`$ was discovered recently  with $`J^{PC}=0^+`$. The correlator (1) is quadratically divergent, so the dispersion relation requires double subtraction:

$$\mathrm{\Pi }^p(q^2)=c_0+q^2c_1+\frac{q^4}{\pi }_0^{\mathrm{}}\frac{\mathrm{Im}\mathrm{\Pi }^p(s+i0)}{s^2(sq^2)}𝑑s$$
(3)
provided the integral in the rhs is convergent, $`c_0`$, $`c_1`$ are unknown constants. In order to suppress the contribution of the higher states in (2) as well as continuum contribution, one considers the derivatives of the polarization operator in euclidean region $`Q^2q^2>0`$, the so-called moments:

$`M_n^p(Q^2)`$ $``$ $`{\displaystyle \frac{8\pi ^2}{3n!}}\left({\displaystyle \frac{d}{dQ^2}}\right)^n\mathrm{\Pi }^p(Q^2)={\displaystyle \frac{8\pi }{3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\mathrm{Im}\mathrm{\Pi }^p(s+i0)}{(s+Q^2)^{n+1}}}𝑑s`$ (4)
$`={\displaystyle \frac{8\pi ^2}{3}}{\displaystyle \underset{\eta }{}}{\displaystyle \frac{\left|0\left|J^p(0)\right|\eta \right|^2}{(m_\eta ^2+Q^2)^{n+1}}}`$
where $`n2`$. The matrix elements $`0\left|J^p(0)\right|\eta `$ are not known experimentally. But if one considers the ratio of some two moments at sufficiently high $`n`$, the contribution of the lightest state $`\eta _c(1S)`$ becomes dominant and

$$\frac{M_n^p(Q^2)}{M_{n+1}^p(Q^2)}=m_{\eta _c}^2+Q^2,n\mathrm{}$$
(5)
This property was exploited to predict the mass of $`\eta _c`$ in . An essential point was noticed in : the QCD corrections to the moments are large at $`Q^2=0`$, so the sum rules should be considered at $`Q^2>0`$. Moreover, huge contribution of the dimension 8 operators $`G^4`$ to the moments at $`Q^2=0`$  becomes tolerable at $`Q^24m_c^2`$.

The subject of this paper is a detailed analysis of the sum rule (5). In the next section the perturbative and OPE corrections to the correlator (1) are described. Section 2 is devoted to the moments both in the pole and $`\overline{\mathrm{MS}}`$ scheme for the charm quark mass. In the Section 3 various contribution to the pseudoscalar sum rule (5) are studied in details for typical values of the charm mass and gluon condensate. The higher dimension $`D=6,8`$ gluon operators are calculated both in the instanton and factorization model. In the final Section the restriction on the $`c`$-quark mass and the gluon condensate $`aG^2`$ are obtained.

## 2 Pseudoscalar correlator in QCD

In QCD the polarization function (1) consists of perturbative part and operator product expansion:

$$\mathrm{\Pi }^p=\mathrm{\Pi }_{\mathrm{PT}}^p+\mathrm{\Pi }_{\mathrm{OPE}}^p$$
(6)
The perturbative part is determined by its imaginary part via dispersion relation (3). The imaginary part is parametrized by the coefficient functions $`R^{(k),p}`$ in the expansion by the running QCD coupling $`a(\mu ^2)\alpha _s(\mu ^2)/\pi `$:

$$\mathrm{Im}\mathrm{\Pi }_{\mathrm{PT}}^p(s+i0)=\frac{3sm^2}{2\pi }\underset{k0}{}R^{(k),p}(s,\mu ^2)a^k(\mu ^2)$$
(7)
It is simpler to parametrize the functions $`R^{(k),p}`$ in terms of the pole masse $`m`$ of $`c`$-quark. The first two terms do not depend on the scale $`\mu ^2`$. They are known analytically :

$`R^{(0),p}`$ $`=`$ $`v`$
$`R^{(1),p}`$ $`=`$ $`{\displaystyle \frac{v}{2}}(7v^2)+4v\left(\mathrm{ln}{\displaystyle \frac{1v^2}{4}}{\displaystyle \frac{4}{3}}\mathrm{ln}v\right)+{\displaystyle \frac{19+2v^2+3v^4}{12}}\mathrm{ln}{\displaystyle \frac{1+v}{1v}}+{\displaystyle \frac{8}{3}}(1+v^2)`$ (8)
$`\times \left[\mathrm{\hspace{0.17em}2}\mathrm{Li}_2\left({\displaystyle \frac{1v}{1+v}}\right)+\mathrm{Li}_2\left({\displaystyle \frac{1v}{1+v}}\right)+\left({\displaystyle \frac{3}{2}}\mathrm{ln}{\displaystyle \frac{1+v}{2}}\mathrm{ln}v\right)\mathrm{ln}{\displaystyle \frac{1+v}{1v}}\right]`$
Here and below $`v=\sqrt{11/z}`$, $`z=s/(4m^2)`$. The function $`R^{(2),p}`$ is usually decomposed into the following gauge invariant parts:

$$R^{(2),p}=C_F^2R_A^{(2),p}+C_AC_FR_{NA}^{(2),p}+C_FTn_lR_l^{(2),p}+C_FTR_F^{(2),p}+C_FTR_S^{(2),p}$$
(9)
where $`C_A=3`$, $`C_F=4/3`$, $`T=1/2`$ are group constants and $`n_l=3`$ is the number of light quarks. The function $`R_l^{(2),p}`$ comes from the diagram with massless quark loop. It was found in  and in our normalization takes the form:

$$R_l^{(2),p}=\left(\frac{1}{4}\mathrm{ln}\frac{\mu ^2}{4s}\frac{5}{12}\right)R^{(1),p}+\delta _P^{(2)}$$
(10)
where the function $`\delta _P^{(2)}`$ is given by equation (110) in ref . The function $`R_F^{(2),p}`$ comes from the diagram with 2 massive quark loops. For $`s<16m^2`$ it contains only the contribution of virtual massive quarks and has the form :

$$R_F^{(2),p}=\mathrm{\hspace{0.17em}2}v\mathrm{Re}P_Q^{(2)}\frac{1}{4}R^{(1),p}\mathrm{ln}\frac{\mu ^2}{m^2},s<16m^2$$
(11)
where $`P_Q^{(2)}`$ is second order correction to the pseudoscalar current vertex from the diagram with massive quark loop; it is given by equation (169) in ref . For $`s>16m^2`$ the 4-particle cut must be included in $`R_F^{(2),p}`$. It is given by the double integral, eq. (97) in ref , which cannot be taken analytically. Here, however, the total function $`R_F^{(2),p}`$ can be replaced by its high-energy asymptotic, available to the terms $`m^8/s^4`$ in .

The functions $`R_A^{(2),p}`$ and $`R_{NA}^{(2),p}`$ correspond to the diagrams with single massive quark loop and various gluon exchanges. They are not known analytically, so one has to use some approximations. It turns out, that the moments, computed by the dispersion relation (4), are sensitive to the choice of these approximations. The accuracy of the moments becomes especially important in $`\overline{\mathrm{MS}}`$ scheme, where there is a sufficient cancellation between large terms (see eq. (26) below), so we describe this point in details.

The first 8 moments $`M_2`$$`M_9`$ at $`Q^2=0`$ are known analytically . We will require, that the approximations for $`R_A^{(2),p}`$ and $`R_{NA}^{(2),p}`$ must reproduce these moments with high accuracy being substituted into the dispersion integral in (4)333We found the approximations proposed in ref , eqs (39), (40) having insufficient accuracy to satisfy these requirement.. As usual, we shall apply the conformal mapping and Pade approximation for the relevant parts of the polarization function $`\mathrm{\Pi }^p`$ and take the imaginary part after then, see Appendix A for details. Although such approximations are constructed so that they reproduce low-$`q^2`$ expansion of the polarization function, they do not give exact values of the first 8 moments at $`Q^2=0`$, computed by taking the dispersion integral in (4). Indeed, the Pade approximations have extra poles away from the cut $`z=[1,\mathrm{})`$ and, strictly speaking, the dispersion relation (3) is not valid for them. The approximated formulas for $`R_A^{(2),p}`$ and $`R_{NA}^{(2),p}`$, used in this paper, are given in the equation (46) of Appendix A.

The last term in (9) $`R_S^{(2),p}`$ is the so-called singlet part with 2 triangle quark loops. This part contains the 2-gluon cut, which is proportional to the 2-photon decay width of the pseudoscalar boson. It is known analytically :

$$R_{gg}^{(2),p}=TC_F\frac{|f(z)|^2}{2z}$$
(12)
where

$$f(z)=\frac{1}{2}_0^1\frac{dx}{x}\mathrm{ln}\left[14zx(1x)\right]=\{\begin{array}{cc}\mathrm{arcsin}^2\sqrt{z}& ,z<1\\ \frac{1}{4}\left(\mathrm{ln}\frac{1+v}{1v}i\pi \right)^2& ,z>1\end{array}$$

The contribution of purely gluonic states to the heavy quark current correlators was discussed in details in  (3-gluon state in case of vector currents). In the narrow width approximation (2) only charmed states are taking into account. Since the 2-gluon state is not associated with any charmonium state, we subtract $`R_{gg}^{(2),p}`$ in the dispersion relation (3) and take the integral from $`s=4m^2`$, in accordance with suggestion of . The approximation for $`R_S^{(2),p}`$ without the 2-gluon cut is given in the eq (46) of Appendix A.

The leading in $`\alpha _s`$ order operator series $`\mathrm{\Pi }_{\mathrm{OPE}}^p`$ (6) for the heavy quark correlator has been computed in  up to operators of dimension $`D=8`$. This series can be compactly expressed in terms of Gauss hypergeometric functions:

$$\mathrm{\Pi }_{\mathrm{OPE}}^p(s)=\frac{1}{2\pi ^2}\underset{k2}{}\underset{j}{}\frac{O_k^{(j)}}{(4m^2)^{k2}}\underset{i=1}{\overset{2k1}{}}\frac{c_{k,i}^{(p,j)}}{k!(1/2)_{k1}}{}_{2}{}^{}F_{1}^{}\left(\genfrac{}{}{0pt}{}{1,k+i2}{k1/2}|\frac{s}{4m^2}\right)$$
(13)
where $`(1/2)_{k1}`$ is Pochhammer symbol, $`O_k^{(j)}`$ are the operators of dimension $`D=2k`$. For the heavy quark correlators there is single operator of dimension $`D=4`$

$$O_2=g^2G_{\mu \nu }^aG_{\mu \nu }^a$$
(14)
and 2 operators of dimension $`D=6`$:

$$O_3^{(1)}=g^3f^{abc}G_{\mu \nu }^aG_{\nu \lambda }^bG_{\lambda \mu }^c,O_3^{(2)}=g^4j_\mu ^aj_\mu ^a$$
(15)
where $`gj_\mu ^a=G_{\mu \nu ;\nu }^a=\frac{g}{2}_{q=u,d,s}\overline{q}\gamma _\mu \lambda ^aq`$. We choose 7 independent operators of dimension $`D=8`$ according to :

$`O_4^{(1)}`$ $`=`$ $`\left(g^2d^{abc}G_{\mu \nu }^bG_{\alpha \beta }^c\right)^2+{\displaystyle \frac{2}{3}}\left(g^2G_{\mu \nu }^aG_{\alpha \beta }^a\right)^2,O_4^{(2)}=\left(g^2f^{abc}G_{\mu \nu }^bG_{\alpha \beta }^c\right)^2,`$
$`O_4^{(3)}`$ $`=`$ $`\left(g^2d^{abc}G_{\mu \alpha }^bG_{\alpha \nu }^c\right)^2+{\displaystyle \frac{2}{3}}\left(g^2G_{\mu \alpha }^aG_{\alpha \nu }^a\right)^2,O_4^{(4)}=\left(g^2f^{abc}G_{\mu \alpha }^bG_{\alpha \nu }^c\right)^2,`$
$`O_4^{(5)}`$ $`=`$ $`g^5f^{abc}G_{\mu \nu }^aj_\mu ^bj_\nu ^c,O_4^{(6)}=g^3f^{abc}G_{\mu \nu }^aG_{\nu \lambda }^bG_{\lambda \mu ;\alpha \alpha }^c,O_4^{(7)}=g^4j_\mu ^aj_{\mu ;\alpha \alpha }^a.`$ (16)
The coefficients $`c_{k,i}^{(p,j)}`$ in (13) can be obtained from :

$$c_{2,i}^{(p)}=(0,\mathrm{\hspace{0.17em}1}/4,1/12),c_{3,i}^{(p,j)}=\begin{array}{cccccc}& & & & & \\ j\backslash i& 1& 2& 3& 4& 5\\ & & & & & \\ 1& 0& 3/2& 5& 22/5& 1\\ 2& 1/3& 8/3& 3& 28/15& 4/3\end{array},$$

$$c_{4,i}^{(p,j)}=\begin{array}{cccccccc}& & & & & & & \\ j\backslash i& 1& 2& 3& 4& 5& 6& 7\\ & & & & & & & \\ 1& 1/3& 31& 258& 1580/3& 365& 70& 0\\ 2& 1/6& 10& 111& 946/3& 4873/14& 144& 10\\ 3& 4/3& 16& 168& 596/3& 106& 140& 0\\ 4& 1/3& 122& 946& 5480/3& 7435/7& 20& 92\\ 5& 10/3& 12& 260& 2888/3& 7290/7& 304& 24\\ 6& 0& 22& 152& 174& 1108/7& 296& 96\\ 7& 2& 36& 100& 128& 1314/7& 208& 72\end{array}.$$
(17)
The $`\alpha _s`$ correction to the $`D=4`$ condensate contribution was obtained analytically in . One could differentiate it $`n`$ times to obtain the moments. However, we prefer to use a dispersion-like relation for this correction, constructed in Appendix B, which is convenient for numerical calculation of the moments, especially for high $`n`$.

## 3 Moments in $`\overline{\mathrm{𝐌𝐒}}`$ scheme

At first let us consider the moments in the pole-mass scheme. In QCD the moments (4) are expanded by the running QCD coupling $`\alpha _s(\mu ^2)`$. The approximation, used in this paper, includes the following ingredients: 1) the perturbative series up to $`\alpha _s^2`$ order, 2) the operator series up to dimension $`D=8`$, and 3) $`\alpha _s`$ correction to the $`G^2`$ operator contribution. Adding all pieces together, we write down the following expression for the pseudoscalar moments:

$$M_n^p(Q^2)=\underset{k=0}{\overset{2}{}}M_n^{(k),p}(Q^2)a^k(m^2)+\underset{k=2}{\overset{4}{}}\underset{j}{}O_k^{(j)}M_n^{(Ok,j),p}(Q^2)+O_2M_n^{(O2)(1),p}(Q^2)a(m^2),$$
(18)
for definiteness the coupling $`a\alpha _s/\pi `$ is taken at the scale $`\mu ^2=m^2`$. As discussed in previous section, the perturbative moments are taken without 2-gluon cut (12):

$$M_n^{(k),p}(Q^2)=\mathrm{\hspace{0.17em}4}m^2_{4m^2}^{\mathrm{}}\frac{s(R^{(k),p}R_{gg}^{(k),p})(s,m^2)}{(s+Q^2)^{n+1}}𝑑s$$
(19)
The leading order can be expressed in terms of Gauss hypergeometric function:

$$M_n^{(0),p}(Q^2)=\frac{1}{(4m^2)^{n2}}\frac{\mathrm{\Gamma }(n1)}{2(1/2)_n}{}_{2}{}^{}F_{1}^{}\left(\genfrac{}{}{0pt}{}{n1,n+1}{n+1/2}|\frac{Q^2}{4m^2}\right)$$
(20)
where $`(a)_n\mathrm{\Gamma }(a+n)/\mathrm{\Gamma }(a)`$ is Pochhammer symbol. The higher order perturbative moments are computed numerically by (19).

The contribution of the operators $`O_k^{(j)}`$ to the moments can be easily obtained by differentiating eq (13):

$$M_n^{(Ok,j),p}(Q^2)=\frac{4}{3}\underset{i=1}{\overset{2k1}{}}\frac{c_{k,i}^{(p,j)}}{(4m^2)^{n+k2}}\frac{(k+i2)_n}{k!(1/2)_{n+k1}}{}_{2}{}^{}F_{1}^{}\left(\genfrac{}{}{0pt}{}{n+1,n+k+i2}{n+k1/2}|\frac{Q^2}{4m^2}\right)$$
(21)
The $`\alpha _s`$-correction to the $`D=4`$ gluon condensate contribution can be obtained by differentiating eq (53) of Appendix B:

$$M_n^{(O2)(1),p}(Q^2)=\frac{2}{3(4m^2)^n}\underset{i=1}{\overset{3}{}}\frac{(n+1)_{i1}}{(i1)!}\left[\frac{\pi ^2f_i^p}{(1+y)^{n+i}}+_1^{\mathrm{}}\frac{F_i^p(z)}{(z+y)^{n+i}}𝑑z\right]$$
(22)
where $`y=Q^2/(4m^2)`$, the constants $`f_i^p`$ and the functions $`F_i^p(z)`$ are given in eqs (50) and (52) of Appendix B. Notice, that eqs (19)–(22) are applicable for noninteger $`n`$ also.

Similarly to the vector case , the $`\alpha _s`$-corrections to the moments are unacceptably large in the pole mass scheme and the series (18) is divergent. The pole mass, in fact, is the mass of free quark. Since the quarks exist only in form of strongly bounded states, the physical meaning of the pole quark mass is rather unclear; it cannot be found from the sum rules with a good accuracy.

Instead of the pole mass one introduces another effective mass parameter, to improve the convergence of the perturbative series. Authors of  used the mass, renormalized at the euclidean point $`p^2=m^2`$. In this paper we shall use the most popular choice for today: the gauge invariant mass in the modified minimal subtraction ($`\overline{\mathrm{MS}}`$) scheme taken at the scale, equal to the mass itself $`\overline{m}\overline{m}(\overline{m}^2)`$. The pole mass $`m`$ is perturbatively expressed in terms of $`\overline{m}`$:

$$\frac{m^2}{\overline{m}^2}=\mathrm{\hspace{0.17em}1}+\underset{n1}{}K_na^n(\overline{m}^2)$$
(23)
The 2-loop factor was found, in particular, in  while the 3-loop factor was recently calculated in :

$`K_1`$ $`=`$ $`{\displaystyle \frac{8}{3}}`$
$`K_2`$ $`=`$ $`28.66462.0828n_l=\mathrm{\hspace{0.17em}22.4162}`$
$`K_3`$ $`=`$ $`417.03956.0871n_l+1.3054n_l^2=\mathrm{\hspace{0.17em}260.526}`$ (24)
We put $`n_l=3`$ in the last column. Now we reexpand the moments (18) by the QCD coupling $`a(\overline{m}^2)`$:

$$M_n^p(Q^2)=\underset{k=0}{\overset{2}{}}\overline{M}_n^{(k),p}(Q^2)a^k(\overline{m}^2)+\underset{k=2}{\overset{4}{}}\underset{j}{}O_k^{(j)}\overline{M}_n^{(Ok,j),p}(Q^2)+O_2\overline{M}_n^{(O2)(1),p}(Q^2)a(\overline{m}^2),$$
(25)
where $`\overline{M}_n^{(k),p}(Q^2)`$ are expressed in terms of $`M_n^{(k),p}(Q^2)`$ :

$`\overline{M}_n^{(0),p}(Q^2)`$ $`=`$ $`M_n^{(0),p}`$
$`\overline{M}_n^{(1),p}(Q^2)`$ $`=`$ $`M_n^{(1),p}K_1(nd/2)M_n^{(0),p}+K_1(n+1)Q^2M_{n+1}^{(0),p}`$
$`\overline{M}_n^{(2),p}(Q^2)`$ $`=`$ $`M_n^{(2),p}K_1(nd/2)M_n^{(1),p}+K_1(n+1)Q^2M_{n+1}^{(1),p}`$
$`+(nd/2)\left[{\displaystyle \frac{K_1^2}{2}}(n+1d/2)K_2\right]M_n^{(0),p}`$
$`+(n+1)\left[K_2K_1^2(n+1d/2)\right]Q^2M_{n+1}^{(0),p}`$
$`+{\displaystyle \frac{K_1^2}{2}}(n+1)(n+2)Q^4M_{n+2}^{(0),p}`$
$`\overline{M}_n^{(Ok,j),p}(Q^2)`$ $`=`$ $`M_n^{(Ok,j),p}`$
$`\overline{M}_n^{(O2)(1),p}(Q^2)`$ $`=`$ $`M_n^{(O2),p}K_1(n+2d/2)M_n^{(G,0),p}+K_1(n+1)Q^2M_{n+1}^{(G,0),p}`$ (26)
where $`d=4`$ is the dimension of the pseudoscalar function $`\mathrm{\Pi }^p(Q^2)`$, all $`M_n^{(i),p}`$ in the rhs are computed with $`\overline{\mathrm{MS}}`$ mass $`\overline{m}`$. The series (25) is much better convergent than (18). The numerical values of the ratios $`\overline{M}^{(1,2)}/\overline{M}^{(0)}`$ and $`\overline{M}^{(O2)(1)}/\overline{M}^{(O2)}`$ for $`Q^2/(4\overline{m}^2)=0,1,2`$ and $`n=230`$ are given in the Table 1 of Appendix C. Notice, that the values of $`\overline{M}^{(2)}`$ are approximate; other approximations for $`R^{(2)}`$ may lead to the moments $`\overline{M}^{(2)}`$, which differ from the numbers of the Table 1 within $`510\%`$.

The expansion (25) goes by $`a(\overline{m}^2)`$. If one takes the QCD coupling at some another scale $`\mu ^2`$, the function $`\overline{M}^{(2),p}`$ changes:

$$a(\overline{m}^2)a(\mu ^2),\overline{M}_n^{(2),p}(Q^2)\overline{M}_n^{(2),p}(Q^2)+\overline{M}_n^{(1),p}(Q^2)\beta _0\mathrm{ln}\frac{\mu ^2}{\overline{m}^2}$$
(27)
so that the series (25) is $`\mu ^2`$-independent at the order $`\alpha _s^2`$.

## 4 Pseudoscalar sum rule

It is convenient to define a dimensionless ratio of the pseudoscalar moments:

$$r_n(Q^2)=\frac{1}{4\overline{m}^2}\frac{M_n^p(Q^2)}{M_{n+1}^p(Q^2)}\frac{m_{\eta _c}^2+Q^2}{4\overline{m}^2},n\mathrm{}$$
(28)
Theoretical ratio depends on the quark mass $`\overline{m}`$, QCD coupling $`\alpha _s`$ and condensates. But if the dimensionless parameters $`Q^2/(4\overline{m}^2)`$, $`aG^2/(4\overline{m}^2)^2`$ etc. are fixed, the l.h.s. of (28) does not depend on the quark mass $`\overline{m}`$ (in fact, the QCD coupling depends on the scale, which itself may depend on $`\overline{m}`$; but this dependence is weak within the range of error of $`\overline{m}`$). So one may use the ratio (28) to find the $`\overline{\mathrm{MS}}`$ charm quark mass $`\overline{m}`$ for given condensates and QCD coupling.

The QCD coupling constant $`\alpha _s`$ is universal value and can be taken from other experiments. As input parameter, it is convenient to take $`\alpha _s`$ at the $`\tau `$-lepton mass :

$$\alpha _s(m_\tau ^2)=\mathrm{\hspace{0.17em}0.33}\pm 0.03,m_\tau =\mathrm{\hspace{0.17em}1.777}\mathrm{GeV}$$
(29)
Using this value as the boundary condition in the renormalization group equation, the QCD coupling can be evaluated at any scale. As argued in , the most natural scale for $`\alpha _s`$ is

$$\mu ^2=Q^2+\overline{m}^2$$
(30)
Indeed, in the limit $`Q^2\overline{m}^2`$ we come to natural massless choice $`\alpha _s(Q^2)`$, while at $`Q^2=0`$ it becomes $`\alpha _s(\overline{m}^2)`$. Later we shall vary the scale (30) to check the stability of results.

The $`\overline{\mathrm{MS}}`$ charm quark mass is determined from vector charmonium sum rules with high accuracy. The analysis of the moments at $`Q^2=0`$ with $`\alpha _s^2`$ corrections leads to $`\overline{m}=1.304\pm 0.027\mathrm{GeV}`$ in  and $`1.23\pm 0.09\mathrm{GeV}`$ in . The authors of  neglected the condensate contribution, while in  the value $`aG^2=0.024\pm 0.012\mathrm{GeV}^4`$ was employed. In fact, the gluon condensate weakly affects on the mass value. But for the condensate determination the mass accuracy is especially important: a small mass variation leads to significant condensate change. As noticed in , the perturbative $`\alpha _s`$ and $`\alpha _s^2`$ corrections to the vector moments $`M_n(Q^2)`$ in $`\overline{\mathrm{MS}}`$ scheme are strongly suppressed for $`Q^2/(4\overline{m}^2)n/51`$ and $`n>5`$. The analysis of  at $`Q^2/(4\overline{m}^2)=1,2`$ allowed to determine the $`c`$-quark mass with high acuracy:

$$\overline{m}(\overline{m})=\mathrm{\hspace{0.17em}1.275}\pm 0.015\mathrm{GeV}$$
(31)
independently of the condensate value. (If the condensate is fixed, the error in (31) can be reduced even further.) This result is close to the recent lattice calculation $`\overline{m}=1.26\pm 0.04(stat)\pm 0.12(syst)\mathrm{GeV}`$ .

For the mass (31) one gets the ratio (28)

$$r_n(Q^2)=\frac{Q^2}{4\overline{m}^2}+\mathrm{\hspace{0.17em}1.37}\pm 0.03$$
(32)
in the limit $`n\mathrm{}`$. Which values of $`Q^2/(4\overline{m}^2)`$ are convenient for the sum rule analysis? The choice $`Q^2=0`$ is not appropriate, since the perturbative corrections to the moments are large for almost all $`n`$, even in $`\overline{\mathrm{MS}}`$ scheme. Large $`Q^212\overline{m}^2`$ are also dangerous: in particular, when one changes the scale of $`\alpha _s`$ in (27), the effective expansion parameter $`a\beta _0\mathrm{ln}(Q^2/\overline{m}^2)`$ becomes large $`0.5`$. In what follows we shall use two choices $`Q^2/(4\overline{m}^2)=1,2`$.

The theoretical ratios $`r_n(4\overline{m}^2)`$ and $`r_n(8\overline{m}^2)`$ are plotted versus $`n`$ in the Fig 1a) and 1b) respectively. The lower shaded curve is purely perturbative, i.e. for $`aG^2=0`$. The central line of the shaded area corresponds to the central value of $`\alpha _s`$ (29), the errorband covers the error of the coupling $`\alpha _s`$ in (29). One sees, that the agreement with (32) is achieved within relatively narrow range of $`n`$: $`n16`$ for $`Q^2=4\overline{m}^2`$ and $`n24`$ for $`Q^2=8\overline{m}^2`$. If we look at the Table 1, the perturbative corrections to the moments in $`\overline{\mathrm{MS}}`$ scheme, as well as $`\alpha _s`$ correction to the condensate contribution, are minimal here. For higher $`n`$ these corrections grow rapidly and the perturbation theory cannot be trusted here. For lower $`n`$ the perturbative corrections are also large, and the leading order of the $`D=4`$ condensate contribution crosses 0 at some point, so the behavior of the $`\alpha _s`$-series is rather unclear here. Moreover, unknown contribution of $`\eta _c(2S)`$ and higher states to the experimental moments could be significant for low $`n`$.

Now we consider nonzero $`D=4`$ condensate. As an illustration, let us fix the ratio $`aG^2/(4\overline{m}^2)^2`$ $`=2\times 10^4`$, which corresponds to $`aG^20.008\mathrm{GeV}^4`$, close to the central value obtained in . The ratio $`r_n`$ with this condensate is shown by the upper shaded curves in Fig 1. The ratio becomes higher for nonzero condensate, which tells in favor of lower mass of $`c`$-quark. At $`Q^2=4\overline{m}^2`$ the ratio is even higher, than (32) for all $`n`$.

Several models were employed to estimate the higher dimension $`D=6,8`$ operators. Here we consider the dilute instanton gas model  and vacuum dominance (or factorization) model.

Instanton model. The vacuum configuration is codsidered as a dilute gas of noniteracting instantons with effective radis $`\rho _c`$ and concentration $`n_0`$. The radius $`\rho _c`$ varies within $`0.31\mathrm{fm}`$ in the literature. Here we shall use the estimation $`\rho _c=0.5\mathrm{fm}`$ obtained in . The instanton concentration $`n_0`$ is fitted to the $`D=4`$ gluon condensate $`aG^2=32n_0`$. Then one obtains the following expressions for the $`D=6,8`$ gluon condensates (15), (16):

$$O_3^{(1)}=\frac{12}{5\rho _c^2}O_2,O_3^{(2)}=\mathrm{\hspace{0.17em}0},$$

$$(O_4^{(1)},\mathrm{},O_4^{(7)})=(4,8,3,4,0,8,0)\frac{16}{7\rho _c^4}O_2$$
(33)
The ratio $`r_n`$ with the operators $`O_{3,4}`$ computed by the instanton model (33) is show by hatched curves in Fig 1. Upper hatched curve includes $`G^3`$ operator only, lower hatched curve includes both $`G^3`$ and $`G^4`$.

The $`G^3`$ operator contribution is small. But the contribution of the $`G^4`$ operators is large in the region of interest. Obviously the place, where the lower hatched curve crosses the perturbative one (lower shaded), the sum rule (28) with the operators (33) is not applicable. Here the $`G^4`$ contribution exceeds the leading order $`G^2`$, and the OPE series diverges.

It is a demonstration, that the higher order operators are essentialy overestimated in the instanton model . Moreover, their values strongly depend on the instanton size $`\rho _c`$, which is not strictly fixed. For this reason we finish the analysis in the instanton model. The main outcome of this analysis is relatively small contribution of the operator $`G^3`$, which will be ignored in what follows.

Factorization hypothesis. In the factorization model the $`G^4`$ operators are proportional to $`\left(G^2\right)^2`$. The operators with the light quark current $`j_\mu ^a`$ in (15), (16) can also be estimated by the factorization, but their size is much smaller. The operator $`O_4^{(6)}`$ with derivatives was taken as $`O_4^{(6)}M^2O_3^{(1)}`$ in , where $`M^20.3\mathrm{GeV}^2`$ characterizes the gluon virtuality in the vacuum. Alternatively, one may express this operator as

$$O_4^{(6)}=\mathrm{\hspace{0.17em}2}g^4f^{abc}G_{\mu \nu }^aG_{\nu \lambda }^bj_{\lambda ;\mu }^c+\mathrm{\hspace{0.17em}2}O_4^{(4)}$$

Since we neglect the operators with $`j_\mu ^a`$, we take $`O_4^{(6)}=2O_4^{(4)}`$ here (both estimations agree in the order of magnitude for typical condensates). Summarizing, we write down the $`D=8`$ operators as:

$$(O_4^{(1)},\mathrm{},O_4^{(7)})=(\frac{65}{144},\frac{5}{16},\frac{19}{72},\frac{1}{16},\mathrm{\hspace{0.17em}0},\frac{1}{8},\mathrm{\hspace{0.17em}0})\left(O_2\right)^2$$
(34)
The accuracy of the factorization is expected to be $`1/N_c^2`$, $`N_c=3`$ is the color number. (The $`1/N_c^2`$ ambiguity of the $`D=8`$ quark-gluon condensate factorization was explicitly demonstrated in .) Another version of the factorization, which employs the heavy quark expansion, was proposed in .

The ratio $`r_n`$ with the operators (34) is shown by the hatched curves in the Fig 2 for the $`D=4`$ gluon condensate $`G^2/(4\overline{m}^2)^2=2\times 10^4`$. Comparing the Figures 1 and 2 one sees, that the contribution of the $`G^4`$ operators in the factorization model is smaller than in the instanton one. It allows to establish certain stability region, where the ratio $`r_n`$ remains almost unchanged. This region is clearly visible for $`Q^2=8\overline{m}^2`$: at $`n=2026`$ the ratio is $`r_n=3.40\pm 0.01`$, which corresponds to the $`c`$-quark mass $`\overline{m}=1.260\pm 0.005\mathrm{GeV}`$. This mass is computed for the condensate $`aG^2=0.008\mathrm{GeV}^4`$.

In the same way the mass can be computed for other values of the condensate. A restriction on the charm quark mass for different condensates is calculated in the next section.

## 5 Restrictions on the $`𝒄`$-quark mass and $`𝑫\mathbf{=}\mathrm{𝟒}`$ gluon condensate

As the main result of the pseudoscalar charmonium sum rule (5), we may establish certain restrictions for the $`c`$-quark mass $`\overline{m}`$ for a given condensate $`aG^2`$.

At first, let us neglect the higher dimension operators $`G^3`$ and $`G^4`$. The calculation goes as follows. For a given $`Q^2/(4\overline{m}^2)`$ one should establish the range of $`n`$, where the perturbation theory as well as operator expansion can be trusted. It is reasonable to require, that the perturbative corrections may not exceed $`3040\%`$ of the leading term. The most dangerous is the $`\alpha _s`$-correction to the gluon condensate contribution $`\overline{M}_n^{(G,1)}`$. Keeping in mind typical size of the QCD couping $`\alpha _s/\pi 0.1`$, let us impose the restriction $`|\overline{M}_n^{(G,1)}/\overline{M}_n^{(G,0)}|<4`$. From the Table 1 we find the following range of $`n`$:

$$n=1419\mathrm{for}\frac{Q^2}{4\overline{m}^2}=1\mathrm{and}n=2230\mathrm{for}\frac{Q^2}{4\overline{m}^2}=2.$$
(35)
The perturbative corrections to the moments are also tolerable in this region: the first correction $`|\overline{M}_n^{(1)}/\overline{M}_n^{(0)}|<2.5`$ and the NNLO correction $`|\overline{M}_n^{(2)}/\overline{M}_n^{(0)}|<17`$. Then, we take some value of $`aG^2/(4\overline{m})^2`$ and find the maximal and minimal value of the ratio $`r_n(Q^2)`$ within this range of $`n`$. From these numbers we find the minimal and maximal values of the charm quark mass $`\overline{m}`$.

The results are shown by the shaded regions in the Fig 3. Fig 3a) and 3b) display the restrictions, obtained from the sum rule (28) at $`Q^2=4\overline{m}^2`$ and $`Q^2=8\overline{m}^2`$ respectively. Since unknown higher order in $`\alpha _s`$ moments are discarded everywhere, the results depend on the choice of the scale, at which $`\alpha _s`$ is taken. The dark area shows the acceptable region for the scale (30). The dashed and dotted lines display the boundaries of the acceptable region, if $`\overline{m}^2`$ is added to or subtracted from this scale. The scale dependence is weaker at $`Q^2=8\overline{m}^2`$.

It is clear from the Fig 3, that the pseudoscalar sum rule prefers lower values of the gluon condensate. In particular, for the mass $`\overline{m}=1.275\pm 0.015\mathrm{GeV}`$  one obtains the upper condensate limit $`aG^2<0.008\mathrm{GeV}^4`$.

Now let us include the higher dimension operators. As follows from the instanton model analisys, the contribution of the $`D=6`$ condensate $`G^3`$ is small in the region of interest. But the $`D=8`$ operators change the ratio $`r_n`$ essentially. At some $`n`$ their contribution exceeds the leading condensate $`G^2`$; at this point the OPE series is divergent. Let us require that the contribution of the $`D=8`$ operators $`O_4^{(i)}`$ to the moments may not exceed $`30\%`$ of the $`D=4`$ condensate contibution. This requirement further reduces the region of $`n`$ (35) depending on the condensate size. In the factorization model $`G^4\left(G^2\right)^2`$, and the region becomes smaller for higher condensate $`G^2`$.

The hatched regions in the Fig 3 display the inclusion of the $`G^4`$ operators in the factorization model. For $`aG^2>0.005\mathrm{GeV}^4`$ the $`D=8`$ operators change the ratio $`r_n`$ drastically. They compensate the leading condensate and the ratio $`r_n`$ becomes almost independent of the $`aG^2`$ condensate in the stability region. From the hatched area in the Fig 3 one gets the following limits of the $`c`$-quark mass:

$$\overline{m}_c(\overline{m}_c)=\mathrm{\hspace{0.17em}1.26}\pm 0.02\mathrm{GeV}$$
(36)
indenpendently on the condensate value. The mass (36) is in agreement with the result $`\overline{m}=1.275\pm 0.015\mathrm{GeV}`$, obtained from the vector charmonium sum rules in .

If the $`G^4`$ operators are included, it becomes rather difficult to obtain certain restrictions on the condensate value. However, for large condensate the stability region is narrow, and the results become unreliable. In particular, for $`aG^2>0.015\mathrm{GeV}^4`$ and $`Q^2=4\overline{m}^2`$ there is no region of $`n`$, where the OPE series looks convergent. This sets the natural limit of the condensate value, at which the pseudoscalar sum rule works.

## Acknowledgement

Author thanks B.L.Ioffe for discussions. The research described in this publication was made possible in part by Award No RP2-2247 of U.S. Civilian Research and Development Foundation for Independent States of Former Soviet Union (CRDF), by the Russian Found of Basic Research grant 00-02-17808 and INTAS grant 00-00334.

## Appendix A: Approximations for $`𝑹^{\mathbf{(}\mathrm{𝟐}\mathbf{)}\mathbf{,}𝒑}`$

Let us define the dimensionless coefficient functions for the perturbative correlator (6) as follows:

$$\mathrm{\Pi }^p(q^2)=4m^2q^2\underset{k0}{}\mathrm{\Pi }^{(k),p}(q^2)a(m^2),\mathrm{\Pi }^{(k),p}(q^2)=\frac{3q^2}{8\pi ^2}_0^{\mathrm{}}\frac{R^{(k),p}(s,m^2)}{s(sq^2)}𝑑s,$$
(37)
for definiteness we take the QCD coupling at the scale $`\mu ^2=m^2`$ and put the constants $`c_0=c_1=0`$ in the dispersion relation (3). The 3-loop function $`\mathrm{\Pi }^{(2),p}`$ is decomposed into 5 gauge invariant parts in the same way as $`R^{(2),p}`$ (9).

At first we consider the nonabelian part $`\mathrm{\Pi }_{NA}^{(2),p}`$. Its expansion near $`zq^2/(4m^2)=0`$ until $`z^8`$ is available in . Then, as usual, we reexpand this series in terms of the variable $`\omega `$, which naturally appears in the perturbative calculations:

$$\omega =\frac{1\sqrt{1z}}{1+\sqrt{1z}}$$
(38)
The expansion of the polarization operator in $`\omega `$ has appropriate analytical properties, namely the cut at $`z=[1,\mathrm{})`$. In many cases the Pade approximation was proved to have better accuracy, than Tailor series. The best results (see the discussion in Section 2) were obtained for the Pade approximation \[5/2\]:

$$\mathrm{\Pi }_{NA}^{(2),p}(\omega )=\frac{3}{16\pi ^2}\times $$

$$\times \frac{7.43220\omega +73.5001\omega ^2+2.69248\omega ^313.3868\omega ^40.91032\omega ^51.61120\omega ^6}{10.18579\omega 0.63849\omega ^2}$$
(39)
The accuracy of the Pade approximated abelian part $`\mathrm{\Pi }_A^{(2),p}`$ is worse because of Coulomb behavior $`1/\sqrt{1z}`$ near the threshold. It turns out, however, that the expansion in $`\omega `$ converges faster, if the multiplier $`1/(1\omega ^2)`$ is separated out:

$`\mathrm{\Pi }_A^{(2),p}(\omega )`$ $`=`$ $`{\displaystyle \frac{3}{16\pi ^2}}[10.8487\omega +90.9348\omega ^2+111.510\omega ^3+73.3467\omega ^4`$
$`+113.363\omega ^5+74.8959\omega ^6+114.366\omega ^7+77.2105\omega ^8+O(\omega ^9)]`$
$`=`$ $`{\displaystyle \frac{3}{16\pi ^2}}{\displaystyle \frac{1}{1\omega ^2}}[10.8487\omega +90.9348\omega ^2+100.662\omega ^317.5881\omega ^4`$
$`+1.85289\omega ^5+1.54919\omega ^6+1.00228\omega ^7+2.31461\omega ^8+O(\omega ^9)]`$
Now we construct the Pade approximation, which well reproduces all asymptotic and first 8 moments from the dispersion relation:

$$\mathrm{\Pi }_A^{(2),p}(\omega )=\frac{3}{16\pi ^2}\times $$

$$\times \frac{10.8487\omega +145.611\omega ^2+507.376\omega ^3+57.361\omega ^4565.406\omega ^5+94.514\omega ^6}{(1\omega ^2)(1+5.03984\omega 4.75471\omega ^2)}$$
(40)
Few more work should be done to construct the approximation for the singlet polarization function $`\mathrm{\Pi }_S^{(2),p}`$. Its expansion near $`z=0`$ until $`z^8`$ is available in . The singlet correlator contains intermediate massless 2-gluon state, so the cut starts from $`z=0`$, the expansion in  has the terms $`\mathrm{ln}(z)`$ and the conformal mapping procedure ($`z\omega `$) is not applicable here. As discussed in Section 2, in our sum rules we use the polarization operator without 2-gluon cut (19), so the correspondent part of the polarization operator should be subtracted from the result of :

$$\mathrm{\Pi }_{gg}^{(2),p}=\frac{3}{16\pi ^2}C_FT(y)_0^{\mathrm{}}\frac{|f(z)|^2}{z^2(z+y)}𝑑z,y=\frac{Q^2}{4m^2}$$
(41)
where the function $`f(z)`$ is given in (12). The integral from $`z=1`$ to $`\mathrm{}`$ is regular at $`y=0`$ and can be expanded by $`y`$ in Tailor series. But the integral from $`z=0`$ to 1 requires special care, since it behaves as $`\mathrm{ln}y`$ at $`y0`$. In order to obtain the expansion for small $`y`$, we suggest to use the following series for the function $`|f(z)|^2`$ for $`0<z<1`$ :

$$\mathrm{arcsin}^4\sqrt{z}=\frac{3}{2}\underset{n=2}{\overset{\mathrm{}}{}}z^n\frac{(n1)!}{n(1/2)_n}\underset{k=1}{\overset{n1}{}}\frac{1}{k^2}$$
(42)
Then we obtain the following expansion:

$$_0^1\frac{\mathrm{arcsin}^4\sqrt{z}}{z^2(z+y)}𝑑z=\frac{\mathrm{arcsh}^4\sqrt{y}}{y^2}\mathrm{ln}y+\underset{n=0}{\overset{\mathrm{}}{}}(y)^nI_n,$$
(43)
where the constants

$$I_n=\frac{3}{2}\underset{k=2,n+2}{\overset{\mathrm{}}{}}\frac{(k1)!}{k(kn2)(1/2)_k}\underset{j=1}{\overset{k1}{}}\frac{1}{j^2}$$

can be computed analytically for any $`n`$ with the help of recursive relations:

$`I_n`$ $`=`$ $`{\displaystyle \frac{\pi ^4}{16(n+2)}}+{\displaystyle \frac{2(n+1)!}{(n+2)(5/2)_n}}\{\pi ^2\mathrm{ln}2{\displaystyle \frac{7}{2}}\zeta _3+{\displaystyle \underset{k=1}{\overset{n+1}{}}}[{\displaystyle \frac{1}{k^2(n+2)}}+{\displaystyle \frac{2}{k^3}}`$ (44)
$`{\displaystyle \frac{(1/2)_k}{k^2(k1)!}}{\displaystyle \frac{\pi ^28}{2}}+{\displaystyle \frac{1}{k^2}}{\displaystyle \underset{j=1}{\overset{n+1}{}}}({\displaystyle \frac{1}{j+1/2}}{\displaystyle \frac{1}{j}}){\displaystyle \frac{1}{k^2}}{\displaystyle \underset{j=1}{\overset{k1}{}}}{\displaystyle \frac{(3/2)_j}{j^2(j1)!}}]\}`$
Now one obtains regular at $`y=0`$ Tailor expansion of the polarization operator $`\mathrm{\Pi }_S^{(2),p}`$ without the 2-gluon cut, applies the conformal mapping and constructs the Pade approximation:

$$\mathrm{\Pi }_{Sgg}^{(2),p}(\omega )\mathrm{\Pi }_S^{(2),p}\frac{1}{C_FT}\mathrm{\Pi }_{gg}^{(2),p}=\frac{3}{16\pi ^2}\times $$

$$\times \frac{7.86155\omega +6.98952\omega ^2+3.24217\omega ^31.68013\omega ^40.31547\omega ^50.00464\omega ^6}{10.05236\omega 0.32689\omega ^2}$$
(45)
Eventually we take the imaginary part of (39), (40), (45) and obtain the correspondent coefficient functions $`R^{(2),p}`$:

$`R_{NA}^{(2),p}`$ $`=`$ $`{\displaystyle \frac{8\pi }{3}}\mathrm{Im}\mathrm{\Pi }_{NA}^{(2),p}(\omega )+{\displaystyle \frac{11}{16}}R^{(1),p}\mathrm{ln}{\displaystyle \frac{\mu ^2}{m^2}},`$
$`R_A^{(2),p}`$ $`=`$ $`{\displaystyle \frac{8\pi }{3}}\mathrm{Im}\mathrm{\Pi }_A^{(2),p}(\omega ),`$
$`R_S^{(2),p}{\displaystyle \frac{1}{C_FT}}R_{gg}^{(2),p}`$ $`=`$ $`{\displaystyle \frac{8\pi }{3}}\mathrm{Im}\mathrm{\Pi }_{Sgg}^{(2),p}(\omega ),\mathrm{where}\omega ={\displaystyle \frac{1+i\sqrt{z1}}{1i\sqrt{z1}}}`$ (46)
## Appendix B: $`𝜶_𝒔`$-correction to the condensate contribution

The $`\alpha _s`$ correction to the $`D=4`$ gluon condensate contribution was found in . Let us parametrize it by dimensionless function $`f^{(1),p}(z)`$:

$$\mathrm{\Pi }_{\mathrm{OPE}}^p(s)=\mathrm{}+\frac{O_2}{4\pi ^2}af^{(1),p}(z),$$
(47)
where dots denote the leading order operator contribution (13). Here we construct a dispersion-like relation for this function, convenient for numerical calculation of the moments. We will follow the method, used in  for the vector current correlator.

The imaginary part is:

$$\mathrm{Im}f^{(1),p}(z+i0)=\frac{\pi }{96z^3v^5}\left[P_2^P(z)+\frac{P_3^P(z)}{zv}\mathrm{ln}\frac{1v}{1+v}+P_4^P(z)(1z)\left(2\mathrm{ln}v+\frac{3}{2}\mathrm{ln}(4z)\right)\right]$$
(48)
where the polynomials $`P_i^P(z)`$ are given in the Table 1 of ref . Taking the contour integral in $`z`$-plane around the cut $`z=[1,\mathrm{})`$, one could write down the following dispersion relation:

$`f^{(1),p}(t)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _{1+ϵ}^{\mathrm{}}}{\displaystyle \frac{\mathrm{Im}f^{(1),p}(z+i0)}{zt}}𝑑z+{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \frac{\pi ^2f_i^p}{(1t)^i}}{\displaystyle \frac{11}{384}}{\displaystyle \frac{ϵ^{3/2}}{1t}}`$ (49)
$`+\left[{\displaystyle \frac{15473}{6912}}+{\displaystyle \frac{35}{36}}\mathrm{ln}(8ϵ)\right]{\displaystyle \frac{ϵ^{1/2}}{1t}}+{\displaystyle \frac{11}{128}}{\displaystyle \frac{ϵ^{1/2}}{(1t)^2}}`$
where $`ϵ0`$ and

$$f_1^p=\frac{11}{256},f_2^p=\frac{69}{256},f_3^p=\frac{197}{2304}.$$
(50)
To simplify the calculation, we represent the imaginary part (48) in the following form:

$$\frac{1}{\pi }\mathrm{Im}f^{(1),p}(z+i0)=F_1^p(z)+F_{2}^{p}{}_{}{}^{}(z)+\frac{1}{2}F_{3}^{p}{}_{}{}^{\prime \prime }(z)$$
(51)
where the functions $`F_p^i(z)`$ grow not faster than $`(z1)^{1/2}`$ at $`z1`$ and have appropriate asymptotic at infinity. Our choice is:

$`F_1^p(z)`$ $`=`$ $`{\displaystyle \frac{1}{96z^2v}}[56+{\displaystyle \frac{49}{3}}z+552z^2+({\displaystyle \frac{207}{4}}60z+904z^24216z^3+3312z^4)\times `$
$`\times {\displaystyle \frac{1}{zv}}\mathrm{ln}{\displaystyle \frac{1v}{1+v}}+({\displaystyle \frac{140}{3}}+{\displaystyle \frac{76}{3}}z{\displaystyle \frac{5120}{3}}z^2+2208z^3)(2\mathrm{ln}v+{\displaystyle \frac{3}{2}}\mathrm{ln}(4z))]`$
$`F_2^p(z)`$ $`=`$ $`{\displaystyle \frac{1}{96z^2v}}[{\displaystyle \frac{11}{4}}+{\displaystyle \frac{645}{2}}z+({\displaystyle \frac{197}{24}}+{\displaystyle \frac{1439}{24}}z){\displaystyle \frac{1}{zv}}\mathrm{ln}{\displaystyle \frac{1v}{1+v}}`$
$`+{\displaystyle \frac{280}{3}}z(2\mathrm{ln}v+{\displaystyle \frac{3}{2}}\mathrm{ln}(4z))]`$
$`F_3^p(z)`$ $`=`$ $`{\displaystyle \frac{1}{96zv}}\left[{\displaystyle \frac{131}{6}}{\displaystyle \frac{197}{12zv}}\mathrm{ln}{\displaystyle \frac{1v}{1+v}}\right]`$ (52)
Then, the r.h.s. of (49) can be integrated by parts twice, all divergent in $`ϵ0`$ terms cancel and the dispersion relation can be brought to the following form:

$$f^{(1),p}(t)=\underset{i=1}{\overset{3}{}}\left[\frac{\pi ^2f_i^p}{(1t)^i}+_1^{\mathrm{}}\frac{F_i^p(z)}{(zt)^i}𝑑z\right]$$
(53)
## Appendix C: $`𝜶_𝒔`$-corrections to the moments

Numerical values of the perturbative corrections to the moments $`\overline{M}_n^{(k),p}/\overline{M}_n^{(0),p}`$ and the $`\alpha _s`$-correction to the condensate contribution $`\overline{M}_n^{(O2)(1),p}/\overline{M}_n^{(O2),p}`$ are given in the Table 1 for $`Q^2/(4\overline{m}^2)=0,1,2`$ and $`n=230`$. The coefficient functions $`\overline{M}_n^p`$ are defined in $`\overline{\mathrm{MS}}`$ scheme according to (25, 26). Remind, that the expansion (25) goes by $`a(\overline{m}^2)`$. If one takes the QCD coupling at another scale, the function $`\overline{M}_n^{(2),p}`$ must be changed according to (27).


format: text

text:
# An adaptive smoothing method for traffic state identification from incomplete information

## 1 Introduction

During the last decades, traffic dynamics and pedestrian flows have been intensively studied. Regarding the observed phenomena and simulation approaches we refer the reader to some recent reviews . Presently, scientists are more and more getting interested in detailed empirical studies. Reasons for this are the availability of better data and the need to verify and calibrate models. Both, theoretical  and empirical  studies indicate that the phenomenology of congested traffic is more complex than originally expected. There seems to be a rich spectrum of traffic states , hysteretic or continuous temporal or spatial transitions among them , and fluctuations or a erratically appearing dynamics play a significant role . In order to make sense out of this, it is increasingly important to have suitable ways of data processing, in order to extract the information relevant for scientific investigations or specific applications. Here, we will propose a three-dimensional data-evaluation method allowing to visualize the spatio-temporal dynamics of traffic patterns along freeways.

The developed “adaptive smoothing method” filters out small-scale fluctuations and adaptively takes into account the main propagation direction of the information flow (i.e. the dominating characteristic line), which have been determined by means of a spatio-temporal correlation analysis in other studies . The temporal length scale of the smoothing procedure can be as small as the sampling interval, while the spatial length scale is related to the distance between successive detectors, which can be up to 3 kilometers long.

By “filter” we just mean a transformation of the data with specific properties. Here, we use a spatio-temporal lowpass filter, i.e., only (Fourier) components of low frequency can pass the filter, while high-frequency contributions are considered as fluctuations and smoothed out. One particular feature of our filter is that it is nonlinear and adaptive to the traffic situation in distinguishing free and congested traffic, as the propagation direction of perturbations differs. The results are three-dimensional visualizations of traffic patterns, which are quite robust with respect to variations of the filter parameters and very helpful in obtaining a clear picture of the systematic spatio-temporal dynamics.

Therefore, this method is suitable for the reconstruction of traffic situations from incomplete information, fast identification of traffic breakdowns (incident detection), and experimental verification of traffic models. First results support the phase diagram of traffic states occuring at bottlenecks , which, apart from free traffic, predicts pinned or moving localized clusters, spatially extended patterns such as triggered stop-and-go waves, and oscillating or homogeneous congested traffic, or a spatial coexistence of some of these states . The respectively occuring spatio-temporal pattern depends on the specific freeway flow and bottleneck strength, but also on the level of fluctuations, as these can trigger transitions from, for example, free traffic to localized cluster states . We do not see sufficient support for one “generalized pattern”  that would always be observed when traffic flow breaks down.

## 2  Description of the Method

The adaptive smoothing method is a data processing method for obtaining traffic variables as smooth functions of space and time out of stationary traffic data. It has following heuristically motivated properties:

1. In case of free traffic, perturbations (of, e.g., velocity or flow) move essentially into the direction of traffic flow . More specifically, they propagate with a characteristic velocity $`c_{\mathrm{free}}`$ at about 80% of the desired velocity $`V_0`$ on empty roads . Therefore, at locations with free traffic, perturbations with propagation velocities near $`c_{\mathrm{free}}`$ should pass the filter.
2. In case of congested traffic, perturbations propagate against the direction of traffic flow with a characteristic and remarkably constant velocity $`c_{\mathrm{cong}}15`$ km . With modern data analysis techniques, it has been shown that such propagation patterns persist even in “synchronized” congested traffic flow, where they are hardly visible in the time series due to a wide scattering of the data in this state . So, for high traffic densities or low velocities, the filter should transmit spatio-temporal perturbations propagating with velocities near $`c_{\mathrm{cong}}`$ more or less unchanged.
3. The filter should smooth out all high-frequency fluctuations in $`t`$ on a time scale smaller than $`\tau `$ and spatial fluctuations in $`x`$ on a length scale smaller than $`\sigma `$. The parameters $`\tau `$ and $`\sigma `$ of the smoothing method can be freely chosen in a wide range (cf. Table 1).

Let us assume that aggregated detector data $`𝒛_{ij}^{\mathrm{in}}`$ are available from $`n`$ cross sections $`i`$ at positions $`x_i`$ where $`j\{j_{\mathrm{min}},\mathrm{},j_{\mathrm{max}}\}`$ denotes the index of the aggregation intervals.

Usually, the aggregation interval

$$\mathrm{\Delta }t=t_jt_{j1}$$
(1)
is fixed (between 20 s and 5 min, depending on the measurement device). A typical value for German highways is $`\mathrm{\Delta }t=1`$ min.

The components of $`𝒛_{ij}^{\mathrm{in}}`$ represent the desired aggregated quantities as obtained from cross section $`i`$ during the time interval $`j`$. Typical examples include the average velocity $`V_{ij}`$, the vehicle flow $`Q_{ij}`$, the occupancy $`O_{ij}`$, or some derived quantity such as the traffic density $`\rho _{ij}`$. The input data $`𝒛_{ij}^{\mathrm{in}}`$ can represent either averages over all lanes or quantities on a given lane.

The adaptive smoothing method provides estimates $`𝒛(x,t)`$ for all locations $`x[x_1,x_n]`$ between the positions of the first and the last detector, and for all times $`t[t_{\mathrm{min}},t_{\mathrm{max}}]`$. Extrapolations (in the sense of a short-term traffic forecast) will also be discussed. Without loss of generality, we assume $`x_1<x_2<\mathrm{}<x_n`$, and a traffic flow in positive $`x`$-direction.

The core of our “adaptive smoothing method” is a nonlinear filter transforming the discrete input detector data $`𝒛_{ij}^{\mathrm{in}}`$ into the smooth spatio-temporal functions $`𝒛(x,t)`$. To satisfy the first two requirements mentioned above, we write the filter as

$$𝒛(x,t)=w(𝒛_{\mathrm{cong}},𝒛_{\mathrm{free}})𝒛_{\mathrm{cong}}(x,t)+\left[1w(𝒛_{\mathrm{cong}},𝒛_{\mathrm{free}})\right]𝒛_{\mathrm{free}}(x,t).$$
(2)
This is a superposition of two linear anisotropic lowpass filters $`𝒛_{\mathrm{cong}}(x,t)`$ and $`𝒛_{\mathrm{free}}(x,t)`$ with an adaptive weight factor $`0w1`$ which itself depends nonlinearly on the output of the linear filters as discussed later on.

The filter $`𝒛_{\mathrm{cong}}(x,t)`$ for congested traffic is given by

$$𝒛_{\mathrm{cong}}(x,t)=\frac{1}{N_{\mathrm{cong}}(x,t)}\underset{i=1}{\overset{n}{}}\underset{j=j_{\mathrm{min}}}{\overset{j_{\mathrm{max}}}{}}\varphi _{\mathrm{cong}}(x_ix,t_jt)𝒛_{ij}^{\mathrm{in}},$$
(3)
with

$$t_j=t_{\mathrm{min}}+j\mathrm{\Delta }t,$$
(4)
and the normalization factor

$$N_{\mathrm{cong}}(x,t)=\underset{i=1}{\overset{n}{}}\underset{j=j_{\mathrm{min}}}{\overset{j_{\mathrm{max}}}{}}\varphi _{\mathrm{cong}}(x_ix,t_jt).$$
(5)
This normalization guarantees that a constant input vector $`𝒛_{ij}^{\mathrm{in}}=𝒛_0`$ for all $`i,j`$ is transformed into a constant output function with the same components: $`𝒛(x,t)=𝒛_0`$. Notice that the normalization depends on both, location $`x`$ and time $`t`$.

The analogous expression for the lowpass filter $`𝒛_{\mathrm{free}}(x,t)`$ for free traffic is

$$𝒛_{\mathrm{free}}(x,t)=\frac{1}{N_{\mathrm{free}}(x,t)}\underset{i=1}{\overset{n}{}}\underset{j=j_{\mathrm{min}}}{\overset{j_{\mathrm{max}}}{}}\varphi _{\mathrm{free}}(x_ix,t_jt)𝒛_{ij}^{\mathrm{in}}$$
(6)
with

$$N_{\mathrm{free}}(x,t)=\underset{i=1}{\overset{n}{}}\underset{j=j_{\mathrm{min}}}{\overset{j_{\mathrm{max}}}{}}\varphi _{\mathrm{free}}(x_ix,t_jt).$$
(7)
The kernels $`\varphi _{\mathrm{cong}}(x,t)`$ and $`\varphi _{\mathrm{free}}(x,t)`$ of the linear homogeneous filters do the required smoothing and are particularly transmissible for perturbations propagating with the typical velocities $`c_{\mathrm{cong}}`$ and $`c_{\mathrm{free}}`$ observed in congested and free traffic, respectively (cf. Fig. 1).

If the propagation velocity of the perturbations was zero, the filter should only perform the smoothing. We chose the (non-normalized) exponential function

$$\varphi _0(x,t)=\mathrm{exp}\left(\frac{|x|}{\sigma }\frac{|t|}{\tau }\right).$$
(8)
Instead, one could apply other localized functions such as a two-dimensional Gaussian. However, it turned out that the exponential had more favourable properties in our application.

The linear filter for nonzero propagation velocities can be mapped to $`\varphi _0(x,t)`$ by the coordinate transformations (cf. Fig. 1)

$$x=x^{},t=t_{\mathrm{cong}}^{}+\frac{x}{c_{\mathrm{cong}}}=t_{\mathrm{free}}^{}+\frac{x}{c_{\mathrm{free}}}.$$
(9)
Thus, we obtain for the kernels of the linear anisotropic filters

$`\varphi _{\mathrm{cong}}(x,t)`$ $`=`$ $`\varphi _0(x^{},t_{\mathrm{cong}}^{})=\varphi _0(x,t{\displaystyle \frac{x}{c_{\mathrm{cong}}}}),`$ (10)
$`\varphi _{\mathrm{free}}(x,t)`$ $`=`$ $`\varphi _0(x^{},t_{\mathrm{free}}^{})=\varphi _0(x,t{\displaystyle \frac{x}{c_{\mathrm{free}}}}).`$ (11)
Figure 2 shows the action of the filters $`𝒛_{\mathrm{cong}}(x,t)`$ and $`𝒛_{\mathrm{free}}(x,t)`$ for the velocity fields $`V_{\mathrm{cong}}(x,t)`$ and $`V_{\mathrm{free}}(x,t)`$.

Finally, we define the nonlinear adaptive weight function $`w(𝒛_{\mathrm{cong}},𝒛_{\mathrm{free}})[0,1]`$. Obviously, we must have $`w1`$ for congested traffic, and $`w0`$ for free traffic, so we need some a priori estimate of the traffic situation at the point $`(x,t)`$.

Congested traffic is characterized by a high traffic density and low average velocity. Since, in contrast to the density, the velocity can be directly measured with stationary detectors, we chose the velocity to determine the a priori estimate. Different possibilities to estimate the velocity at point $`(x,t)`$ are:

* The measured velocity $`V_{ij}`$ of the detector cross section whose position $`x_i`$ is nearest to $`x`$ in the time interval $`j`$ containing the actual time $`t`$,
* the velocity $`V_{\mathrm{cong}}(x,t)`$ as calculated with the “congested-traffic” filter $`𝒛_{\mathrm{cong}}`$ according to Eq. (3) with kernel (10) (assuming a priori congested traffic),
* the velocity $`V_{\mathrm{free}}(x,t)`$ as calculated with the “free-traffic” filter $`𝒛_{\mathrm{free}}`$ according to Eq. (6) with the kernel (11),
* or some combination of the above estimates.

The first way to estimate the velocity is subject to errors, if the typical length scale $`\lambda `$ of occuring stop-and-go structures is not larger than the distance $`\mathrm{\Delta }x_i`$ between two neighbouring detectors. At this point, it is crucial that propagating structures in congested traffic, especially stop-and-go waves, are very persistent. It has been shown  that they can propagate through freeway intersections or other inhomogeneities nearly unchanged, passing all perturbations of free traffic on their way (see, e.g., Fig. 9). Therefore, whenever at least one of the estimates indicates congested traffic, the weight function should favour the filter for congested traffic. Specifically, we assume

$$w(𝒛_{\mathrm{cong}},𝒛_{\mathrm{free}})=w(V_{\mathrm{cong}}(x,t),V_{\mathrm{free}}(x,t))=\frac{1}{2}\left[1+\mathrm{tanh}\left(\frac{V_cV^{}(x,t)}{\mathrm{\Delta }V}\right)\right],$$
(12)
where

$$V^{}(x,t)=\mathrm{min}(V_{\mathrm{cong}}(x,t),V_{\mathrm{free}}(x,t)).$$
(13)
$`V_c`$ and $`\mathrm{\Delta }V`$ are parameters that can be varied in a wide range.

If not explicitely stated otherwise, for all simulations in the following section we will use the parameters specified in Table 1. We will also show that the four parameters $`c_{\mathrm{cong}}`$, $`c_{\mathrm{free}}`$, $`V_c`$ and $`\mathrm{\Delta }v`$ can be varied in a wide range without great differences in the output. In this way, we show that the proposed adaptive smoothing method does not need to be calibrated to the respective freeway. One can take the values from Table 1 as a global setting.

The smoothing parameters $`\sigma `$ and $`\tau `$ have the same meaning and the same effect as in standard smoothing methods, e.g., Eqs. (22) and (23) in Ref. .

In summary, the proposed “adaptive smoothing method” is given by Eq. (2) with the nonlinear weight function (12), the filter (3) with normalization (5) and kernel (10) for congested traffic, the filter (6) with normalization (7) and kernel (11) for free traffic, and the smoothing filter (8). Table 1 gives an overview of the six parameters involved and typical values for them. The adaptive smoothing method includes the following special cases:

* Isotropic smoothing resulting in the limits $`c_{\mathrm{free}}\mathrm{}`$ and $`c_{\mathrm{cong}}\mathrm{}`$ (for practical purposes, one may chose $`c_{\mathrm{free}}=c_{\mathrm{cong}}=10^6`$ km/h);
* only filtering for structures of congested traffic in the limit $`V_c\mathrm{}`$ (for practical purposes, one may set $`V_cV_0`$ with $`V_0`$ being the desired velocity);
* only filtering for structures of free traffic for $`V_c=\mathrm{\Delta }V=0`$;
* consideration of the data of the nearest detector only, if $`\sigma =0`$;
* application of the actual sampling interval of the detectors, if $`\tau =0`$, $`c_{\mathrm{free}}\mathrm{}`$ and $`c_{\mathrm{cong}}\mathrm{}`$. (To avoid divisions by zero, zero values of $`\sigma `$, $`\tau `$, and $`\mathrm{\Delta }V`$ are replaced by very small positive values in the software).

## 3 Application to German Freeways

We will now discuss data from the German freeways A8-East and A9-South near Munich, and of the freeway A5-North near Frankfurt. In all cases, the traffic data were obtained from several sets of double-induction-loop detectors recording, separately for each lane, the passage times and velocities of all vehicles. Only aggregate information was stored with an aggregation interval of $`\mathrm{\Delta }t=1`$ min. We will use the following input data $`𝒛_{ij}^{\mathrm{in}}`$:

* The lane-averaged vehicle flow

$$Q_{ij}=\underset{l=1}{\overset{L}{}}\frac{n_{ij}^l}{L\mathrm{\Delta }t},$$
(14)
where $`n_{ij}^l`$ is the vehicle count at cross section $`i`$ during time interval $`j`$ on lane $`l`$. The considered sections of the freeways have $`L=3`$ lanes in most cases.
* The lane-averaged mean velocity

$$V_{ij}=\underset{l=1}{\overset{L}{}}\frac{Q_{ij}^lV_{ij}^l}{Q_{ij}},$$
(15)
where $`V_{ij}^l`$ is the average velocity at cross section $`i`$ during time interval $`j`$ on lane $`l`$.
* The traffic density determined via the formula

$$\rho _{ij}=\frac{Q_{ij}}{V_{ij}},$$
(16)
as occupancies were not available for all freeways.

Figure 3 shows an example of a complex traffic breakdown that occurred on the freeway A8 East from Munich to Salzburg during the evening rush hour on November 2, 1998. Two different kinds of bottlenecks were relevant, (i) a relatively steep uphill gradient from $`x=38`$ km to $`x=40`$ km (the “Irschenberg”), and (ii) an incident leading to the closing of one of the three lanes between the cross sections D23 and D24 between $`t=`$ 17:40 h and $`t=`$ 18:10 h. For further details, see Ref. . Figure 3(c) shows the traffic density using the proposed data processing method with the parameters specified in Table 1. For comparison, (b) shows the result of conventional isotropic smoothing ignoring the speed of information propagation, cf. Ref. . While the structures in the free-traffic regions are nearly the same in both plots, the new smoothing method gives a better reconstruction of the congested traffic patterns. In contrast to the conventional method, our proposed method can resolve oscillations with wavelengths $`\lambda `$ comparable to or smaller than the distance $`\mathrm{\Delta }x_i`$ between two neighbouring cross sections as illustrated for the three stop-and-go waves propagating through the whole displayed section ($`\lambda 2`$ km). Notice that the conventional method generates artefacts for this example.

The congestion patterns caused by the uphill gradient (the “wings” in Fig. 3(b), (c)) have larger characteristic wavelengths and are resolved, at least partially, also with the conventional method.

For comparison, we plot time series of some cross sections in Figure 4. Let us consider more closely the jammed traffic patterns caused by the temporary bottleneck associated with the incident. While the fluctuations are irregular at detector D22 (about 1.2 km upstream of the aforementioned incident), they grow to three stop-and-go waves further upstream, cf. the arrows in Fig. 5. However, due to measurement errors and other errors in determining the density , the time series of the density is very noisy, especially that of cross section D20. The adaptive smoothing method suppresses this noise without suppressing the structure of the pattern.

To check, whether there are significant patterns in the data of detector D20 at all, we plot the velocity data at D20 as well (Fig. 5). In contrast to the density, the velocity is measured directly resulting in less noise. One clearly sees three dips corresponding to three stop-and-go waves which confirms the results of the “adaptive smoothing method”.

As a further check that the structures shown in Fig. 3(c) are not artefacts of the data processing, we tested the method for a reduced data set. Figure 6 shows the result using only the even-numbered detectors as input, whose time series are plotted in Fig. 4. The main structures remain nearly unchanged, particularly, the three stop-and-go waves propagating through the whole section. Notice that the whole plot is based on data of only four cross sections, and that the data are extrapolated in the upstream direction over a distance of nearly 2 km.

We now apply the processing method to data from a section of the freeway A9-South near Munich. There are two major intersections I1 and I2 with other freeways. Virtually every weekday, traffic broke down to oscillatory congested traffic upstream of each of these intersections. (For details, see Ref. ).

Figure 7 (b) shows the inverse of the velocity, $`1/V(x,t)`$ for a typical congested situation. When plotting $`1/V(x,t)`$, the structures of congested traffic comes out more clearly than for the density. The method resolves small density clusters in the region 508 km $`x510`$ km between about 8:30 am and 9:30 am which disappear around $`x=508`$ km. These structures cannot be identified with the conventional smoothing method (see Fig. 7(c)), which just shows a hilly pattern.

Figure 7(b) is an example for the spatial coexistence of homogeneous congested traffic (a short stretch around $`x=510`$ km), oscillating congested traffic (around $`x=509`$ km), and stop-and-go waves (around $`x=507`$ km). This may be a three-dimensional illustration of the so-called “pinch effect” , but we do not observe a merging of narrow clusters to form wide moving jams. The narrow structures of short wavelengths rather disappear. This may be an artefact of our smoothing method, so that video data are required to get a more detailed picture. Simulations of this spatial coexistence, however, indicate that narrow clusters rather disappear than merge, namely when they do not exceed the critical amplitude in an area of metastable traffic .

Notice that the bottleneck causing this congestion is located at about $`x=510`$ km implying that it is caused by weaving traffic and slowing down at an off-ramp, not by an on-ramp as is most often the case on German highways.

Now we demonstrate that the method is robust with respect to reasonable parameter changes. Figure 8(a) shows the result for an assumed propagation velocity of $`c_{\mathrm{cong}}=12`$ km/h instead of -15 km/h. In plot (b), we assumed $`c_{\mathrm{cong}}=20`$ km/h instead. In both cases, the results are systematically better than with the isotropic procedure. Since the propagation velocity of perturbations in congested traffic is always about $`c_{\mathrm{cong}}=15`$ km , one can take this value as a global setting.

Due to the comparatively high magnitude of the propagation velocity $`c_{\mathrm{free}}`$, the related time shifts in the transformation (9) for free traffic are small. Therefore, the method is insensitive to changes with respect to this parameter. It turned out that taking the isotropic limit $`c_{\mathrm{free}}\mathrm{}`$ nearly gives the same results.

The crossover parameters $`V_c`$ and $`\mathrm{\Delta }V`$ can be varied within reasonable limits as well. As an example, Fig. 8(c) shows the result after changing the crossover velocity from $`V_c=60`$ km/h to $`V_c=40`$ km/h.

Finally, we applied the method to isolated moving localized clusters and pinned localized clusters observed on the A5-North near Frankfurt. This freeway is particularly well equipped with detectors, so that the typical length scales of the observed structures are always larger than the distance between neighbouring detectors. However, we obtained nearly the same results when using only the even-numbered detectors or only the odd-numbered detectors as input. Figure 9 shows that even a further reduction of the information to only 35% of the detectors yields good results.

## 4 Summary and Outlook

We have proposed a new “adaptive smoothing method” for the three-dimensional visualization of spatio-temporal traffic patterns, which takes into account the characteristic propagation velocities observed in free and congested traffic. The method is robust with respect to variations of its parameters, so that it can be applied to new freeway sections without calibration. In principle, it would be possible to determine the parameters (such as the propagation velocity of perturbations) locally (e.g., by means of a correlation analysis). However, the results would look less smooth and regular, as the small number of data to determine the local parameters would be associated with considerable statistical errors. Consequently, a large part of the variations in the local parameters would not reflect systematic variations of the parameters. So, both the local and global parameter calibration may produce artefacts, but it is advantageous to use global parameter settings. Moreover, there is empirical support for surprisingly constant propagation velocities $`c_{\mathrm{free}}`$ and $`c_{\mathrm{cong}}`$ of perturbations in free and congested traffic. The parameters $`V_c`$ and $`\mathrm{\Delta }V`$ are related to the transition from free to congested traffic and, therefore, can also be well determined. The spatial and temporal smoothing parameters $`\sigma `$ and $`\tau `$ can be specified according to the respective requirements. Suitable parameters allow a good representation of traffic patterns even when the distances between successive detectors are about 3 kilometers.

We point out that the suggested “adaptive smoothing method” itself can be applied to calibrate the characteristic propagation velocities $`c_{\mathrm{free}}`$ and $`c_{\mathrm{cong}}`$. For this purpose, the ranges $`\sigma `$ and $`\tau `$ of spatial and temporal smoothing are chosen small. The optimal propagation velocities minimize the offsets in the propagation patterns.

The aim of the “adaptive smoothing method” is to reconstruct the spatio-temporal traffic data from incomplete information as good as possible to allow a better understanding of the complex traffic dynamics. Potential applications are, for example, traffic state visualization, incident detection, or the experimental verification of traffic models. Our method could be further improved by taking into account information about the traffic dynamics such as the continuity equation or a suitable equation for the average vehicle velocity as a function of space and time. It could, then, be used to determine short-term traffic forecasts.

### Acknowledgments

The authors would like to thank the German Research Foundation (DFG) for financial support through the grant He 2789/2-1. They are also grateful to the Autobahndirektion Südbayern, and the Hessisches Landesamt für Straßen und Verkehrswesen for providing the traffic data.


text:
# Ground state degeneracy of the Pauli-Fierz Hamiltonian with spin

## 1 Introduction and main results

Nonrelativistic quantum electrodynamics predicts the anomalous magnetic moment ($`g`$-factor) of the electron with an error of less than $`1\%`$. One definition of the $`g`$-factor comes from the motion of the electron in a weak uniform external magnetic field. $`g/2`$ is then the ratio of the spin precession relative to the orbital precession. To provide a theoretical analysis of such an experiment one has to show that the ground state band of the Hamiltonian is adiabatically decoupled from the spectrum of excitations and has to determine the effective Hamiltonian governing the motion in the ground state band. The details are given elsewhere . They result in a nonperturbative microscopic definition of the $`g`$-factor, which, when computed to order $`e^2`$, yields $`g/2=1.0031`$ as compared to the accurate relativistic value of $`g/2=1.0012`$. One crucial input of the analysis is the microscopic Hamiltonian at fixed total momentum to have an exactly two-fold degenerate ground state. In our paper we will, under suitable conditions, establish such a spectral property.

Let us first introduce the Pauli-Fierz Hamiltonian for a single electron coupled to the quantized radiation field. By translation invariance the total momentum is conserved and the operator of interest will be the Pauli-Fierz Hamiltonian fibered with respect to the total momentum. The Hilbert space for the coupled system is

$$=L^2(R^3)C^2.$$

$`L^2(R^3)`$ is the Hilbert space for the translational degrees of freedom of the electron, multiplication by $`x`$ stands for the position, $`i_x`$ for the momentum of the electron. $`C^2`$ is the Hilbert space for the spin, $`\sigma `$, of the electron, where $`\sigma =(\sigma _1,\sigma _2,\sigma _3)`$ are the spin $`1/2`$ Pauli spin matrices,

$$\sigma _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).$$

$``$ is the symmetric Fock space for the photons given by

$$=_{n=0}^{\mathrm{}}\left(L^2(R^3)C^2\right)_{\mathrm{sym}}^n,$$

where $`(\mathrm{})_{\mathrm{sym}}^n`$ denotes the $`n`$-fold symmetric tensor product of $`(\mathrm{})`$. $`L^2(R^3)C^2`$ will be identified with $`L^2(R^3\times \{1,2\})`$. The Fock vacuum is denoted by $`\mathrm{\Omega }`$ and a vector $`\psi `$ by $`\psi =(\psi ^{(0)},\psi ^{(1)},\mathrm{})`$, $`\psi ^{(0)}=c\mathrm{\Omega }`$, $`cC`$. The photons live in $`R^3`$ and have helicity $`\pm 1`$. The photon field is thus represented by the two-component Bose field $`a(k,j)`$ on $``$ with commutation relations

$$[a(k,j),a^{}(k^{},j^{})]=\delta _{jj^{}}\delta (kk^{}),$$

$$[a(k,j),a(k^{},j^{})]=0,[a^{}(k,j),a^{}(k^{},j^{})]=0,$$

$`k,k^{}R^3`$, $`j,j^{}=1,2`$. The energy of the photons is given by

$$H_\mathrm{f}=\underset{j=1,2}{}\omega (k)a^{}(k,j)a(k,j)𝑑k,$$
(1)
i.e. $`H_\mathrm{f}`$ restricted to $`(L^2(R^3)C^2)_{\mathrm{sym}}^n`$ is multiplication by $`_{j=1}^n\omega (k_j)`$. Throughout units are such that $`\mathrm{}=1`$, $`c=1`$. Physically $`\omega (k)=|k|`$. This case is somewhat singular. To regularize, the photons are assumed to have a small mass as

$$\omega (k)=\sqrt{|k|^2+m_{\mathrm{ph}}^2},m_{\mathrm{ph}}>0.$$
(2)
In fact we can allow for a more general class of dispersion relations. We assume $`\omega :R^3R`$ to be continuous with the properties

$`inf_{kR^3}\omega (k)\omega _0>0,`$

$`\omega (k_1)+\omega (k_2)\omega (k_1+k_2),`$

$`\omega (k)=\omega (Rk)`$ for an arbitrary rotation $`R`$.

The quantized transverse vector potential is defined through

$$A_\phi (x)=\underset{j=1,2}{}\frac{\widehat{\phi }(k)}{\sqrt{2\omega (k)}}e_j(k)\left(e^{ikx}a^{}(k,j)+e^{ikx}a(k,j)\right)𝑑k.$$

Here $`e_1`$ and $`e_2`$ are polarization vectors which together with $`\widehat{k}=k/|k|`$ form a standard basis in $`R^3`$. $`\phi :R^3R`$ is the form factor which ensures an ultraviolet cutoff. It is assumed to be rotational invariant, $`\phi (Rx)=\phi (x)`$ for an arbitrary rotation $`R`$, continuous, bounded with some decay at infinity, and normalized as $`\phi (x)𝑑x=1`$. We will mostly work with the Fourier transform $`\widehat{\phi }(k)=(2\pi )^{3/2}\phi (x)e^{ikx}𝑑x`$. It satisfies (1) $`\widehat{\phi }(Rk)=\widehat{\phi }(k)`$ for an arbitrary rotation $`R`$, (2) $`\overline{\widehat{\phi }}=\widehat{\phi }`$ for notational simplicity, (3) the normalization $`\widehat{\phi }(0)=(2\pi )^{3/2}`$, and (4) the decay

$$\left(\omega (k)^2+\omega (k)^1+1+\omega (k)\right)|\widehat{\phi }(k)|^2𝑑k<\mathrm{}.$$

With these preparations the Pauli-Fierz Hamiltonian, including spin, is defined by

$$H=\frac{1}{2m}\left\{\sigma \left(i_xeA_\phi (x)\right)\right\}^2+H_\mathrm{f},$$
(3)
where $`m`$ is the bare mass and $`e`$ the charge of the electron.

Translations for the electron are generated by $`i_x`$ and translations for the photon field by the field momentum

$$P_\mathrm{f}=\underset{j=1,2}{}ka^{}(k,j)a(k,j)𝑑k.$$
(4)
By translation invariance of $`H`$ the total momentum

$$p=i_x+P_\mathrm{f}$$

is thus conserved,

$$[H,p]=0,$$

as can be checked also directly. We unitarily transform $`H`$ such that the fibering with respect to $`p`$ becomes transparent. In momentum representation, momentum multiplication by $`k`$, position $`i_k`$, an element $`\psi `$ is written as $`\psi ^{(n)}(k,k_1,j_1,\mathrm{},k_n,j_n)`$ with values in $`C^2`$. For each $`n`$ let

$$U\psi ^{(n)}(k,k_1,j_1,\mathrm{},k_n,j_n)=\psi ^{(n)}(k\underset{j=1}{\overset{n}{}}k_j,k_1,j_1,\mathrm{},k_n,j_n)$$
(5)
with inverse

$$U^1\psi ^{(n)}(k,k_1,j_1,\mathrm{},k_n,j_n)=\psi ^{(n)}(k+\underset{j=1}{\overset{n}{}}k_j,k_1,j_1,\mathrm{},k_n,j_n).$$
(6)
Clearly $`U`$ is unitary. We set

$$A_\phi =A_\phi (0)$$

as an operator on $``$. Similarly for the quantized magnetic field,

$$B_\phi (x)=i\underset{j=1,2}{}\frac{\widehat{\phi }(k)}{\sqrt{2\omega (k)}}(ke_j(k))\left(e^{ikx}a^{}(k,j)e^{ikx}a(k,j)\right)𝑑k,$$

and we denote

$$B_\phi =B_\phi (0)$$

as an operator on $``$. Then, first working out the square in (3), one obtains

$$UHU^1=\frac{1}{2m}(pP_\mathrm{f}eA_\phi )^2\frac{e}{2m}\sigma B_\phi +H_\mathrm{f}.$$
(7)
In (7) $`p`$ is multiplication by $`k`$ in the representation from (5) and (6). Thus the fibering is

$$=_{R^3}^{}_p𝑑p$$

with $`_p`$ isomorphic to $`C^2`$ and

$$UHU^1=_{R^3}^{}H_p𝑑p.$$

In the following we will regard $`pR^3`$ simply as a parameter. We also choose units such that $`m=1`$. Then the operator under study is

$$H_p=\frac{1}{2}\left(pP_\mathrm{f}eA_\phi \right)^2\frac{e}{2}\sigma B_\phi +H_\mathrm{f}$$
(8)
acting on

$$_p=C^2.$$

For simplicity the index $`p`$ of $`_p`$ will be omitted.

For $`e=0`$ (8) reduces to the noninteracting Hamiltonian

$$H_{p0}=\frac{1}{2}(pP_\mathrm{f})^2+H_\mathrm{f}.$$
(9)
Clearly, with the definitions (1) and (4), $`H_{p0}`$ is self-adjoint with the domain $`D(H_\mathrm{f}+P_{\mathrm{f}}^{}{}_{}{}^{2})=D(H_\mathrm{f})D(P_{\mathrm{f}}^{}{}_{}{}^{2})`$. The rest of (8) is regarded as the interaction part of the Hamiltonian,

$$H_{p\mathrm{I}}=H_pH_{p0}=e(pP_\mathrm{f})A_\phi +\frac{e^2}{2}A_\phi ^2\frac{e}{2}\sigma B_\phi .$$
(10)
For sufficiently small $`e`$, $`|e|<e^{}`$, $`H_{p\mathrm{I}}`$ is bounded relative to $`H_{p0}`$ with a bound less than 1, which by a theorem of Kato and Rellich implies that $`H_p`$ is a self-adjoint operator for every $`pR^3`$. In addition $`H_p`$ is bounded from below. To see this let $`a^{\mathrm{}}(f)=_{j=1,2}f(k,j)a^{\mathrm{}}(k,j)𝑑k`$. Using the inequalities

$$a^{\mathrm{}}(f)\psi c_1\left\{\underset{j=1,2}{}\left(\omega (k)^1+1\right)|f(k,j)|^2𝑑k\right\}^{1/2}(H_\mathrm{f}+1)^{1/2}\psi ,$$

$$a^{\mathrm{}}(f)a^{\mathrm{}}(f)\psi c_2\underset{j=1,2}{}\left(\omega (k)^1+\omega (k)\right)|f(k,j)|^2𝑑k(H_\mathrm{f}+1)\psi ,$$

with some constants $`c_1`$ and $`c_2`$, one has

$$H_{p\mathrm{I}}\psi c_0(e)(H_{p0}+1)\psi .$$
(11)
Here

$$c_0(e)=c_0\{|e|\{(\omega (k)^2+\omega (k))|\widehat{\phi }(k)|^2dk\}^{1/2}$$

$$+|e|^2(\omega (k)^2+1)|\widehat{\phi }(k)|^2dk\}$$

with some numerical constant $`c_0`$ of order one. Thus $`|e|<e^{}`$ with a suitable $`e^{}>0`$ implies $`c_0(e)<1`$.

The goal of our paper is to study some ground state properties of $`H_p`$. The ground state energy of $`H_p`$ is

$$E(p)=inf\mathrm{Spec}(H_p)=\underset{\psi D(H_p),\psi =1}{inf}(\psi ,H_p\psi ).$$

It is easily seen that resolvent $`(H_pz)^1`$ with $`z\mathrm{Spec}(H_p)`$ is continuous in both $`p`$ and $`e`$ in the operator norm. Thus $`E(p)`$ is continuous in both $`p`$ and $`e`$. If $`E(p)`$ is an eigenvalue, the corresponding spectral projection is denoted by $`P_\mathrm{g}`$. $`\mathrm{Tr}P_\mathrm{g}`$ is the degeneracy of the ground state. The bottom of the continuous spectrum is denoted by $`E_\mathrm{c}(p)`$. Under the assumptions (A.1)–(A.3) one knows that

$$E_\mathrm{c}(p)=\underset{kR^3}{inf}\left\{E(pk)+\omega (k)\right\},$$

see . Thus it is natural to set

$$\mathrm{\Delta }(p)=E_\mathrm{c}(p)E(p)=\underset{kR^3}{inf}\left\{E(pk)+\omega (k)E(p)\right\}.$$

Let $`e=0`$. Then $`\mathrm{\Delta }(p)>0`$ for $`|p|<p_0`$ with some $`p_0`$ by (A.1). Thus, by the continuity of $`E(p)`$ mentioned above, $`\{(p,e)R^3\times R|\mathrm{\Delta }(p)>0\}\mathrm{}.`$ As our main result we state

###### Theorem 1.1.

Suppose $`|e|<e_0`$ and $`\mathrm{\Delta }(p)>0`$. Then $`\mathrm{Tr}P_\mathrm{g}=2.`$

We briefly comment on our assumptions. $`\mathrm{\Delta }(p)>0`$ means that we assume, rather than prove, a spectral gap. If instead one would merely impose the implicit but rather natural condition that

$$E(p)E(0),$$
(12)
then, following the arguments in , one concludes in case of the dispersion relation (2) that $`\mathrm{\Delta }(p)>0`$ for $`|p|<\sqrt{3}1`$ and arbitrary $`eR`$. For a general dispersion relation satisfying (A.1) to (A.3) a corresponding bound can be established. For the spinless Pauli-Fierz Hamiltonian (12) is proven through a suitable variant of a diamagnetic inequality. We did not succeed to establish (12) for the Hamiltonian (8). Note that on physical grounds one expects $`\mathrm{\Delta }(p)>0`$ for $`|p|<p_\mathrm{c}`$ and $`\mathrm{\Delta }(p)=0`$ for $`|p|>p_\mathrm{c}`$ with a simultaneous loss of the ground state. For the Hamiltonian (8) with $`\omega `$ given by (2) one has $`p_\mathrm{c}<\mathrm{}`$ at $`e=0`$ and expects $`p_\mathrm{c}<\mathrm{}`$ to persist for all $`e0`$ provided the spatial dimension $`d3`$.

The constant $`e_0`$ does not depend on $`m_{\mathrm{ph}}`$. Thus, together with (12), the domain of validity of Theorem 1.1 is independent of $`m_{\mathrm{ph}}`$ and it may seem that one could take the limit $`m_{\mathrm{ph}}0`$. Of course, thereby $`\mathrm{\Delta }(p)0`$. Unfortunately, this will not work, since the Pauli-Fierz Hamiltonian is infrared divergent and the number of photons increases without bound as $`m_{\mathrm{ph}}0`$ . The physical ground states are no longer in Fock space. The only exception is $`p=0`$, and one might hope to apply our method directly to the case $`p=0`$ and $`m_{\mathrm{ph}}=0`$. This problem will be taken up in Section 3.

Let us indicate the strategy for the proof of Theorem 1.1. From a pull-through formula one estimates the overlap of a ground state with the subspace $`P_0=C^2\{C\mathrm{\Omega }\}`$. If $`|e|`$ is sufficiently small, the overlap is large which implies $`\mathrm{Tr}P_\mathrm{g}<3`$. For a lower bound we will derive the algebraic relation $`P_0P_\mathrm{g}P_0=aP_0`$, $`a>0`$, which implies $`\mathrm{Tr}P_\mathrm{g}2`$.

The Hamiltonian (8) is invariant under rotations with axis $`\widehat{p}=p/|p|`$. To understand the implications let us define the field angular momentum relative to the origin by

$$J_\mathrm{f}=i\underset{j=1,2}{}a^{}(k,j)(k_k)a(k,j)𝑑k$$

and the helicity by

$$S_\mathrm{f}=i\widehat{k}\left\{a^{}(k,2)a(k,1)a^{}(k,1)a(k,2)\right\}𝑑k.$$

Let $`\stackrel{}{n}R^3`$ be an arbitrary unit vector, $`\theta R`$, and let $`R=R(\stackrel{}{n},\theta )`$ be the rotation around $`\stackrel{}{n}`$ through the angle $`\theta `$. We note that

$$e^{i\theta \stackrel{}{n}(J_\mathrm{f}+S_\mathrm{f})}H_\mathrm{f}e^{i\theta \stackrel{}{n}(J_\mathrm{f}+S_\mathrm{f})}=H_\mathrm{f},$$

$$e^{i\theta \stackrel{}{n}(J_\mathrm{f}+S_\mathrm{f})}P_\mathrm{f}e^{i\theta \stackrel{}{n}(J_\mathrm{f}+S_\mathrm{f})}=RP_\mathrm{f},$$

$$e^{i\theta \stackrel{}{n}(J_\mathrm{f}+S_\mathrm{f})}A_\phi e^{i\theta \stackrel{}{n}(J_\mathrm{f}+S_\mathrm{f})}=RA_\phi .$$

Define the total angular momentum through

$$J=J_\mathrm{f}+S_\mathrm{f}+\frac{1}{2}\sigma .$$

Then it follows that

$$e^{i\theta \stackrel{}{n}J}H_pe^{i\theta \stackrel{}{n}J}=H_{R^1p}.$$

In particular, the ground state energy is rotation invariant,

$$E(p)=E(Rp).$$
(13)
If one sets $`\stackrel{}{n}=\widehat{p}`$, then

$$e^{i\theta \widehat{p}J}H_pe^{i\theta \widehat{p}J}=H_p,$$
(14)
which expresses the rotation invariance relative to $`\widehat{p}`$.

Since $`\mathrm{Spec}(\widehat{p}(J_\mathrm{f}+S_\mathrm{f}))=Z`$ and $`\mathrm{Spec}(\widehat{p}\sigma )=\{1,+1\}`$, it follows that

$$\mathrm{Spec}(\widehat{p}J)=Z_{1/2},$$

where $`Z_{1/2}`$ is the set of half integers $`\{\pm 1/2,\pm 3/2,\pm 5/2,\mathrm{}\}`$. By virtue of (14), $``$ and $`H_p`$ are decomposable as

$$=\underset{zZ_{1/2}}{}(z),$$

$$H_p=\underset{zZ_{1/2}}{}H_p(z).$$

Therefore, if $`\mathrm{Tr}P_\mathrm{g}=2`$, the ground states $`\psi _{\mathrm{g}\pm }P_\mathrm{g}`$ can be chosen such that $`\psi _{\mathrm{g}+}(z)`$ and $`\psi _\mathrm{g}(z^{})`$ for some $`z,z^{}Z_{1/2}`$.

###### Theorem 1.2.

Suppose $`|e|<e_0`$ and $`\mathrm{\Delta }(p)>0`$. Then $`H_p`$ has two orthogonal ground states, $`\psi _{\mathrm{g}\pm }`$, with the property $`\psi _{\mathrm{g}\pm }(\pm 1/2)`$.

## 2 Spectral properties

### 2.1 Upper bound

Let us denote the number operator by

$$N_\mathrm{f}=\underset{j=1,2}{}a^{}(k,j)a(k,j)𝑑k.$$

In what follows $`\psi _\mathrm{g}`$ denotes an arbitrary normalized ground state of $`H_p`$. We note that $`a(k,j)\psi `$, $`\psi D(N_{\mathrm{f}}^{}{}_{}{}^{1/2})`$, is well defined and

$$(a(k,j)\psi )^{(n)}(k_1,j_1,\mathrm{},k_n,j_n)=\sqrt{n+1}\psi ^{(n+1)}(k,j,k_1,j_1,\mathrm{},k_n,j_n).$$

Moreover it follows that $`a(k,j)\psi N_{\mathrm{f}}^{}{}_{}{}^{1/2}\psi `$ and

$$(\psi ,N_\mathrm{f}\varphi )=\underset{j=1,2}{}(a(k,j)\psi ,a(k,j)\varphi )𝑑k.$$

###### Lemma 2.1.

Suppose $`\mathrm{\Delta }(p)>0`$. Then

$$(\psi _\mathrm{g},N_\mathrm{f}\psi _\mathrm{g})2e^2\frac{(|k|^2/4)+6E(p)}{(E(pk)+\omega (k)E(p))^2}\frac{|\widehat{\phi }(k)|^2}{\omega (k)}𝑑k.$$

###### Proof.

By the pull-through formula

$$a(k,j)H_p\psi _\mathrm{g}=\left(H_{pk}+\omega (k)\right)a(k,j)\psi _\mathrm{g}$$

$$e\frac{\widehat{\phi }(k)}{\sqrt{2\omega (k)}}\left\{(pP_\mathrm{f}eA_\phi )e_j(k)+\frac{1}{2}\sigma (ike_j(k))\right\}\psi _\mathrm{g}.$$

Hence we have

$$a(k,j)\psi _\mathrm{g}=e\frac{\widehat{\phi }(k)}{\sqrt{2\omega (k)}}\left(H_{pk}+\omega (k)E(p)\right)^1$$

$$\times \left\{(pP_\mathrm{f}eA_\phi )e_j(k)+\frac{1}{2}\sigma (ike_j(k))\right\}\psi _\mathrm{g}.$$

Thus

$$(\psi _\mathrm{g},N_\mathrm{f}\psi _\mathrm{g})=\underset{j=1,2}{}a(k,j)\psi _\mathrm{g}^2𝑑k$$

$$\underset{j=1,2}{}\left|\frac{\widehat{\phi }(k)}{\sqrt{2\omega (k)}}\frac{(pP_\mathrm{f}eA_\phi )e_j(k)\psi _\mathrm{g}+(1/2)\sigma (ike_j(k))\psi _\mathrm{g}}{E(pk)+\omega (k)E(p)}\right|^2𝑑k.$$

We note that

$$(pP_\mathrm{f}eA_\phi )e_j(k)\psi _\mathrm{g}^2=\underset{\mu =1,2,3}{}\sigma _\mu (pP_\mathrm{f}eA_\phi )_\mu \sigma _\mu e_j^\mu (k)\psi _\mathrm{g}^2$$

$$3\underset{\mu =1,2,3}{}\sigma _\mu (pP_\mathrm{f}eA_\phi )_\mu \psi _\mathrm{g}^26E(p),$$

and

$$\sigma (ike_j(k))\psi _\mathrm{g}^2|k|^2\underset{\mu =1,2,3}{}(i\widehat{k}e_j(k))_\mu \psi _\mathrm{g}^2|k|^2.$$

Thus

$$(pP_\mathrm{f}eA_\phi )e_j(k)\psi _\mathrm{g}+(1/2)\sigma (ike_j(k))\psi _\mathrm{g}^22\left\{(|k|^2/4)+6E(p)\right\},$$

which leads to

$$(\psi _\mathrm{g},N_\mathrm{f}\psi _\mathrm{g})2e^2\frac{(|k|^2/4)+6E(p)}{(E(pk)+\omega (k)E(p))^2}\frac{|\widehat{\phi }(k)|^2}{\omega (k)}𝑑k$$

and the lemma follows. ∎We set

$$\theta (p)=\theta (p,e)=2\frac{(|k|^2/4)+6E(p)}{(E(pk)+\omega (k)E(p))^2}\frac{|\widehat{\phi }(k)|^2}{\omega (k)}𝑑k.$$

Note that $`\theta (p)`$ is rotation invariant, i.e. $`\theta (Rp)=\theta (p)`$ for an arbitrary rotation $`R`$. Let $`P_0=1P_\mathrm{\Omega }`$ be the projection onto $`C^2\{C\mathrm{\Omega }\}`$.

###### Lemma 2.2.

Let $`\mathrm{\Delta }(p)>0`$. Suppose $`|e|<1/\sqrt{3\theta (p)}`$. Then $`\mathrm{Tr}P_\mathrm{g}2`$.

###### Proof.

By Lemma 2.1 we have $`\mathrm{Tr}(P_\mathrm{g}N_\mathrm{f})e^2\theta (p)\mathrm{Tr}P_\mathrm{g}.`$ Therefore

$$\mathrm{Tr}P_\mathrm{g}\mathrm{Tr}(P_\mathrm{g}P_0)=\mathrm{Tr}P_\mathrm{g}(IP_0)\mathrm{Tr}(P_\mathrm{g}N_\mathrm{f})e^2\theta (p)\mathrm{Tr}P_\mathrm{g},$$

and

$$(1e^2\theta (p))\mathrm{Tr}P_\mathrm{g}\mathrm{Tr}(P_\mathrm{g}P_0)2,$$

which implies

$$\mathrm{Tr}P_\mathrm{g}\frac{2}{1e^2\theta (p)}<3.$$

Thus the lemma follows. ∎

### 2.2 Lower bound

We say that $`\psi `$ is real, if for all $`n0`$, $`\psi ^{(n)}`$ is a real-valued function on $`L^2(R^{3n}\times \{1,2\})`$. The set of real $`\psi `$ is denoted by $`_{\mathrm{real}}`$. We define the set of reality-preserving operators $`𝒪_{\mathrm{real}}`$ as follows:

$$𝒪_{\mathrm{real}}=\left\{A\right|A:_{\mathrm{real}}D(A)_{\mathrm{real}}\}.$$

It is seen that $`H_\mathrm{f}`$ and $`P_\mathrm{f}`$ are in $`𝒪_{\mathrm{real}}`$. Since for all $`kR`$ and $`zR`$ with $`z\mathrm{Spec}(H_{p0})`$,

$$((H_{p0}z)^k\psi )^{(n)}(k_1,j_1,\mathrm{},k_n,j_n)$$

$$=\left(\frac{1}{2}\left(p\underset{i=1}{\overset{n}{}}k_i\right)^2+\underset{i=1}{\overset{n}{}}\omega (k_i)+z\right)^k\psi ^{(n)}(k_1,j_1,\mathrm{},k_n,j_n),$$

$`(H_{p0}z)^k`$ is also in $`𝒪_{\mathrm{real}}`$. Moreover $`a(f)`$ and $`a^{}(f)`$ are in $`𝒪_{\mathrm{real}}`$ for real $`f`$’s. In particular $`A_\phi `$ and $`iB_\phi `$ are in $`𝒪_{\mathrm{real}}`$. Note that, if $`\psi D(B_\phi )_{\mathrm{real}}`$, then $`(\psi ,B_\phi \psi )=0.`$ Let $`_{\mathrm{fin}}=_{N=0}^{\mathrm{}}_{n=0}^N\left(L^2(R^3\times \{1,2\})\right)_{\mathrm{sym}}^n`$ denote the finite particle subspace of $``$. Note that $`A_\phi `$, $`B_\phi `$, $`\left(zH_{p0}\right)^{1/2}`$ and $`\left(zH_{p0}\right)^1`$ leave $`_{\mathrm{fin}}`$ invariant.

###### Lemma 2.3.

Suppose $`|e|<e^{}`$. Let $`xC^2`$. Then there exists $`a(t)R`$ independent of $`x`$ such that for $`t0`$

$$(x\mathrm{\Omega },e^{t(H_pE(p))}x\mathrm{\Omega })_{}=a(t)(x,x)_{C^2}.$$
(15)
###### Proof.

Note that $`H_{p\mathrm{I}}(1+H_{p0})^1<1`$ for $`|e|<e^{}`$ by (11). Then, by spectral theory, one has

$$e^{t(H_pE(p))}=\underset{n\mathrm{}}{lim}\left(1+\frac{t}{n}(H_pE(p))\right)^n$$

$$=\underset{n\mathrm{}}{lim}\underset{m\mathrm{}}{lim}\left\{\underset{k=0}{\overset{m}{}}\left(1+\frac{t}{n}H_{p0}\right)^1\left\{\left(\frac{t}{n}(H_{p\mathrm{I}}E(p))\right)\left(1+\frac{t}{n}H_{p0}\right)^1\right\}^k\right\}^n$$

$$=\underset{n\mathrm{}}{lim}\underset{m\mathrm{}}{lim}\left\{\left(1+\frac{t}{n}H_{p0}\right)^{1/2}\left(\underset{k=0}{\overset{m}{}}\left(\frac{t}{n}\stackrel{~}{H_\mathrm{I}}\right)^k\right)\left(1+\frac{t}{n}H_{p0}\right)^{1/2}\right\}^n.$$

Here

$$\stackrel{~}{H_\mathrm{I}}=\stackrel{~}{H}_{p\mathrm{II}}+i\sigma \stackrel{~}{B}_\phi ,$$

$$\stackrel{~}{H}_{p\mathrm{II}}=\left(1+\frac{t}{n}H_{p0}\right)^{1/2}(H_{p\mathrm{II}}E(p))\left(1+\frac{t}{n}H_{p0}\right)^{1/2},$$

$$\stackrel{~}{B}_\phi =\left(1+\frac{t}{n}H_{p0}\right)^{1/2}(iB_\phi )\left(1+\frac{t}{n}H_{p0}\right)^{1/2},$$

$$H_{p\mathrm{II}}=e(pP_\mathrm{f})A_\phi +\frac{e^2}{2}A_\phi ^2.$$

It is seen that

$$\stackrel{~}{H_\mathrm{I}}^2=\stackrel{~}{H}_{p\mathrm{II}}\stackrel{~}{H}_{p\mathrm{II}}\stackrel{~}{B}_\phi \stackrel{~}{B}_\phi +i\sigma (\stackrel{~}{H}_{p\mathrm{II}}\stackrel{~}{B}_\phi +\stackrel{~}{B}_\phi \stackrel{~}{H}_{p\mathrm{II}}\stackrel{~}{B}_\phi \stackrel{~}{B}_\phi )=M+i\sigma L.$$

Here both of $`M=\stackrel{~}{H}_{p\mathrm{II}}\stackrel{~}{H}_{p\mathrm{II}}\stackrel{~}{B}_\phi \stackrel{~}{B}_\phi `$ and $`L=\stackrel{~}{H}_{p\mathrm{II}}\stackrel{~}{B}_\phi +\stackrel{~}{B}_\phi \stackrel{~}{H}_{p\mathrm{II}}\stackrel{~}{B}_\phi \stackrel{~}{B}_\phi `$ are in $`𝒪_{\mathrm{real}}`$. Moreover

$$\stackrel{~}{H_\mathrm{I}}^3=\stackrel{~}{H}_{p\mathrm{II}}M\stackrel{~}{B}_\phi L+i\sigma (\stackrel{~}{B}_\phi M+\stackrel{~}{H}_{p\mathrm{II}}L\stackrel{~}{B}_\phi L),$$

where both of $`\stackrel{~}{H}_{p\mathrm{II}}M\stackrel{~}{B}_\phi L`$ and $`\stackrel{~}{B}_\phi M+\stackrel{~}{H}_{p\mathrm{II}}L\stackrel{~}{B}_\phi L`$ are also in $`𝒪_{\mathrm{real}}`$. Thus, repeating above procedure, one obtains

$$\underset{k=0}{\overset{m}{}}\left(\frac{t}{n}\stackrel{~}{H_\mathrm{I}}\right)^k=a_m+i\sigma b_m,$$

where $`a_m`$ and $`b_m`$ are in $`𝒪_{\mathrm{real}}`$. Hence there exist $`a_{nm}𝒪_{\mathrm{real}}`$ and $`b_{nm}𝒪_{\mathrm{real}}`$ such that

$$\left\{\left(1+\frac{t}{n}H_{p0}\right)^{1/2}\left(\underset{k=0}{\overset{m}{}}\left(\frac{t}{n}\stackrel{~}{H_\mathrm{I}}\right)^k\right)\left(1+\frac{t}{n}H_{p0}\right)^{1/2}\right\}^n=a_{nm}+i\sigma b_{nm}.$$

Finally

$$(x\mathrm{\Omega },e^{t(H_pE(p))}x\mathrm{\Omega })$$

$$=\underset{n\mathrm{}}{lim}\underset{m\mathrm{}}{lim}(x,x)(\mathrm{\Omega },a_{nm}\mathrm{\Omega })+i\underset{n\mathrm{}}{lim}\underset{m\mathrm{}}{lim}(x,\sigma x)(\mathrm{\Omega },b_{nm}\mathrm{\Omega }).$$

Since the left-hand side is real, the second term of the right-hand side vanishes and $`a(t)=lim_n\mathrm{}lim_m\mathrm{}(\mathrm{\Omega },a_{nm}\mathrm{\Omega })`$ exists, which establishes the desired result. ∎

###### Lemma 2.4.

Suppose that $`\mathrm{\Delta }(p)>0`$ and $`|e|<1/\sqrt{\theta (p)}`$. Then $`(\psi _\mathrm{g},P_0\psi _\mathrm{g})0.`$

###### Proof.

Since $`P_\mathrm{\Omega }+N_\mathrm{f}1`$, we have from Lemma 2.1

$$(\psi _\mathrm{g},P_0\psi _\mathrm{g})\psi _\mathrm{g}^2(1N_{\mathrm{f}}^{}{}_{}{}^{1/2})\psi _\mathrm{g}^2>1e^2\theta (p)>0.$$

Thus the lemma follows. ∎

###### Lemma 2.5.

Suppose $`|e|<e^{}`$ and $`|e|<1/\sqrt{\theta (p)}`$. Then there exists $`a>0`$ such that

$$P_0P_\mathrm{g}P_0=aP_0.$$
(16)
###### Proof.

Note that

$$P_\mathrm{g}=s\underset{t\mathrm{}}{lim}e^{t(H_pE(p))}.$$

Thus by Lemma 2.3,

$$(x\mathrm{\Omega },P_\mathrm{g}x\mathrm{\Omega })=\underset{t\mathrm{}}{lim}(x\mathrm{\Omega },e^{t(H_pE(p))}x\mathrm{\Omega })=\underset{t\mathrm{}}{lim}a(t)(x,x)$$

for all $`xC^2`$. Since by Lemma 2.4, $`(x\mathrm{\Omega },P_\mathrm{g}x\mathrm{\Omega })0`$ for some $`xC^2`$, $`lim_t\mathrm{}a(t)`$ exists and it does not vanish. For arbitrary $`\varphi _1,\varphi _2`$, the polarization identity leads to $`(\varphi _1,P_0P_\mathrm{g}P_0\varphi _2)=a(\varphi _1,P_0\varphi _2).`$ The lemma follows. ∎

###### Lemma 2.6.

Suppose $`|e|<e^{}`$ and $`|e|<1/\sqrt{\theta (p)}`$. Then $`\mathrm{Tr}P_\mathrm{g}2`$.

###### Proof.

Suppose $`\mathrm{Tr}P_\mathrm{g}=1`$. Then $`P_0P_\mathrm{g}P_0/\mathrm{Tr}P_0P_\mathrm{g}P_0`$ is a one-dimensional projection which contradicts (16). Thus the lemma follows. ∎

### 2.3 Proofs of Theorems 1.1 and 1.2

We define

$$e_0=inf\{|e|||e|<1/\sqrt{3\theta (p)},|e|<e^{}\}.$$

Proof of Theorem 1.1:
Since $`\mathrm{\Delta }(p)>0`$ and $`|e|<e_0`$, we conclude $`\mathrm{Tr}P_\mathrm{g}2`$ by Lemma 2.6 and $`\mathrm{Tr}P_\mathrm{g}2`$ by Lemma 2.2. Hence the theorem follows. ∎Proof of Theorem 1.2:
Without restriction in generality we may suppose $`\widehat{p}=(0,0,1)`$. Let $`\psi _\pm `$ be ground states of $`H_p`$ such that $`\psi _+(z)`$ and $`\psi _{}(z^{})`$ with some $`z,z^{}Z_{1/2}`$. $`\mathrm{\Omega }_+=\left(\begin{array}{c}\mathrm{\Omega }\\ 0\end{array}\right)`$ and $`\mathrm{\Omega }_{}=\left(\begin{array}{c}0\\ \mathrm{\Omega }\end{array}\right)`$ are ground states of $`H_{p0}`$ and $`\mathrm{\Omega }_\pm (\pm 1/2)`$. Let $`P_\mathrm{g}\mathrm{\Omega }_+=c_1\psi _++c_2\psi _{}`$ and $`P_\mathrm{g}\mathrm{\Omega }_{}=c_3\psi _++c_4\psi _{}`$ with some $`c_jC`$ and $`Q_{\pm 1/2}`$ be the projection of $``$ onto $`(\pm 1/2)`$. Since $`P_0P_\mathrm{g}P_0=aP_0`$ we conclude that

$$(\mathrm{\Omega }_+,P_\mathrm{g}\mathrm{\Omega }_+)=a>0,$$
(17)
$$(\mathrm{\Omega }_{},P_\mathrm{g}\mathrm{\Omega }_{})=a>0.$$
(18)
Then $`Q_{1/2}P_\mathrm{g}\mathrm{\Omega }_+0`$ and $`Q_{1/2}P_\mathrm{g}\mathrm{\Omega }_{}0`$. The alternative $`Q_{1/2}\psi _+0`$ or $`Q_{1/2}\psi _{}0`$ holds by (17), the alternative $`Q_{1/2}\psi _+0`$ or $`Q_{1/2}\psi _{}0`$ by (18). We may set $`Q_{1/2}\psi _+0`$. Then $`\psi _+(1/2)`$ and $`\psi _{}(1/2)`$. ∎

## 3 Zero total momentum

### 3.1 Spinless Hamiltonian

In the spinless case $`H_p`$ simplifies to

$$H_p=\frac{1}{2}(pP_\mathrm{f}eA_\phi )^2+H_\mathrm{f}$$

acting on $``$. The bound (12) is available and $`H_p`$ has at least one ground state for $`|p|<p_\mathrm{c}`$ with some $`p_\mathrm{c}>0`$ and arbitrary $`e`$. To show the uniqueness of the ground state only the pull-through argument seems to be available. The details of Section 2.1 remain unchanged and one concludes that if

$$|e|^2\frac{1}{2}\left\{\frac{E(p)}{(E(pk)+\omega (k)E(p))^2}\frac{|\widehat{\phi }(k)|^2}{\omega (k)}𝑑k\right\}^1,$$

then $`\mathrm{Tr}P_\mathrm{g}1`$, which implies that the ground state of $`H_p`$ is unique.

### 3.2 Confining potentials

For $`p=0`$ no infrared divergence is expected and for the remainder of this section we set

$$\omega (k)=|k|$$

as the physical dispersion relation. We did not succeed to apply the methods of Section 2 to this case. Therefore rather than considering $`p=0`$ directly we add to the Hamiltonian (3) a confining potential, $`V:R^3R`$, which in spirit amounts to the same physical situation. The Hamiltonian under study is

$$H_V=\frac{1}{2}\left\{\sigma \left(i_xeA_\phi (x)\right)\right\}^2+V(x)+H_\mathrm{f}.$$
(19)
Let

$$H_{\mathrm{el}}=\frac{1}{2}\mathrm{\Delta }+V$$

on $`L^2(R^3)C^2`$,

$$\mathrm{\Sigma }_{\mathrm{el}}=inf\mathrm{Spec}_{\mathrm{ess}}(H_{\mathrm{el}}),$$

$$E_{\mathrm{el}}=inf\mathrm{Spec}(H_{\mathrm{el}})$$

and

$$E=inf\mathrm{Spec}(H_V).$$

Let $`V`$ be relatively bounded with respect to $`\mathrm{\Delta }`$ with a bound less than 1. Then $`H_V`$ is self-adjoint on $`D(\mathrm{\Delta })D(H_\mathrm{f})`$ as established in . For arbitrary $`e`$ the existence of ground states has been proven by Griesemer, Lieb, and Loss  under the condition that $`H_{\mathrm{el}}`$ has a ground state separated by a gap from the continuous spectrum, i.e.

$$\mathrm{\Sigma }_{\mathrm{el}}E_{\mathrm{el}}>0.$$
(20)
We also refer to  for prior results, where in particular it is proven that the charge density of an arbitrary ground state $`\psi _\mathrm{g}`$ is localized, i.e.

$$e^{c|x|}\psi _\mathrm{g}c_1\psi _\mathrm{g}$$
(21)
with some constant $`c`$. In the spinless case the uniqueness of the ground state, $`\mathrm{Tr}P_\mathrm{g}=1`$, would follow from a positivity argument . Let $`P_{\mathrm{el}}`$ be the projection of the subspace spanned by ground states of $`H_{\mathrm{el}}`$. Suppose (20). Take $`e`$ such that

$$\mathrm{\Sigma }_{\mathrm{el}}E>0,$$
(22)
which can be satisfied by the continuity of $`E`$ in $`e`$. Then a pull-through argument and (21) yield that

$$(\psi _\mathrm{g},(1N_\mathrm{f}+P_{\mathrm{el}}^{}P_0)\psi _\mathrm{g})\eta (e),$$
(23)
where $`lim_{ϵ0}\eta (e)=0`$. Hence in the similar way as Lemma 2.2 with $`P_0`$ and $`N_\mathrm{f}`$ replaced by $`P_{\mathrm{el}}P_0`$ and $`1N_\mathrm{f}+P_{\mathrm{el}}^{}P_0`$, respectively, we see that if, in addition to (22), $`e`$ satisfies $`\eta (e)<1/3`$, then

$$\mathrm{Tr}P_\mathrm{g}<3.$$
(24)
The realness argument of Section 2.2 requires some extra conditions on $`V`$.

###### Theorem 3.1.

Suppose $`V(x)=V(x)`$, $`\mathrm{\Sigma }_{\mathrm{el}}E_{\mathrm{el}}>0`$, and $`\mathrm{Tr}P_{\mathrm{el}}=2`$. Then there exists a positive constant $`e_{00}`$ such that, if $`|e|<e_{00}`$, then $`H_V`$ has a two-fold degenerate ground state.

###### Proof.

Suppose that $`e`$ satisfies (22) and $`\eta (e)<1`$. (23) yields

$$(\psi _\mathrm{g},P_{\mathrm{el}}P_0\psi _\mathrm{g})\psi _\mathrm{g}^2(\psi _\mathrm{g},(1N_\mathrm{f}+P_{\mathrm{el}}^{}P_0)\psi _\mathrm{g})1\eta (e)>0.$$
(25)
Let $`F`$ denote the Fourier transform of $`L^2(R^3)`$ and define the unitary operator of $``$ by $`T=Fe^{ixP_\mathrm{f}}`$. Then we have

$$TH_VT^1=\frac{1}{2}\left\{\sigma (xP_\mathrm{f}eA_\phi (0))\right\}^2+FVF^1+H_\mathrm{f}.$$

The assumption $`V(x)=V(x)`$ implies that $`FVF^1`$ is a reality preserving operator on $`L^2(R^3)`$. Let $`\phi _{\mathrm{el}}`$ be the ground state of $`H_{\mathrm{el}}`$. In the similar way as Lemmas 2.3 and 2.5 with $`\mathrm{\Omega }`$ and $`P_0`$ replaced by $`\phi _{\mathrm{el}}\mathrm{\Omega }`$ and $`P_{\mathrm{el}}P_0`$, respectively, it is established that by (25) there exists a positive constant $`c^{}`$ such that if $`|e|<c^{}`$, then (22) and $`\eta (e)<1`$ hold, and $`\mathrm{Tr}P_\mathrm{g}2`$. Take $`e_{00}=sup\left\{|e||\eta (e)<1/3,|e|<c^{}\right\}.`$ Then the theorem follows from (24). ∎

Acknowledgment

We thank Volker Bach for explaining to us that the overlap with the vacuum yields an upper bound on the degeneracy. F.H. gratefully acknowledges the kind hospitality at Technische Universität München. This work is in part supported by the Graduiertenkolleg “Mathematik in ihrer Wechselbeziehung zur Physik” of the LMU München and Grant-in-Aid 13740106 for Encouragement of Young Scientists from the Ministry of Education, Science, Sports, and Culture.