format: text

text:
# Towards Order 𝛼_𝑠⁴ Accuracy in 𝜏-decays

## I Introduction

The dependence of the $`\tau `$-decay rate on the strong coupling $`\alpha _s`$ has been used for a determination of $`\alpha _s`$ at lower energies, with the most recent results of $`0.334\pm 0.007_{\mathrm{exp}}\pm 0.021_{\mathrm{theo}}`$ and $`0.3478\pm 0.009_{\mathrm{exp}}\pm 0.019_{\mathrm{theo}}`$ by the ALEPH exp1  and OPAL exp2  collaborations. After evolution up to higher energies these results agree remarkably well with determinations based on the hadronic $`Z`$ decay rate. In view of the relatively large value of $`\alpha _s(M_\tau )`$ estimates for the yet unknown terms of higher orders play an important role in current determinations of $`\alpha _s`$ at low energies. This is in contrast to high energy measurements, where the uncertainty from terms of order $`\alpha _s^4`$ of 0.001 to 0.002 phys\_report  is somewhat less or at most comparable to the present experimental errors.

The situation is even more problematic for the determination of the strange quark mass from the Cabbibo suppressed $`\tau `$ decays. Perturbative QCD corrections affecting the $`m_s^2`$ term are extremely large and contributions from increasing powers of $`\alpha _s`$ are barely decreasing, which casts doubts on our ability to extract a reliable result for $`m_s`$ from this (in principle) clean and straightforward measurement Maltman:1998qz ; Chetyrkin93 ; Chetyrkin93-update ; Chetyrkin:1998ej ; pichprades ; krajms ; chenpich ; Maltman:2002wb .

Partial results of order $`\alpha _s^4`$ for the absorptive part of the massless vector and scalar correlators have been obtained recently Baikov:2001aa , namely terms proportional to $`n_f^2`$, where $`n_f`$ denotes the number of massless fermion species. These allow to test two popular optimization schemes — based on the principles of “minimal sensitivity” (PMS) and of fastest apparent convergence (FAC) ste ; gru ; krasn  — which have been used to predict yet uncalulated higher order terms KS ; Kataev:1994wd .

It will be demonstrated that the predictions of both schemes (coinciding at $`\alpha _s^4`$) for the coefficient of order $`n_f^2\alpha _s^4`$ are in reasonable agreement with our calculations, which are then used to predict the complete fixed order (FO) and the ”contour improved” (CI) tau:resum ; DP  $`𝒪(\alpha _s^4)`$ contributions to the $`\tau `$ decay rate. Employing the four-loop QCD beta function in combination with improved $`\alpha _s^4`$ terms a rough estimate even for $`𝒪(\alpha _s^5)`$ terms can be obtained. The results lead to fairly stable values for $`\alpha _s`$ consistent with current analysises.

Unfortunately, no essential decrease of the difference between the central values $`\alpha _s(M_\tau )`$ as obtained with FI and CI approaches is observed after including order $`\alpha _s^4`$ and $`\alpha _s^5`$ corrections, assuming for the moment that these estimates are indeed correct. On the other hand, the theoretical uncertainty assigned to $`\alpha _s(M_\tau )`$ within each method according to standard techniques does decrease significantly.

The implication of this approach for the extraction of $`m_s`$ from Cabbibo suppressed decays will be investigated along the same lines. New results will be presented for the terms of order $`n_f\alpha _s^3m_s^2`$ in the total rate. In this case the agreement between FAC/PMS predictions and our results is quite encouraging and naturally suggests the use of the former as a reliable prediction for the complete $`\alpha _s^3m_s^2`$ term. Following an approach discussed in Kataev:1994wd , even a rough estimate $`m_s^2\alpha _s^4`$ terms can be obtained from these considerations.

However, the rapid increase of the coefficients indicates that the inherent uncertainty of the present $`m_s`$ determinations will not necessarily decrease with inclusion of the higher orders. As we will see, the situation is somewhat better for the spin one contribution if considered separately.

## II Generalities

We start with the well-known representation SchTra84 ; Bra88 ; NarPic88 ; Bra89 ; BraNarPic92  of the tau-lepton hadronic rate as the contour integral along a circle C of radius $`|s|=M_\tau ^2`$

$$R_\tau =6i\pi _{|s|=M_\tau ^2}\frac{ds}{M_\tau ^2}\left(1\frac{s}{M_\tau ^2}\right)^2\left[\mathrm{\Pi }^{[2]}(s)\frac{2}{M_\tau ^2}\mathrm{\Pi }^{[1]}(s)\right].$$
(1)
Here $`\mathrm{\Pi }^{[1]}`$ and $`\mathrm{\Pi }^{[2]}`$ are proper flavour combinations of the polarization operators appearing in the decomposition of the correlators of vector and axial vector currents of light quarks

$$\begin{array}{cc}\mathrm{\Pi }_{\mu \nu ,ij}^{V/A}(q,m_i,m_j,m,\mu ,\alpha _s)\hfill & =i𝑑xe^{iqx}T[j_{\mu ,ij}^{V/A}(x)(j_{\nu ,ij}^{V/A})^{}(0)]\hfill \\ & =g_{\mu \nu }\mathrm{\Pi }_{ij,V/A}^{[1]}(q^2)+q_\mu q_\nu \mathrm{\Pi }_{ij,V/A}^{[2]}(q^2)\hfill \end{array}$$
(2)
with $`m^2=_{f=u,d,s}m_f^2`$ and $`j_{\mu ,ij}^{V/A}=\overline{q}_i\gamma _\mu (\gamma _5)q_j`$. The two (generically different) quarks with masses $`m_i`$ and $`m_j`$ are denoted by $`q_i`$ and $`q_j`$ respectively.

For the case of the $`\tau `$-lepton the relevant combinations of quark flavours are $`ij=ud`$ and $`ij=us`$. The polarization functions $`\mathrm{\Pi }_{V/A}^{[l]},l=1,2`$ are conveniently represented in the form ($`Q^2=q^2`$)

$$(Q^2)^{(l2)}\mathrm{\Pi }_{us,V/A}^{[l]}(q^2)=\frac{3}{16\pi ^2}\mathrm{\Pi }_{V/A,0}^{[l]}(\frac{\mu ^2}{Q^2},\alpha _s)+\frac{3}{16\pi ^2}\underset{D2}{}Q^D\mathrm{\Pi }_{V/A,D}^{[l]}(\frac{\mu ^2}{Q^2},m_s^2,\alpha _s).$$
(3)
Here the first term on the rhs corresponds to the massless limit while the first term in the sum stands for quadratic mass corrections. We neglect the masses of $`u`$ and $`d`$ quarks. Therefore in perturbative QCD $`\mathrm{\Pi }_V^{[l]}=\mathrm{\Pi }_A^{[l]}`$ and we will often omit the subscript $`V/A`$ in the following. Current conservation implies: $`\mathrm{\Pi }_0^{[1]}=\mathrm{\Pi }_0^{[2]}`$.

The full tau-lepton hadron rate $`R_\tau `$ can be presented as a sum of spin 1 and spin 0 parts, viz.

$$\begin{array}{cc}R_\tau ^{(1)}\hfill & =6i\pi _{|s|=M_\tau ^2}\frac{ds}{M_\tau ^2}\left(1\frac{s}{M_\tau ^2}\right)^2\left[\left(1+2\frac{s}{M_\tau ^2}\right)\mathrm{\Pi }^{(1)}(s)+\mathrm{\Pi }^{[1]}(0)/s\right],\hfill \\ R_\tau ^{(0)}\hfill & =6i\pi _{|s|=M_\tau ^2}\frac{ds}{M_\tau ^2}\left(1\frac{s}{M_\tau ^2}\right)^2\left[\mathrm{\Pi }^{(0)}(s)\mathrm{\Pi }^{[1]}(0)/s\right].\hfill \end{array}$$
(4)
where

$$\mathrm{\Pi }^{(1)}=\mathrm{\Pi }^{[1]}/q^2,\mathrm{\Pi }^{(0)}=\mathrm{\Pi }^{[2]}+\mathrm{\Pi }^{[1]}/q^2.$$
(5)
and the contribution of the singularity at the origin (proportional to $`\mathrm{\Pi }^{[1]}(0)`$) has to be included. A nonvanishing value of $`\mathrm{\Pi }^{[1]}(0)`$ is a nonperturbative constant.

On the other hand, the unknown constant drops out if one considers moments

$$R_\tau ^{(1,0)k,l}(s_0)=_0^{s_0}\frac{ds}{M_\tau ^2}\left(1\frac{s}{M_\tau ^2}\right)^k\left(\frac{s}{M_\tau ^2}\right)^l\frac{dR_\tau ^{(1,0)}}{ds},$$
(6)
with $`k0,l1`$. (Note that the moments introduced in DP  are related to ours as $`R_\tau ^{kl}=R_\tau ^{(1)k,l}+R_\tau ^{(0)k,l}`$.)

The decay rate $`R_\tau `$ may be expressed as the sum of different contributions corresponding to Cabibbo suppressed or allowed decay modes, vector or axial vector contributions and the mass dimension of the corrections

$$R_\tau =R_{\tau ,V}+R_{\tau ,A}+R_{\tau ,S}$$
(7)
with

$$\begin{array}{cc}R_V=\hfill & \frac{3}{2}|V_{ud}|^2\left(1+\delta _0+\underset{D=2,4,\mathrm{}}{}\delta _{V,ud,D}\right),\hfill \\ R_A=\hfill & \frac{3}{2}|V_{ud}|^2\left(1+\delta _0+\underset{D=2,4,\mathrm{}}{}\delta _{A,ud,D}\right),\hfill \\ R_S=\hfill & 3|V_{us}|^2\left(1+\delta _0+\underset{D=2,4,\mathrm{}}{}\delta _{us,D}\right).\hfill \end{array}$$
(8)
Here $`D`$ indicates the mass dimension of the fractional corrections $`\delta _{V/A,ij,D,}`$ and $`\delta _{ij,D}`$ denotes the average of the vector and the axial vector contributions: $`\delta _{ij,D}=(\delta _{V,ij,D}+\delta _{A,ij,D})/2`$. If a decomposition into different spin/parity contributions is made or a particular pattern of moments is considered then we will use the corresponding obvious generalization of (8). For instance,

$$R_{S,V}^{(1)kl}=a_{kl}|V_{us}|^2\left(1+\delta _0^{kl}+\underset{D=2,4,\mathrm{}}{}\delta _{V,us,D}^{(1),kl}\right),$$
(9)
and

$$R_{S,V}^{(0)kl}=|V_{us}|^2\left(\underset{D=2,4,\mathrm{}}{}\delta _{V,us,D}^{(0),kl}\right).$$
(10)
Thus, in our notation we have the relation

$$\delta _{V,us,2}^{kl}=a_{kl}\delta _{V,us,2}^{(1),kl}+\delta _{V,us,2}^{(0),kl}.$$
(11)
The integral in eq. (1) is, obviously, insensitive, to the $`Q^2`$-independent terms in the polarization functions $`\mathrm{\Pi }_0^{[1]}`$ and $`\mathrm{\Pi }_2^{[1]}`$. This means that without loss of generality we may deal with the corresponding (Adler) $`D`$-functions, viz.

$$D_0^{[1]}(Q^2)\frac{3}{4}Q^2\frac{\mathrm{d}}{\mathrm{d}Q^2}\mathrm{\Pi }_0^{[1]},D_2^{[1]}(Q^2)\frac{1}{2}Q^2\frac{\mathrm{d}}{\mathrm{d}Q^2}\mathrm{\Pi }_2^{[1]}.$$

An important property of the functions $`D_0^{[1]},D_2^{[1]}`$ and $`\mathrm{\Pi }_2^{[2]}`$ is their scale independence, which implies that they are directly related to measurements.

The Adler functions $`D_0^{[1]}`$ and $`D_2^{[1]}`$ have been calculated with $`𝒪(\alpha _s^3)`$ accuracy, the polarization function $`\mathrm{\Pi }_2^{[2]}`$, however, to $`𝒪(\alpha _s^2)`$ only (see Chetyrkin:1998ej  and references therein).

The (apparent) convergence of the perturbative series for $`D_0^{[1]}`$ is acceptable, the one for $`D_2^{[1]}`$ and $`\mathrm{\Pi }_2^{[2]}`$ is at best marginal. This has led to significant theoretical uncertainties in extracting $`\alpha _s`$ and to a fairly unstable behaviour in extraction of $`m_s`$ from $`\tau `$-decays.

To improve the situation, we have computed the two leading terms in the large $`n_f`$ expansion of the next order, i.e. terms of order $`𝒪(n_f^3\alpha _s^4)`$, $`𝒪(n_f^2\alpha _s^4)`$ to $`D_0^{[1]}`$, $`D_2^{[1]}`$ and $`𝒪(n_f^2\alpha _s^3)`$, $`𝒪(n_f\alpha _s^3)`$ to $`\mathrm{\Pi }_2^{[2]}`$. Our results are described in the following section.

## III Fixed order results in $`\alpha _s^3`$ and $`\alpha _s^4`$

Using the technique described in bai1 ; bai2 ; Baikov:2001aa  and the parallel version of FORM Fliegner:2000uy ; Vermaseren , the leading and subleading (in $`n_f`$) terms of the next order in the perturbative series for $`D_0^{[1]}`$, $`D_2^{[1]}`$ and $`\mathrm{\Pi }_2^{[2]}`$ have been obtained in the standard $`\overline{\text{MS}}`$ renormalization scheme ms ; msbar :

$`D_0^{[1]}`$ $`=`$ $`1+a_s+a_s^2\left\{\left[{\displaystyle \frac{11}{12}}+{\displaystyle \frac{2}{3}}\zeta _3\right]n_f+{\displaystyle \frac{365}{24}}11\zeta _3\right\}`$ (12)
$`+`$ $`a_s^3\{[{\displaystyle \frac{151}{162}}{\displaystyle \frac{19}{27}}\zeta _3]n_f^2+[{\displaystyle \frac{7847}{216}}+{\displaystyle \frac{262}{9}}\zeta _3{\displaystyle \frac{25}{9}}\zeta _5]n_f`$
$`+{\displaystyle \frac{87029}{288}}{\displaystyle \frac{1103}{4}}\zeta _3+{\displaystyle \frac{275}{6}}\zeta _5\}`$
$`+`$ $`a_s^4\{[{\displaystyle \frac{6131}{5832}}+{\displaystyle \frac{203}{324}}\zeta _3+{\displaystyle \frac{5}{18}}\zeta _5]n_f^3+[{\displaystyle \frac{1045381}{15552}}+{\displaystyle \frac{5}{6}}\zeta _3^2`$
$`{\displaystyle \frac{40655}{864}}\zeta _3{\displaystyle \frac{260}{27}}\zeta _5]n_f^2+d_{0,1}^{[1]4}n_f+d_{0,0}^{[1]4}\}`$
$`=`$ $`1+a_s+a_s^2\left\{0.1153n_f+1.986\right\}`$
$`+`$ $`a_s^3\left\{0.08621n_f^24.216n_f+18.24\right\}`$
$`+`$ $`a_s^4\left\{0.01009n_f^3+1.875n_f^2+d_{0,1}^{[1]4}n_f+d_{0,0}^{[1]4}\right\},`$
$`D_2^{[1]}=m_s^2(\text{}1`$ $`+`$ $`{\displaystyle \frac{5}{3}}a_s+a_s^2\left\{\left[{\displaystyle \frac{11}{8}}+{\displaystyle \frac{2}{3}}\zeta _3\right]n_f+{\displaystyle \frac{5185}{144}}{\displaystyle \frac{39}{2}}\zeta _3\right\}`$
$`+`$ $`a_s^3\{[{\displaystyle \frac{8671}{11664}}{\displaystyle \frac{13}{27}}\zeta _3]n_f^2+[{\displaystyle \frac{44273}{972}}+{\displaystyle \frac{3257}{81}}\zeta _3`$
$`{\displaystyle \frac{5}{6}}\zeta _4{\displaystyle \frac{1265}{81}}\zeta _5]n_f+{\displaystyle \frac{2641517}{5184}}{\displaystyle \frac{131275}{216}}\zeta _3+{\displaystyle \frac{12845}{36}}\zeta _5\}`$
$`+`$ $`a_s^4\{[{\displaystyle \frac{396781}{559872}}+{\displaystyle \frac{461}{1296}}\zeta _3{\displaystyle \frac{1}{48}}\zeta _4+{\displaystyle \frac{5}{18}}\zeta _5]n_f^3+[{\displaystyle \frac{61913567}{1119744}}{\displaystyle \frac{59}{54}}\zeta _3^2`$
$`{\displaystyle \frac{352549}{7776}}\zeta _3+{\displaystyle \frac{67}{96}}\zeta _4+{\displaystyle \frac{22859}{3888}}\zeta _5]n_f^2+d_{2,1}^{[1]4}n_f+d_{2,0}^{[1]4}\}\text{})`$
$`=m_s^2(\text{}1`$ $`+`$ $`1.667a_s+a_s^2\left\{0.5736n_f+12.57\right\}`$ (13)
$`+`$ $`a_s^3\{0.1646n_f^214.31n_f+149.\}`$
$`+`$ $`a_s^4\{0.01563n_f^3+6.067n_f^2+d_{2,1}^{[1]4}n_f+d_{2,0}^{[1]4}\}\text{}),`$
$`\mathrm{\Pi }_2^{[2]}=4m_s^2(\text{}1`$ $`+`$ $`{\displaystyle \frac{7}{3}}a_s+a_s^2\left\{\left[{\displaystyle \frac{25}{24}}{\displaystyle \frac{2}{9}}\zeta _3\right]n_f+{\displaystyle \frac{15331}{432}}+{\displaystyle \frac{359}{54}}\zeta _3{\displaystyle \frac{520}{27}}\zeta _5\right\}`$
$`+`$ $`a_s^3\{[{\displaystyle \frac{2131}{11664}}+{\displaystyle \frac{19}{81}}\zeta _3]n_f^2+[{\displaystyle \frac{68135}{1944}}{\displaystyle \frac{52}{27}}\zeta _3^2`$
$`{\displaystyle \frac{3997}{486}}\zeta _3{\displaystyle \frac{5}{6}}\zeta _4+{\displaystyle \frac{3875}{243}}\zeta _5]n_f+k_{2,0}^{[2]3}\}\text{})`$
$`=4m_s^2(\text{}1`$ $`+`$ $`2.333a_s+a_s^2\left\{1.309n_f+23.51\right\}`$ (14)
$`+`$ $`a_s^3\{0.4647n_f^232.08n_f+k_{2,0}^{[2]3}\}\text{}).`$
Here we have used $`a_s=\frac{\alpha _s(Q^2)}{\pi },m_s=m_s(Q^2)`$ and set the normalization scale $`\mu ^2=Q^2`$; results for generic values of $`\mu `$ can be easily recovered with the standard renormalization group techniques. The result for the $`\alpha _s^4`$ terms in $`D_0^{[1]}`$ has been already presented in Baikov:2001aa , the coefficients of the $`\alpha _s^4`$ and $`\alpha _s^3`$ terms in in $`D_2^{[1]}`$ and $`\mathrm{\Pi }_2^{[2]}`$ respectively are new.

## IV Implications for the $`\alpha _s^4`$ predictions and phenomenological analysis

### IV.1 Massless case

FAC (Fastest Apparent Convergence) and PMS (Principle of Minimal Sensitivity) methods are both based eventually on the concept of scheme-invariant properties and the idea of the choice of an ”optimal” scheme to provide better convergence of the resulting perturbative series. For both methods the optimal scheme depends on the physical observable we are dealing with. With FAC it should be a scheme which minimizes (set to zero by construction) all the terms of order $`\alpha _s^2`$ and higher, while the PMS scheme is fixed by the requirement that the perturbative expansion for the observable is as insensitive as possible to a change in the scheme. The assumption that a renormalization scheme is in a sense optimal sets certain constraints on not yet computed higher order corrections in any other scheme. These constraints can be used to ”predict” (at least roughly) the magnitude of these corrections.

For the function $`D_0^{[1]}`$ the result is known since long from Ref. KS  (see Table 1, column 5). From the three entries corresponding to $`n_f=3,4`$ and $`5`$ one easily restores the FAC/PMS prediction for the $`n_f`$ dependence of the $`\alpha _s^4`$ term in $`D_0^{[1]}`$ (the term of order $`n_f^3`$ was fixed to its computed value):

$$d_0^{[1]4}=127.644.2n_f+3.64n_f^20.0100928n_f^3d_0^{[1]4}=127.5844.211n_f+3.6439n_f^20.0100928n_f^3.$$
(15)
(Note that FAC and PMS predictions happen to coincide for the $`\alpha _s^4`$ term.)

It is interesting to compare the FAC/PMS predictions for the $`n_f`$ dependence of coefficient $`d_0^{[1]3}`$ with the exact result given in eq. 12. The results of both estimates are

$$d_0^{[1]3}(FAC)=8.541.013n_f+0.0116n_f^2,$$
(16)
$$d_0^{[1]3}(PMS)=9.931.23n_f+0.0125n_f^2.$$
(17)
The comparison of the complete $`\alpha _s^3`$ and partial $`\alpha _s^4`$ results with FAC and PMS estimates leads to the following observations:

* Starting from $`\alpha _s^3`$, the leading in $`\alpha _s`$ and $`n_f`$ terms of order $`\alpha _s^3n_f^2`$ and $`\alpha _s^4n_f^3`$ are numerically quite small (at least for $`n_f6`$) and, thus, should have a negligibly small influence on the coefficients of the $`\alpha _s`$ expansion. On the other hand, the term subleading in $`n_f`$, say, of order $`\alpha _s^3n_f`$ is comparable in size with the term of order $`\alpha _s^3n_f^0`$. Similarly, the $`\alpha _s^4n_f^3`$ term is significantly smaller than the $`\alpha _s^4n_f^2`$ one, whereas the $`\alpha _s^4n_f`$ and $`\alpha _s^4n_f^0`$ terms are expected to be of similar magnitude.
* In general FAC/PMS methods correctly reproduce sign and order of magnitude of the higher order coefficients. The agreement is getting better when the coefficients happen to be large, as is the case for the quadratic mass corrections (see below).
* Taken separately, the FAC/PMS estimates of the coefficients of the $`n_f`$ expansion could deviate rather strongly from the true result. However, for a given $`n_f`$, the deviation in the predicted value of the full $`O(\alpha _s^n)`$ term tends to be significantly smaller than what could be expected from summing individual terms of the $`n_f`$ expansion. In addition, for the particular point $`n_f=3`$ very good agreement is observed.

To illustrate this feature, let us consider the worst case: the FAC/PMS prediction for the $`\alpha _s^3n_f`$ term in the function $`d_0^{[1]}`$. Here the ratio of the exact result relative to the predicted one is quite large (about 4). Without knowledge of the $`\alpha _s^3n_f^0`$ contribution one would expect that the uncertainty of the prediction for the full $`\alpha _s^3`$ coefficient should be at least around

$$(d_{0,1}^{[1]3}|_{exact}d_{0,1}^{[1]3}|_{FAC})n_f.$$

For $`n_f=3,4,5`$ this amounts to 9, 12 and 15, which should be compared to the corresponding differences of the full order $`O(\alpha _s^3)`$ (summed over all contributing power of $`n_f`$) coefficients, viz. $`1,2,4`$. This example demonstrates that the deviation of FAC/PMS predictions for subleading in $`n_f`$ terms may well serve as conservative estimate of the accuracy in the prediction of the complete terms of orders $`\alpha _s^3`$ and $`\alpha _s^4`$.

These observations motivate the assumption that the prediction for the coefficient $`d_0^{[1]4}(d_{0,0}^{[1]4}+d_{0,1}^{[1]4}n_f+d_{0,2}^{[1]4}n_f^2+d_{0,3}^{[1]4}n_f^3)`$ should also be correct within

$$\pm (d_0^{[1]4}|_{exact}d_0^{[1]4}|_{FAC})n_f^2$$

Thus, in our phenomenological analysis of the $`\tau `$-lepton decays the estimate

$$d_0^{[1]4}|_{n_f=3}=27\pm 16$$
(18)
will be used. On the basis of this improved estimate and the four loop $`\beta `$ function vanRitbergen:1997va  one may even speculate about the $`\alpha _s^5`$ term (whose exact evaluation is completely out of reach in the foreseeable future). Following the discussion of Kataev and Starshenko Kataev:1994wd , one obtains

$$d_0^{[1]5}|_{n_f=3}=145\pm 100,$$
(19)
not far from the previous estimates of Ref. Kataev:1994wd . The variation of $`d_0^{[1]4}`$ by $`\pm 16`$ leads to the variation of $`d_0^{[1]5}`$ by $`\pm 100`$. For other values of $`n_f`$, the corresponding predictions can be obtained in the same way. They are listed in Table 1.

The FAC/PMS prediction (15) for the $`n_f`$ dependence of the coefficient $`d_0^{[1]4}`$ does not take into account the available knowledge of the corresponding $`n_f^2`$ part. One can easily include this by fitting the FAC/PMS predictions for only $`\mathrm{𝑡𝑤𝑜}`$ values of $`n_f`$ with a $`\mathrm{𝑙𝑖𝑛𝑒𝑎𝑟}`$ function of $`n_f`$. As a result one obtains 111We thank Matthias Steihnauser and Robert Harlander for a useful discusson of this point. (we have boxed the predicted coefficients in order to separate them clearly from the input)

$`d_0^{[1]4}(\mathrm{FAC}/\mathrm{PMS},n_f=3,4)`$ $`=`$ $`\overline{)105.731.8n_f}+1.875n_f^20.01009n_f^3,`$ (20)
$`d_0^{[1]4}(\mathrm{FAC}/\mathrm{PMS},n_f=4,5)`$ $`=`$ $`\overline{)107.732.3n_f}+1.875n_f^20.01009n_f^3,`$ (21)
$`d_0^{[1]4}(\mathrm{FAC}/\mathrm{PMS},n_f=3,5)`$ $`=`$ $`\overline{)106.432.0n_f}+1.875n_f^20.01009n_f^3.`$ (22)
One could now perform a self-consistency check of Table I by predicting, say, $`d_0^{[1]4}`$ for $`n_f=5`$ from eq. (20). The result – (-7.5) – is compared successfully to the value listed in the table, viz. -8. The corresponding predictions from eqs. (21) and (22) are also in very good agreement to Table I.

An instructive example of how knowledge and inclusion of the subleading $`n_f`$ term can improve FAC/PMS predictions is provided by the (exactly known) coefficient $`d_0^{[1]3}`$. Indeed, assuming the knowledge of $`d_{0,2}^{[1]3}`$ and $`d_{0,1}^{[1]3}`$ and using the values of $`d_3^{\mathrm{FAC}}`$ and $`d_3^{\mathrm{PMS}}`$ (at $`n_f=3`$) as given by Table 1, one easily arrives at

$$d_0^{[1]3}(\mathrm{FAC})=\text{17.48}4.216n_f+0.08621n_f^2$$
(23)
and

$$d_0^{[1]3}(\mathrm{PMS})=\text{18.27}4.216n_f+0.08621n_f^2,$$
(24)
which should be compared to the exact value $`d_{0,0}^{[1]3}=18.24`$. Repeating the same exercise for $`n_f=4,5`$ we get 20.16 (FAC) and 20.75 (PMS) for $`n_f=4`$ as well as $`22.69`$ (FAC) and 23.08 (PMS) for $`n_f=5`$ respectively.

Bearing in mind that in many cases FAC/PMS predictions made for $`n_f=3`$ are in better agreement with the exact results (see, e.g. Table I and tables from Ref. Kataev:1994wd ), we suggest eq. (20) as the best FAC/PMS prediction for the constant term and the term linear in $`n_f`$ in the coefficient $`d_0^{[1]4}`$.

Eqs. (18,19) can be used to predict $`R_\tau `$ in the massless limit first in fixed order perturbation theory (FOPT)

$`R_\tau ^{\mathrm{FOPT}}=`$
$`3\left(1+a_s+5.202a_s^2+26.37a_s^3+a_s^4(105\pm 16)+a_s^5(138\pm 230\pm 100)\right)`$ .
Here $`a_s=\alpha _s(M_\tau )/\pi `$. The first uncertainty in the $`\alpha _s^5`$ term comes from that of $`d_0^{[1]4}`$ while the second is our estimation of the error in the very coefficient $`d_0^{[1]5}`$. (Of course, within this approach they are strongly correlated).

Similarly, we can use ”contour improved” (CI) formulae tau:resum ; DP  (assuming as reference value $`\alpha _s(M_\tau )=0.334`$ exp1 ) to get

$$R_\tau ^{CI}=3\left(1+1.364a_s+2.54a_s^2+9.71a_s^3+1.31a_s^4d_0^{[1]4}+0.95a_s^5d_0^{[1]5}\right)$$
(26)
or, equivalently,

$$R_\tau ^{CI}=3\left(1+1.364a_s+2.54a_s^2+9.71a_s^3+a_s^4(35\pm 20)+a_s^5(138\pm 95)\right).$$
(27)
Let us compare our new value for the coefficient $`d_0^{[1]4}`$ in (18) with the ones used in extracting $`\alpha _s(M_\tau )`$ from $`\tau `$ data by the OPAL exp2  and ALEPH exp1  collaborations, namely

$$d_0^{[1]4}|_{n_f=3}=25\pm 25\text{(OPAL)},50\pm 50\text{(ALEPH)}.$$
(28)
While OPAL’s central value is basically the same as ours, their error bar is somewhat larger. In the case of ALEPH both the central value and the uncertainty assigned to it are significantly larger than our numbers. In this connection, we would like to stress that eq. (18) utilizes completely new non-trivial information given in (12): the subleading term in $`n_f`$ term of order $`\alpha _s^4`$.

It is of interest to see in detail which accuracy in the determination of $`\alpha _s(M_\tau )`$ one could achieve assuming eqs. (18,19). Let us introduce the quality $`\delta _P`$ as follows:

$$R_{\tau S=0}=\frac{\mathrm{\Gamma }(\tau h_{S=0}\nu )}{\mathrm{\Gamma }(\tau l\overline{\nu }\nu )}=|V_{ud}|^2S_{EW}R_\tau ,$$
(29)
with

$$R_\tau =3(1+\delta _P+\delta _{EW}+\delta _{NP}).$$

The first term here is the parton result, the second stands for pQCD effects. The non-perturbative corrections represented by $`\delta _{NP}`$ happens to be rather small $`\delta _{NP}=.003\pm 0.003`$ (see, e.g. NarPic88 ). Here the flavour mixing matrix element $`|V_{ud}|^2=0.9475\pm 0.0016`$ PDG . The factor $`S_{EW}=1.0194`$ is the electroweak correction which collects the large logarithmic terms ewcorr1 , while $`\delta _{EW}=0.001`$ is an additive electroweak correction ewcorr2 . Using for definiteness the result of ALEPH exp1

$$R_{\tau S=0}=3.492\pm 0.016,$$
(30)
one arrives at

$$\delta _P^{\mathrm{exp}}=0.207\pm 0.007.$$
(31)
To get a value for $`\alpha _s(M_\tau )`$ one should simply fit $`\delta _P^{\mathrm{exp}}`$ against $`R_\tau /31`$ as given by eq. (IV.1) or by eq. (27) to get a result corresponding FOPT or CIPT (Fixed Order or Contour Improved PT).

Unfortunately, there is no unique way to assign a theoretical uncertainty $`\delta \alpha _s`$ to the obtained value of $`\alpha _s(M_\tau )`$. In the literature one finds several suggestions. Let us consider them in turn.

1. $`\delta \alpha _s`$ is a half of the shift in $`\alpha _s`$ induced by the last fully computed term in the PT (that is by the one of order $`\alpha _s^3`$ at present).
2. $`\delta \alpha _s`$ is equal to the change in $`\alpha _s`$ caused by varying the normalization point $`\mu `$ around $`M_\tau `$, typically within the range of 1.1–2.5 GeV.

The corresponding results read (terms of order $`\alpha _s^4`$ and higher in eq. (IV.1) and in the $`D`$-function have been set to zero)

$`\alpha _s^{\mathrm{FOPT}}(M_\tau )`$ $`=`$ $`0.345\pm (0.025|0.037),`$ (32)
$`\alpha _s^{\mathrm{CIPT}}(M_\tau )`$ $`=`$ $`0.364\pm (0.012|0.021).`$ (33)
Here the first(second) value in brackets correspond to the use of the first(second) suggestion for the error estimation. After evolution from $`M_\tau `$ to $`M_Z`$ this corresponds to

$`\alpha _s^{\mathrm{FOPT}}(M_Z)`$ $`=`$ $`0.1209\pm (0.0024|0.0037),`$ (34)
$`\alpha _s^{\mathrm{CIPT}}(M_Z)`$ $`=`$ $`0.1229\pm (0.0011|0.0020).`$ (35)
3. $`\delta \alpha _s`$ is equal to the change in $`\alpha _s`$ caused by the uncertainty in the predicted (that is not yet completely known) higher order terms in the perturbative series for $`R_\tau `$;
4. $`\delta \alpha _s`$ is a half of the difference in the $`\alpha _s(M_\tau )`$ as obtained within FOPT and CIPT. This difference comes from different handling of higher order terms.

In order to quantify the error estimates according to 3. and 4. we have shown in Table 2 the results for $`\alpha _s(M_\tau )`$ obtained with various choices for $`d_0^{[1]4}`$, $`d_0^{[1]5}`$ and $`\mu `$. The entries with the choices $`\pm 100`$ for the coefficients illustrate the large change in $`\alpha _s`$ which would result from a failure of PMS and FAC once higher order terms are included. For the plausible values of $`d_0^{[1]4}`$ and $`d_0^{[1]5}`$ we observe an significant decrease of the $`\mu `$ dependence after inclusion of additional terms in the $`\alpha _s`$ series.

Our final predictions for $`\alpha _s^{\mathrm{FOPT}}(M_\tau )`$ and $`\alpha _s^{\mathrm{CIPT}}(M_\tau )`$ are given in the first column of Table 3, together with experimental error222 According to eq. (31) the latter includes (small) uncertainties in the values of $`|V_{ud}|`$ and the nonperturbative correction $`\delta _{NP}`$ in addition to the experimental error per se as displayed in eq. (30). and the combined theory uncertainty. The values of theory uncertainties are listed separately in columns 3,4 and 5. The corresponding values at the scale of $`M_Z`$ are

$`\alpha _s^{\mathrm{FOPT}}(M_Z)`$ $`=`$ $`0.1192\pm 0.0007\pm 0.002,`$ (36)
$`\alpha _s^{\mathrm{CIPT}}(M_Z)`$ $`=`$ $`0.1219\pm 0.001\pm 0.0006.`$ (37)
Thus we observe that the total uncertainty based on a combination of not yet calculated higher order terms, $`\mu `$-dependence and scheme dependence is reduced, once $`\alpha _s^4`$ terms are available. However, the difference between FOPT and CIPT results of roughly 0.02 is a remaining, at the moment irreducible uncertainty333 This is in agreement to the analysis of Ref. Groote:1997cn , where it has been concluded that “the resummed values of $`\alpha _s`$ from $`\tau `$ decay lie outside the convergence radii and can therefore not be obtained from a power series expansion. Regular perturbation series do not converge to their resummed counterparts. The experimental value of $`R_\tau `$ appears to be too large for a fixed order perturbation analysis to apply”. See also krajals ..

It is thus of interest to study this difference as a function of $`m_\tau `$. In practice, this could be applied to sum rules for spectral functions as determined in $`e^+e^{}`$ annihilation.

Therefore, let us consider a hypothetical case of $`\tau `$ lepton with the mass equal 3 GeV. Assuming $`\alpha _s(1.77\mathrm{GeV})=0.334`$ and running this value to 3 GeV via standard 4-loop evolution equation one gets $`\alpha _s(3\mathrm{GeV})=0.2558`$ and predicts

$$\delta _P^{\mathrm{exp}}=0.1353$$

which corresponds to eq. (IV.1) with the $`\alpha _s^4`$ and $`\alpha _s^5`$ terms fixed to their central values (see eqs. (1819)). Let us now investigate the results for $`\alpha _s`$ and the theory error that would result from $`\delta _P^{\mathrm{exp}}=0.1353`$ as a starting point. The corresponding analogs of eqs. (3235) and Table (2) are displayed below as eqs. ((38,39)) and Table (4) correspondingly. The difference between FOPT and CIPT decreases significantly, and this remains true even after extrapolating to $`\alpha _s(M_Z)`$.

$`\alpha _s^{\mathrm{FOPT}}(3\text{GeV})`$ $`=`$ $`0.263\pm (0.013|0.014),`$ (38)
$`\alpha _s^{\mathrm{CIPT}}(3\text{GeV}`$ $`=`$ $`0.265\pm (0.005|0.008),`$ (39)
$`\alpha _s^{\mathrm{FOPT}}(M_Z)`$ $`=`$ $`0.1198\pm (0.002|0.003),`$ (40)
$`\alpha _s^{\mathrm{CIPT}}(M_Z)`$ $`=`$ $`0.1203\pm (0.009|0.0016).`$ (41)
The same coefficients $`d_0^{[1]4}`$ and $`d_0^{[1]5}`$ can be used to predict KS ; Kataev:1994wd  (the non-singlet part of) corrections of orders $`\alpha _s^4`$ and $`\alpha _s^5`$ to the $`R`$-ratio in Z decays:

$$R(n_f=5)=1+a_s+1.409a_s^212.77a_s^3+(97\pm 44)a_s^4+(76\pm 230)a_s^5.$$

It is also of interest to display $`R(s)`$ for $`n_f=3`$ and 4, which is accessible at $`e^+e^{}`$ colliders at lower energies.

$$R(n_f=4)=1+a_s+1.525a_s^211.52a_s^3+(112\pm 30)a_s^4+(245\pm 160)a_s^5.$$

$$R(n_f=3)=1+a_s+1.640a_s^210.28a_s^3+(129\pm 16)a_s^4+(635\pm 100)a_s^5.$$

Our results are close to those of KS ; Kataev:1994wd , they employ, however, the additional information from Baikov:2001aa ; vanRitbergen:1997va .

These formulae demonstrate rather good convergency for $`n_f=5`$ and a reasonably good one for $`n_f=3`$ and $`4`$ if our predictions for the coefficients $`d_0^{[1]4}`$, $`d_0^{[1]5}`$ deviate from the true values within the assumed error margins.

### IV.2 Quadratic mass corrections

Let us first discuss the function $`D_2^{[1]}`$. The FAC/PMS predictions can be easily obtained following CheKniSir97 ; they are listed in Table 5.

We again restore the $`n_f`$ dependence of the coefficient $`d_2^{[1]4}`$ as predicted by FAC/PMS:

$$d_2^{[1]4}(\text{FAC/PMS})=1931.44281.956n_f+9.0294n_f^20.0156289n_f^3$$
(42)
as well as that of $`d_2^{[1]3}`$:

$$d_2^{[1]3}(\text{FAC})=123.65411.6638n_f+0.133293n_f^2,$$
(43)
$$d_2^{[1]3}(\text{PMS})=125.97512.0028n_f+0.134769n_f^2$$
(44)
to be compared with

$$d_2^{[1]3}(\text{exact})=148.97814.3097n_f+0.16463n_f^2.$$

The comparison of estimates and exact results reveals a picture qualitatively similar to the massless case but with some modifications. A few important observations are in order.

1. All three terms of the $`n_f`$ expansion of $`d_2^{[1]3}`$ are successfully predicted within about 20 % accuracy.

2. Unlike the massless case the agreement between FAC/PMS predictions for the coefficient $`d_2^{[1]3}`$ for $`n_f=3,4,5`$ and the corresponding exact numbers is within the range of 15%–20%. On the other hand, the estimation of accuracy of the $`\alpha _s^3`$ fixed $`n_f`$ predictions obtained exclusively from the knowledge of the subleading contribution of $`O(\alpha _s^3n_f)`$ is of the right order of magnitude but somewhat less. All this is probably a consequence of a significantly larger $`n_f`$ independent contribution.

3. At $`𝒪(\alpha _s^4)`$ the exact result for the full coefficient

$$d_2^{[1]4}=d_{2,0}^{[1]4}+d_{2,1}^{[1]4}n_f+d_{2,2}^{[1]4}n_f^2+d_{2,3}^{[1]4}n_f^3$$

is unknown, apart from its leading and subleading terms in $`n_f`$. The prediction for the subleading coefficient $`d_{2,2}^{[1]4}=9.0`$ is by 50% larger than the exact value 6.07. The predicted values for $`d_{2,0}^{[1]4}`$ and $`d_{2,1}^{[1]4}`$ are very large. In view of this largeness, the estimate $`(d_{2,2}^{[1]4}|_{exact}d_{2,2}^{[1]4}|_{FAC})n_f^2`$ (= 30, 50, 80 for $`n_f=3,4,5`$ respectively) looks somewhat too optimistic. Therefore we assign a conservative 30% uncertainty to the fixed $`n_f`$ predictions listed in the fifth column of the Table 5.

At last, we repeat the analysis for the function $`\mathrm{\Pi }_2^{[2]}`$. The results of FAC/PMS optimization methods are given in Table 6. Using the values of $`k_3^{\mathrm{PMS}}`$ from the table for $`n_f=3,4,5`$ we reconstruct the corresponding full $`n_f`$ dependence:

$$k_2^{[1]3}=294.47233.2429n_f+0.696598n_f^2.$$

The comparison with the known terms of order $`n_f`$ and $`n_f^2`$ ($`32.0843n_f,0.464663nf^2`$) demonstrates a remarkably good agreement for the subleading $`n_f`$ contribution. The 50% error in the predicted value of the $`n_f^2`$ contribution looks natural as the corresponding coefficient is small. Following the same line of reasoning as above we have assigned a 30% uncertainty to the $`O(\alpha _s^3)`$ fixed $`n_f`$ result.

To get a general idea about the size of the $`\alpha _s^4`$ contribution to eq. (14) we used FAC/PMS and the predicted $`\alpha _s^3`$ coefficient. The results are listed in the forth column of Table 6.

Let us now consider the effect of $`\alpha _s^3`$ and $`\alpha _s^4`$ corrections on the determination of the strange quark mass444For a recent review of various attempts to extract the strange quark mass from $`\tau `$ data see Ref. Maltman:2002wb .. The mass correction to the $`R_\tau `$ depend on both functions $`D_2^{[1]}`$ and $`\mathrm{\Pi }_2^{[2]}`$. Let us use the central ALEPH value of $`\alpha _s(M_\tau )=0.334`$ when estimating the size of the perturbative corrections. For fixed order one finds the mass correction to the total rate:

$`\delta _{us,2}^{00}`$ $`=`$ $`8{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.+5.33a_s+46.0a_s^2+284a_s^3`$ (45)
$`+`$ $`0.75a_s^3k_2^{[2]3}+a_s^4(723.+0.25d_2^{[1]4}+9.84k_2^{[2]3}+0.75k_2^{[2]4}))`$
$`=`$ $`8{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.+0.567+0.520+0.521\pm 0.05+0.593)`$
$`=`$ $`8{\displaystyle \frac{m_s^2}{M_\tau ^2}}(3.2\pm 0.6),`$
where in the last line we have assumed the (maximal!) value of the $`𝒪(\alpha _s^4)`$ term as an estimate of the theoretical uncertainty (this convention will be used also below).

For the “contour improved” series one obtains

$`\stackrel{~}{\delta }_{us,2}^{00}`$ $`=`$ $`8{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.44+3.65a_s+30.9a_s^2+72.2a_s^3+1.18a_s^3k_2^{[2],3}`$ (46)
$`+a_s^4(0.678d_2^{[1]4}+1.06k_2^{[2]4}))`$
$`=`$ $`8{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.44+0.389+0.349+0.371\pm 0.09+0.403)`$
$`=`$ $`8{\displaystyle \frac{m_s^2}{M_\tau ^2}}(2.95\pm 0.4).`$
Now we consider the contributions of spin 1 and spin 0 separately. The lowest moments ($`L=0`$) of the spin-dependent functions depend on a nonperturbative quantity and, thus, can not be treated perturbatively in principle Chetyrkin:1998ej .

For the spin one part and for $`(k,l)=(0,1)`$ we find

$`\delta _{us,2}^{(1)01}`$ $`=`$ $`5{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.+4.83a+35.7a_s^2+276.a_s^3+a_s^4(1350+d_2^{[1]4})))`$ (47)
$`=`$ $`5{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.+0.514+0.404+0.331+0.326)`$
$`=`$ $`5{\displaystyle \frac{m_s^2}{M_\tau ^2}}(2.58\pm 0.33)`$
and

$`\stackrel{~}{\delta }_{us,2}^{(1)01}`$ $`=`$ $`5{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.37+2.55a+16.1a_s^2+135a_s^3+0.895a_s^4d_2^{[1]4})`$ (48)
$`=`$ $`5{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.37+0.271+0.182+0.163+0.137)`$
$`=`$ $`5{\displaystyle \frac{m_s^2}{M_\tau ^2}}(2.12\pm 0.14).`$
Note that spin 1 contribution is determined by the component $`\mathrm{\Pi }^{[1]}`$ alone and is known up to third order. Clearly, this series is decreasing in a reasonable way (comparable to the behaviour of $`\stackrel{~}{\delta }_{us,2}^{00}`$) and, at the same time, only moderately dependent on the improvement prescription with $`\stackrel{~}{\delta }_{us,2}^{(1)01}/\delta _{us,2}^{(1)01}=0.82`$. On the basis of Eq. (48) this moment might well serve for a reliable $`m_s`$ determination, with a sufficiently careful interpretation of the theoretical uncertainty.

The corresponding spin zero part is, per se, proportional to $`m_s^2`$ (not counting non-perturbative, so-called “condensate” contributions) and thus could be considered as ideal for a measurement of $`m_s`$. However, the behaviour of the perturbative series

$`\delta _{us,2}^{(0)01}`$ $`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.+9.33a+110a_s^2+1323a_s^3`$ (49)
$`+a_s^4(12200+d_2^{[1]4}+17.5k_2^{[2]3})`$
$`=`$ $`{\displaystyle \frac{m_s^2}{M_\tau ^2}}(1.+0.992+1.24+1.59+2.16)`$
$`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{m_s^2}{M_\tau ^2}}(7.0\pm 2)`$
and

$`\stackrel{~}{\delta }_{us,2}^{(0)01}`$ $`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{m_s^2}{M_\tau ^2}}(3.19+11.2a+126.a_s^2+289.a_s^3+6.63a_s^3k_2^{[2]3}`$ (50)
$`+a_s^4(2.71d_2^{[1]4}+7.76k_2^{[2]4}))`$
$`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{m_s^2}{M_\tau ^2}}(3.19+1.19+1.42+1.94+2.6\pm 1.31)`$
$`=`$ $`{\displaystyle \frac{3}{2}}{\displaystyle \frac{m_s^2}{M_\tau ^2}}(10.3\pm 2.6)`$
shows a rapid growth of the coefficients. The series is not expected to provide an accurate prediction for the mass effects.

## V Summary

Implications of the newly calculated $`\alpha _s^4n_f^2`$ terms together with the $`\alpha _s^4n_f^3`$ terms for an improved extraction of $`\alpha _s`$ from the $`\tau `$ decay are presented. Arguments are presented in support of predictions for the remaining terms of order $`\alpha _s^4n_f`$ and $`\alpha _s^4n_f^0`$ which are based on FAC or PMS optimization. The complete calculation will lead to a reduction of the theory uncertainty within the frameworks of FOPT or CIPT down to a negligible amount. However, an irreducible difference between the results from these two schemes of $`\delta \alpha _s(M_\tau )0.02`$ corresponding to $`\delta \alpha _s(M_Z)0.002`$ persists even after inclusion of $`𝒪(\alpha _s^4)`$ (and even $`𝒪(\alpha _s^5)`$ terms). Similar investigations, based on data up to higher energies (e. g. for fictitious heavy lepton of 3 GeV of for sum rules based on $`e^+e^{}`$ data) would lead to significantly smaller errors.

New contributions of orders $`𝒪(m_s^2\alpha _s^4n_f^2)`$ and $`𝒪(\alpha _s^4n_f^2)`$ to (axial)vector correlators relevant for the QCD description of the semileptonic $`\tau `$ decay into hadrons are obtained. The moments $`R^{(0,0)}`$ and $`R^{(0,1)}`$ are evaluated separately for spin zero and spin one final states Chetyrkin:1998ej . The results are tested against predictions of FAC/PMS optimization methods. Good agreement is found. This has motivated us to take the full of FAC/PMS predictions as the basis for a new extraction of $`\alpha _s`$ and $`m_s`$ from $`\tau `$-decays with $`𝒪(\alpha _s^4)`$ accuracy. Using $`\delta _P^{\mathrm{exp}}=0.207`$ we find 0.330 and 0.1192 for $`\alpha _s^{\mathrm{FOPT}}(M_\tau )`$ and $`\alpha _s^{\mathrm{FOPT}}(M_Z)`$ respectively. In the framework of contour improved evaluation these values increase to 0.354 and 0.1219 respectively

In contrast to the massless result, the PT series contributing to the $`m_s^2`$ dependent part seem to be barely convergent and the additional higher order terms seemingly do not lead to any significant improvement of the theoretical accuracy in the determination of the strange quark mass from the $`\tau `$ decays.

A slightly more favorable pattern of convergence is observed for the moments of the spin one contribution separately.

## VI Acknowledgments

The authors are grateful to A.A. Pivovarov for a lot of useful discussions. One of us (K.G. Ch.) acknowledges useful correspondence to K. Maltman and A.L. Kataev. This work was supported by the DFG-Forschergruppe “Quantenfeldtheorie, Computeralgebra und Monte-Carlo-Simulation” (contract FOR 264/2-1), by INTAS (grant 00-00313), by RFBR (grant 01-02-16171), by Volkswagen Foundation and by the European Union under contract HPRN-CT-2000-00149.


text:
# Does the Sun Shine by 𝐩𝐩 or CNO Fusion Reactions?

## Abstract

We show that solar neutrino experiments set an upper limit of 7.8% (7.3% including the recent KamLAND measurements) to the fraction of energy that the Sun produces via the CNO fusion cycle, which is an order of magnitude improvement upon the previous limit. New experiments are required to detect CNO neutrinos corresponding to the 1.5% of the solar luminosity that the standard solar model predicts is generated by the CNO cycle.

In 1939, Hans Bethe described in an epochal paper bethe39  two nuclear fusion mechanisms by which main sequence stars like the Sun could produce the energy corresponding to their observed luminosities. The two mechanisms have become known as the $`pp`$ chain and the CNO cycle book . For both the $`pp`$ chain and the CNO cycle the basic energy source is the burning of four protons to form an alpha particle,two positrons, and two neutrinos. Thus

$$4p{}_{}{}^{4}\mathrm{He}+2e^++2\nu _e+25\mathrm{MeV}(\mathrm{to}\mathrm{the}\mathrm{star}).$$
(1)
In the $`pp`$ chain, fusion reactions among elements lighter than $`A=8`$ produce a characteristic set of neutrino fluxes, whose spectral energy shapes are known but whose fluxes must be calculated with a detailed solar model. In the CNO chain, with 12C as a catalyst, 13N and 15O beta decays are the primary source of neutrinos.

The first sentence in Bethe’s paper reads: “It is shown that the most important source of energy in ordinary stars is the reactions of carbon and nitrogen with protons.” Bethe’s conclusion about the dominant role of the CNO cycle relied upon a crude model of the Sun. Over the next two and a half decades, the results of increasingly more accurate laboratory measurements of nuclear fusion reactions and more detailed solar model calculations led to the theoretical inference that the Sun shines primarily by the $`pp`$ chain rather than the CNO cycle. Currently, solar model calculations imply bp00  that 98.5% of the solar luminosity is provided by the $`pp`$ chain and only $`1.5`$% is provided by CNO reactions.

In recent years, there have been many analyses of solar neutrino oscillations, essentially all of which assumed that the CNO neutrino fluxes were equal to their predicted standard solar model values (see Refs. postsno02 ; snodaynight ; bargernc  and references cited therein). However, from the earliest days of solar neutrino research, a primary goal of the field was to test the solar model prediction that the Sun shines by the $`pp`$ chain and not by the CNO cycle bahcall69 . This goal has largely been ignored in the last decade or so as solar neutrino experiments concentrated on the more accessible, higher-energy 8B neutrinos. In this paper, we return to the question of how well we can measure, or set an upper limit to, the CNO neutrino fluxes.

Unfortunately, the standard solar model prediction for the CNO fluxes is difficult to test. Radiochemical experiments with chlorine chlorine  and gallium sage ; gallex ; gno  do not measure the energy of the neutrinos detected; they measure the rate of neutrino induced events above a fixed energy threshold. The neutrino-electron scattering experiments, Kamiokande kamiokande  and Super-Kamiokande superk , provide information about neutrinos but only those that have energies well above the maximum energies of the 13N ($`E_{\mathrm{max}}=1.2`$ MeV) and 15O ($`E_{\mathrm{max}}=1.7`$ MeV) neutrinos. The radiochemical experiments are only sensitive to electron type neutrinos, and the neutrino-electron scattering experiments are primarily sensitive to electron type neutrinos. The heavy water experiment, SNO snodaynight ; snoccnc , measures higher energy neutrinos. The goal of uniquely identifying CNO neutrinos is made even more difficult by the fact that neutrino oscillations can change in an energy dependent way the probability that electron type neutrinos created in the Sun reach the Earth as electron type neutrinos pontecorvo ; msw .

Because of these complications, it was possible to find neutrino oscillation solutions in which 99.95 % of the Sun’s luminosity is supplied by the CNO cycle cno . These ‘large CNO’ oscillation solutions describe well all of the measurements from the chlorine chlorine , SAGE sage , GALLEX gallex , and Kamiokande kamiokande  solar neutrino experiments. Modern solar models do not predict a large CNO contribution to the solar luminosity, but the goal is to test experimentally-not just assume-this prediction.

In this paper, we use data from the chlorine, SAGE, GALLEX, GNO, Super-Kamiokande, and SNO solar neutrino experiments, and from the recent KamLAND kamland  reactor measurements, to set an experimental limit on the CNO contribution to the solar luminosity that is an order of magnitude more stringent than the previous best limit cno . Although individual experiments do not constrain well the CNO fluxes, a global solution to all the available neutrino data provides a powerful upper limit. We also discuss how well future experiments can do in detecting the CNO neutrinos.

Here is our strategy. For each value of the CNO luminosity fraction, $`L_{\mathrm{CNO}}/L_{}`$, we search the two-component neutrino oscillation parameter space with a dense mesh corresponding to the neutrino mass difference $`10^{12}\mathrm{eV}^2<\mathrm{\Delta }m^2<10^3\mathrm{eV}^2`$ (721 mesh points) and mixing angles $`0.0001<\mathrm{tan}^2\theta <10`$ (401 mesh points), as well as the solar neutrino fluxes (see below). \[We verify later that our approximation of only two neutrinos does not limit the validity of the upper bound we derive. See discussion following Eq. (5).\] We calculate the global $`\chi ^2`$ by fitting to all the available data,

$$\chi ^2=\chi _{\mathrm{solar}}^2+\chi _{\mathrm{KamLAND}}^2.$$
(2)
We carry out a global analysis postsno02  of the solar neutrino data letting the neutrino fluxes be free variables and using data from 80 measurements: 44 data points from the Super-Kamiokande zenith-angle energy distribution, 34 data points from the SNO day-night energy spectrum snoccnc ; snodaynight , and 2 radio-chemical rates from Cl chlorine  and Ga sage ; gallex ; gno . We define the $`3\sigma `$ upper limit for the CNO neutrino fluxes by determining when $`\chi ^2=\chi _{min}^2+9`$ after marginalizing over the oscillation parameters $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$ and over the other neutrino fluxes.

Using the data provided in Ref. kamland , we calculate the positron spectrum in the KamLAND detector with the procedures described in Refs. kamland ; prekamlandus ; prekamlandthem . In the absence of neutrino oscillations, we find (in agreement with Ref. kamland ) 86.8 expected neutrino events above 2.6 MeV visible energy for the stated experimental conditions. The positron energy spectrum that we calculate is in excellent agreement with the energy spectrum presented by the KamLAND collaboration. Further details of our analysis of the KamLAND and solar data can be found in Ref. postkamland .

We impose the ‘luminosity constraint’ on the solar neutrino fluxes, i.e., we require that the sum of the thermal energy generation rates associated with each of the solar neutrino fluxes be equal to the solar luminosity spiro . The fraction of the sun’s luminosity that arises from CNO reactions can be written as luminosity  :

$$\frac{L_{CNO}}{L_{}}=\underset{i=N,O,F}{}(\frac{\alpha _i}{10\mathrm{MeV}})a_i\varphi _i,$$
(3)
where the constant $`\alpha _i`$ is the energy provided to the star by nuclear fusion reactions associated with the $`i^{th}`$ neutrino flux, $`a_i`$ is the ratio of the neutrino flux $`\mathrm{\Phi }_i(\mathrm{BP00})`$ of the standard solar model to the characteristic solar photon flux defined by $`L_{}/[4\pi (A.U.)^2(10\mathrm{M}\mathrm{e}\mathrm{V})]`$, and $`\varphi _i`$ is the ratio of the true solar neutrino flux to the neutrino fluxes predicted by the BP00 standard solar model bp00 . Ref. luminosity  presents a detailed derivation of Eq. (3) and the numerical values for the coefficients $`\alpha _i`$ and $`a_i`$.

We treat as free parameters all the solar neutrino fluxes that are normally reported in solar neutrino calculations. There are then 10 free parameters: the two oscillation parameters, $`\mathrm{\Delta }m^2`$ and $`\mathrm{tan}^2\theta `$,and the 8 neutrino fluxes, $`pp`$, $`pep`$, 7Be,8B, and hep (from the $`pp`$ chain) and 13N 15O, and 17F (from the CNO cycle). To speed up the calculations, we made some approximations that we have checked do not affect the accuracy of our search. Two of the solar neutrino fluxes, hep and 17F, are small as a result of nuclear physics considerations. In the initial search calculations, we set hep equal to its solar model value and $`{}_{}{}^{17}\mathrm{F}=\left[\left({}_{}{}^{17}\mathrm{F}/^{13}\mathrm{N}\right)_{\mathrm{solar}\mathrm{model}}\right]^{13}\mathrm{N}`$. We checked that our results are unchanged if the hep solar model flux is multiplied by eight (present experimental bound from the high energy bins at Super-Kamiokande smy ) or if we set the 17F flux equal to the 17N flux. Also, the ratio of the $`pep`$ neutrino flux to the $`pp`$ neutrino flux is fixed to high accuracy because they have the same nuclear matrix element. We have set the ratio equal to the standard solar model value and have checked that our results are unchanged if this ratio is varied by 10% (an enormous change). We set the 13N flux equal to the 15O flux, which is expected in the limit that the CNO contribution to the luminosity is dominant. We also verified that the results of our search are unchanged if we set the ratio of 13N to 15O neutrino fluxes equal to the standard solar model value, the ratio expected if the $`pp`$ contribution is dominant. Finally, we checked several intermediate values of this ratio to see that the upper limit we quote here is robust and valid in all cases.

We find the minimum value of $`\chi ^2`$ for each assumed value of the CNO luminosity fraction by marginalizing over the neutrino oscillation parameters and over the non-CNO neutrino fluxes. We performed the calculations in two stages: first using only the solar neutrino data and second using both the solar neutrino and the KamLAND data. We carried out calculations for oscillations to purely active neutrinos, to purely sterile neutrinos, and to active-sterile admixtures as described in Ref. four  (see also last reference in Ref. prekamlandus ). We considered sterile admixtures that range from the maximum allowed by the recent KamLAND data, $`\mathrm{sin}^2\eta =13`$postkamland , to 25%, 50%, 75%, as well as the extremes of 0% and 100%. For all values of the CNO luminosity fraction, the minimum $`\chi ^2`$ was, as expected, achieved for purely active oscillations.

Figure 1 summarizes our main results. The figure shows $`\mathrm{\Delta }\chi ^2`$ as a function of the CNO luminosity fraction when only solar neutrino data are used (denoted by dotted curves) and when solar and KamLAND data are used (denoted by solid curves). The minimum value of $`\chi ^2`$, relative to which $`\mathrm{\Delta }\chi ^2`$ is measured, is reached in both cases for a zero value of the CNO flux. However, as is apparent from Figure 1, the global $`\chi ^2`$ is essentially flat for all values of $`L_{\mathrm{CNO}}/L_{}<5`$%. Current experiments are not sensitive to CNO neutrino fluxes that correspond to less than 5% of the solar luminosity. For example, the $`\mathrm{\Delta }\chi ^2`$ is only 0.5 between the best fit at $`L_{\mathrm{CNO}}/L_{}=0.0`$ and the standard solar model value of 1.5%.

In all cases, the best fit is a LMA (large mixing angle) MSW solution msw . In principle, new solutions might have been found by allowing $`L_{\mathrm{CNO}}/L_{}`$ to be a free parameter. In practice, no preferred new solutions are found for any value of $`L_{\mathrm{CNO}}/L_{}`$.

For completeness, we explored the entire oscillation parameter space specified just preceding Eq. (2), including the parameter ranges of all previously recognized oscillation solutions postsno02 ; snodaynight ; bargernc . The best-fit results for these searches for other oscillation cases are labeled ‘non-LMA’ in Figure 1. For $`L_{\mathrm{CNO}}/L_{}`$ less than about 5%, the best-fit non-LMA solution is the well-known LOW solution. For the larger values of $`L_{\mathrm{CNO}}/L_{}`$ shown in Figure 1, the best-fit LMA solution is a vacuum solution with $`\mathrm{\Delta }m^2=7.9\times 10^{11}\mathrm{eV}^2`$ and $`\mathrm{tan}^2\theta =0.22`$. We conclude, using both solar neutrino and KamLAND experimental data (using just solar neutrino data), that

$$\frac{L_{\mathrm{CNO}}}{L_{}}<7.3\%(7.8\%)\mathrm{at}3\sigma .$$
(4)
The recent KamLAND measurements reduce the upper limit by 0.5%.

The order of magnitude improvement between the previous limit of 99.95%cno  and the present limit, Eq. (4), is due to the Super-Kamiokande and SNO measurements. The earlier large CNO oscillation solution was confined to small mixing angles, SMA, which cannot fit simultaneously the flat recoil energy spectrum measured by Super-Kamiokande superk  and the total event rates measured by Super-Kamiokande and SNO superk ; snoccnc .

The maximum CNO neutrino flux allowed by the existing experiments is

$$\varphi _{\mathrm{CNO},\mathrm{max}}<3.41\times 10^{10}\mathrm{cm}^2\mathrm{s}^1\left(L_{\mathrm{CNO}}/L_{}\right),$$
(5)
where $`\varphi _{\mathrm{CNO},\mathrm{max}}=\varphi ({}_{}{}^{13}\mathrm{N})_{\mathrm{max}}=\varphi ({}_{}{}^{15}\mathrm{O})_{\mathrm{max}}`$.

We have verified that the upper limits given in Eq. (4) and Eq. (5) are not affected by the approximation of assuming that there is just one mass scale (i. e., two neutrinos). We repeated the analysis assuming the standard three-neutrino mixing scenario invoked to explain both solar and atmospheric data and assumed values for $`\theta _{13}`$ values below the CHOOZ bound chooz , $`\mathrm{tan}^2(\theta _{13})=0.0`$, $`0.03`$ and $`0.06`$ (the CHOOZ bound at 3$`\sigma `$ concha13 ). The minimum $`\chi ^2`$ was, as expected, achieved for $`\mathrm{tan}^2(\theta _{13})=0.0`$.

New solar neutrino experiments are required to measure the CNO contribution to the solar luminosity.

How much can a future 7Be neutrino-electron scattering experiment, BOREXINO borexino  or KamLAND kamland , improve the limit given in Eq. (4)? We find an approximate answer to this question by computing a global $`\chi ^2`$ including the existing solar neutrino data, the KamLAND reactor data, and a simulated BOREXINO rate measurement (simulations guided by Ref. borexino ). We assume that the BOREXINO rate will be consistent with the predicted best fit point from the solar plus KamLAND global fit with a total error of 10% (5%) for the rate measurement. If these assumptions are valid, one will be able to either measure $`L_{\mathrm{CNO}}/L_{}`$ or conclude that $`L_{\mathrm{CNO}}/L_{}<5.6\%(4.9)`$%.

In order to measure the CNO contribution at the 1.5% level predicted by the standard solar model, one must be able to distinguish the continuum 13N and 15O neutrinos from the 7Be and $`pep`$ neutrino lines, as well as from all the sources of background. The appropriate analyses of proposed low energy neutrino-electron scattering detectors have not yet been carried out, so one cannot say for sure whether or not this will be possible. But, it seems very difficult. The energy resolution required to measure the energy of the CNO neutrinos and determine their flux, may, however, be within the reach of low-energy CC experiments lowe .

The solar model predictions for CNO neutrino fluxes are not precise because the CNO fusion reactions are not as well studied as the $`pp`$ reactions adelberger  and because the Coulomb barrier is higher for the CNO reactions, implying a greater sensitivity to details of the solar model. For the standard solar model CNO neutrino fluxes, the $`1\sigma `$ errors vary between 17% and 25% bp00 . A measurement of the CNO neutrino fluxes would constitute a stringent test of the theory of stellar evolution and provide unique information about the solar interior.

###### Acknowledgements.

JNB and CPG acknowledge support from NSF grant No. PHY0070928. MCG-G is supported by the European Union Grant No HPMF-CT-2000-00516, by Spanish Grants No PB98-0693 and CTIDIB/2002/24.


text:
# Universal properties of shortest paths in isotropically correlated random potentials

## 1 Introduction

Optimal paths have been a subject of intensive studies during the recent years. Besides being one of the simplest problems involving disorder, this interest can be traced back to the relevance of this problem to various fields, such as polymer models ref10 ; ref16 ; ref26 , surface growth ref6 , random bond ferromagnets ref2 ; ref3 ; ref4 , spin glasses ref1 , and the traveling salesman problem ref5 .

The model under consideration is easily sketched: given an arbitrary weighted graph, each edge has a particular cost. The optimal or shortest path connecting two sites is the one of minimal weight, which is the sum of all costs along that path. We do not restrict to a particular geometry yet as well as we do not specify the costs more precisely so far. In the simplest case, we choose them from a set of random numbers that are uniformly distributed. In this context, the directed polymer model has drawn the most significant attention during the past years ref10 ; ref16 ; ref7 ; ref8 ; ref9 , where one assumes in $`d=2`$ a simple square lattice being cut along its diagonal and oriented as a triangle with the diagonal as its base. One allows only paths in the direction to the base, i.e., the path cannot turn backwards. The costs of the edges that belong to the shortest path are interpreted as potential energies for a polymer configuration that passes through these edges (or bonds).

We may now ask whether the properties (scaling exponents) of the shortest path are either influenced by the distribution of the random numbers or the geometry of the lattice. The former is still discussed ref10 ; ref11 ; ref12 ; ref13 ; ref14 . As far as the latter is concerned, it seems to be clear that the universal properties are not changed if the randomness is uncorrelated. For this case Schwartz *et al.* ref15  investigated directed and undirected paths in $`d=2,3`$ using Dijkstra’s algorithm to find the shortest path and Marsili and Zhang ref27  used a transfer matrix method approach considering directed and undirected paths up to $`d=6`$. Both state that overhangs exist but nevertheless they suggest that both problems belong to the same universality class, even in high dimensions where overhanging configurations are entropically favored. It is not a priori clear that this observation remains true for correlated disorder. In fact, as we will point out below, the average number of overhangs increases for strongly correlated disorder indicating that they might become relevant for strong enough correlations.

In the present study, we study the universal properties of shortest paths and focus on the effect of isotropically correlated random weights on the scaling exponents. To this end, we consider directed and undirected lattice-graphs in two and three dimensions with bond weights $`\eta _𝐣`$, where the $`d`$-dimensional index vector $`𝐣=i_1,\mathrm{},i_d^d`$ denotes the position of a particular bond in the lattice. The total energy or cost of a path $`𝒫`$ from one end of the lattice (e.g. from one special site or node in the top layer) to the opposite end (e.g. to an arbitrary site or node in the bottom layer) is simply the sum of these bond weights

$$E=\underset{𝐣𝒫}{}\eta _𝐣.$$
(1)
The weights $`\eta _𝐣`$ are correlated positive random variables, which we define below. We choose the index vector $`𝐣`$ of the bonds in such a way that it is identical with the position of the center of the bonds in an Euclidean lattice. E.g. for the square lattice ($`2d`$) in which we define the upper left corner as the origin, all indices take on odd values. We refer to isotropic correlations if the connected correlation function of the costs $`\eta _𝐣`$ decays with a power law with Euclidean distance of two bonds, viz.

$$G(𝐣,𝐣^{})=\eta _𝐣\eta _𝐣^{}\eta _𝐣\eta _𝐣^{}|𝐣𝐣^{}|^{2\rho 1},$$
(2)
where $`\rho <1/2`$. Here $`\mathrm{}`$ denotes the disorder average, i.e. an average of the probability distribution of $`\eta _𝐣`$.

For $`\rho =0`$ model (1) with (2) represents an effective single vortex line model for interacting vortex lines in a random vector potential natter . We would like to emphasize that isotropic correlations have not been investigated so far: For historical reasons one discriminates between $`d1`$ transverse or spatial directions and 1 longitudinal or time direction. This can be traced back to the relation between directed polymers and the KPZ-equation. If the randomness is only correlated in $`d1`$ dimensions and uncorrelated in the remaining direction one refers to spatially correlated randomness, if the correlations are only present in time direction one refers to temporal correlations. Let $`𝐣=(𝐱,t)=((i_1,\mathrm{},i_{d1}),t)`$. Then we can describe both cases by $`G|𝐱𝐱^{}|^{2\rho 1}\delta (tt^{})`$ (spatial correlations) and $`G|tt^{}|^{2\theta 1}\delta (𝐱𝐱^{})`$ (temporal correlations).

The first systematic numerical work done on spatially correlated noise is due to Amar, Lam and Family ref22 . Their results are in agreement with the dynamical ref23  and functional RG ref24 ; ref25  predictions. Subsequent careful work by Peng *et al.* ref21  and Pang *et al.* ref17  was done. As far as temporal correlations are concerned, we refer to the simulation of a ballistic deposition model by Lam, Sander and Wolf ref20 .

In this paper we study isotropically correlated disorder, it is organized as follows: In section II we introduce the models, the numerical method and the quantities that we are interested in. In section III we present our results and in section IV we summarize our findings.

One remark on the notation: We will use the words graph and lattice as well as costs and energy, node and site and edge and bond synonymously throughout the paper.

## 2 Model

The *undirected graph* can be described as follows (see Fig.1b)): we choose a simple lattice structure and define one longitudinal and $`d1`$ transversal directions. We will refer to them by means of $`t`$ and $`𝐱`$ respectively. We assume the lattice to be periodic in space ($`𝐱`$) with period $`H`$ and $`L`$ to be the longitudinal size. We choose the origin of the coordinate system as being the starting point (the source) of the path. We assign a particular amount of energy to each bond, whereby these energies are isotropically correlated. Generating these random numbers we follow the method introduced by Pang *et al.* ref17 . We infer periodicity and symmetry of the correlator $`g_\rho `$ in any direction, where

$$g_\rho (\mathrm{\Delta }𝐣):=G(𝐣,𝐣+\mathrm{\Delta }𝐣)=|\mathrm{\Delta }𝐣|^{2\rho 1}$$
(3)
with $`g_\rho (\mathrm{},(\mathrm{\Delta }𝐣)_i+p_i,\mathrm{}):=g_\rho (\mathrm{},(\mathrm{\Delta }𝐣)_i,\mathrm{})`$ and
$`g_\rho (\mathrm{},(\mathrm{\Delta }𝐣)_i,\mathrm{})=g_\rho (\mathrm{},p_i/2|(\mathrm{\Delta }𝐣)_ip_i/2|,\mathrm{})`$ in case of a period $`p_i`$ in direction $`i`$. In contrast to pure spatial and temporal correlations where the correlator is taken to be the product of two separate ones, one for the time direction and one for the remaining spatial coordinates, here it is due to a generic vector $`\mathrm{\Delta }𝐣`$. Adapting the correlator to the lattice, we require a period $`H`$ in the transversal and $`2L`$ in the longitudinal direction. By means of the factor $`2`$ we guarantee that $`g_\rho (t)`$ has the required form (3) in the range $`1tL`$. Its Fourier transform yields $`S_\rho (𝐤)`$ such that the correlator in the $`𝐤`$-space is given by $`\eta _𝐤\eta _𝐤^{}\delta _{𝐤+𝐤^{},0}S_\rho (𝐤)`$. Choosing $`\eta _𝐤\sqrt{S_\rho (𝐤)}(r_𝐤1/2)\mathrm{exp}(2\pi i\varphi _𝐤)`$, that relation can be fulfilled, where $`r_𝐤`$ and $`\varphi _𝐤`$ are random variables uniformly distributed between $`0`$ and $`1`$. A transformation back to real space provides the random numbers correlated according to the power law rule (see also Appendix A).

From the set of $`L^{d1}`$ optimal paths that connect the source $`(𝐱,t)=(\mathrm{𝟎},0)`$ with nodes of the bottom layer with coordinates $`(𝐱,L)`$, we select the shortest one. Technically this is achieved not by repeating the same calculation $`L^{d1}`$ times (i.e. once for each end-point) but by introducing an extra target node connected to the bottom layer by zero-weight bonds. The algorithm that computes the optimal paths in polynomial time is Dijkstra’s algorithm that works in any graph with non-negative weights.

Our study is focused on the universal characteristics of the optimal path, i.e. on the scaling exponents $`\nu `$, the roughness exponent, and $`\omega `$, the energy fluctuation exponent. $`\nu `$ describes the fluctuations of the path with regard to a line parallel to the $`t`$-axes that is shifted to the origin by an amount matching to the mean position $`𝐱`$ of the polymer. We refer to these fluctuations by $`D`$. Due to the direct mapping between the model of shortest paths and growing interfaces ref6 , it is immediately seen that the energy of the polymer also fluctuates, where we consider several realizations of disorder: we have

$$Dt^\nu \mathrm{\Delta }Et^\omega $$
(4)
There is no need to modify these relations when we consider *directed paths*. For this purpose we introduce directed bonds, e.g., we restrict to the positive axes. In order to determine the roughness of the polymer, we refer to a line parallel to the bisecting line $`(𝐱,t)(1^{d1},1)`$ that crosses the origin. For each size $`L`$ we can determine the distance of the target node to that line by considering the projection of its position vector onto the bisector. That is, a directed path connects the source $`(𝐱,t)=(\mathrm{𝟎},0)`$ and the sites on the line between $`(0,L)`$ and $`(L,0)`$, what has to be extended to the notion of a plane in $`d=3`$.

As far as undirected paths are concerned one has to make clear how $`𝐱(t)`$ can be defined if overhangs appear. In that case $`𝐱(t)`$ is not a single valued quantity anymore. We have checked the relevance of several choices but finally we took $`𝐱(t)=\text{max}\{𝐱_i(t)\}`$ with $`t`$ constant for all $`𝐱_i`$. The numerical results that we present are independent of this choice.

As there are $`𝒪(2^n)`$ undirected paths across the lattice, where $`n`$ is the number of nodes, an efficient method is needed that terminates within a reasonable time, such as Dijkstra’s algorithm ref15 ; ref18 , that is able to generate the shortest path in polynomial time. In this algorithm the path is successively constructed, i.e., one obtains not only a single path but a cluster of them with energy labels smaller than a certain limit, so that a growth front is established. This cluster contains the so called permanently labeled sites. The algorithm proceeds by extending the front by this site that is the nearest neighbor to it with respect to the energy. If the growth front reaches the base of the lattice we are enabled to reconstruct the shortest path. The growth process of that front can directly be mapped onto a growing Eden cluster ref6  providing the same scaling behavior as directed polymers.

If we assume uncorrelated costs the roughness of the growth front scales like the energy fluctuations of the polymer. Since the properties of the surface are described by the KPZ-equation and, consequently, the height-height correlation increases according to a power law with exponent $`\omega =1/3`$, we expect in $`d=2`$ the scaling exponents of the optimal path of size:

$$\omega _{\text{OP}}=\beta _{\text{KPZ}}=1/3\nu _{\text{OP}}=1/z_{\text{KPZ}}=2/3$$

Exact values are only accessible in $`d=2`$. The value of the exponent $`\nu `$ can be extracted from $`\omega `$ if we take into account the exponent relation

$$2\nu \omega =1$$

which holds for uncorrelated noise. This relation follows from the Galilean invariance of the KPZ-equation and the transfer to shortest paths afterwards. In $`d=3`$ there exist no analytical predictions. Nevertheless, the estimates of numerous numerical studies ref16 ; ref10 ; ref15  yield

$$\omega 0.19\nu 0.62$$

As long as the Galilean invariance holds, the scaling relation remains unchanged. This invariance is not altered in the presence of spatial (transversal) correlations but it is broken in the presence of temporal (longitudinal) correlations.

A Flory type argument ref8  leads to the following estimate $`\nu _F`$ of the roughness exponent as a function of $`\rho `$. A continuum Hamiltonian of the energy (1) has to include an elastic part, since this is generated in a coarse graining procedure.

$$=𝑑t\left[\frac{\lambda }{2}(_t𝐱)^2+\eta (𝐱,t)\right].$$
(5)
Rescaling $`t`$ with a factor $`b`$ and $`x`$ with a factor $`b^\nu `$ the elastic term scales with $`b^{2\nu 2}`$ and the disorder term with $`b^{\rho 1/2}`$ (because of the power law decay of the disorder correlations), resulting in

$$\nu _F=\frac{1}{2}\rho +\frac{3}{4}.$$
(6)
In case of uncorrelated random numbers the disorder term scales with $`b^{\nu (d1)/21/2}`$, resulting in

$$\nu _{uncorr.}=\frac{3}{d+3}.$$
(7)
Hence we expect

$$\nu =\{\begin{array}{ccc}\nu _{uncorr.}\hfill & \mathrm{for}& \rho \frac{3(1d)}{2(3+d)}\hfill \\ \nu _F\hfill & \mathrm{for}& \rho >\frac{3(1d)}{2(3+d)}\hfill \end{array}$$
(8)
This simplified scaling picture should yield at least a lower bound for the roughness exponent.

## 3 Numerical results

In addition to the relations in (4), we define the two exponents $`\gamma `$ and $`\delta `$ by $`lL^\gamma `$ and $`BL^\delta `$, where $`B`$ is the number of backward bonds with respect to time, and $`l`$ is the total length of the path. We also determine the fractal dimension $`d_c`$ of the shortest path cluster $`ML^{d_c}`$, where $`M`$ is the mass of its surface. The shortest path cluster consists of all nodes with labels smaller than a maximal one given by the shortest path weight from the source to the base. As far as Dijkstra’s algorithm is concerned, its surface is constituted by all the sites that are part of the growth front with at least one nearest neighbor that is not yet permanently labeled.

The scaling exponents are extracted from a set of data that reproduces the simulation of several lattice sizes (Fig.2).

The statistical error is usually smaller than the symbol size. We adapt the transversal expansion $`H`$ to the size $`L`$ in such a way that we eliminate further effects on our data by increasing $`H`$, even if $`\rho =\rho _{max}`$. Finally, we consider lattices of size $`H4L`$ if $`L2048`$ and $`H=4096`$ if $`L=2048`$ in $`d=2`$, as well as $`H=128`$ if $`L32`$ and $`H=256`$ if $`L64`$ in $`d=3`$. At any time, we forbid the shortest path cluster spanning the lattice in order to avoid saturation effects. This requirement cannot be satisfied for all $`L=128`$ in case of $`d=3`$ and $`\rho `$ close to $`\rho _{max}`$.

Our results are averages over more than $`10000`$ disorder realizations per size $`L`$. The generation of correlated random numbers is the most time consuming part. For each sample we have to perform a Fourier transformation of $`N=2^dLH^{d1}`$ numbers twice. On the average approximately $`90\text{\%}`$ of the CPU time is needed for this execution. Some more information about the generation is given in Appendix A.

As a first check we studied the directed polymer problem in $`d=2`$ with only spatially correlated bond weights as it was done by Pang *et al.* ref17 . It can be seen from Fig.3 that we are in very good agreement with their results. In case of strongest correlations ($`\rho =0.40`$) we obtain $`\nu =0.71\pm 0.01`$ and $`\omega =0.43\pm 0.01`$. It is quite evident that correlations are only relevant in the regime $`\rho >0`$.

### Two dimension ($`d=2`$)

In $`d=2`$ we find that the roughness exponent $`\nu `$ does not depend on the directedness of the lattice. Overhangs do not play an important role and the undirected path can be regarded as a directed one (Fig.4). Independent from the strength of the correlations all exponents $`\nu `$ are smaller than the critical value 1. As far as undirected shortest paths are concerned, exceeding that critical value leads to fractal objects that cannot become directed, even on large scales. The errorbars depicted in Fig.4 are not the result of the least square fit but estimates of the minimal and maximal slopes being in nearly perfect agreement with our data.

In contrast, we obtain a less significant data collapse with respect to the energy fluctuation exponent $`\omega `$, if $`\rho >0`$. This may be affected by a statistics that has room of improvement but it indicates a tendency that is especially noticeable in $`d=3`$: the stronger the correlation the more significant variations in $`\omega `$ occur. The reason for this relation is not clear to us. In both cases the exponent relation can be satisfied (Fig.5) where we are in better agreement with $`2\nu \omega =1`$ for undirected paths. We learn from Fig.3 that isotropic correlations in the randomness are relevant for $`\rho >0`$. The scaling of the energy fluctuations is much more sensitive to passing from white noise to weak correlations ($`\rho =0.5`$) in comparison to the roughness. Whereas $`\nu `$ keeps constant ($`\nu =0.66\pm 0.01`$), $`\omega `$ reduces from $`\omega =0.32\pm 0.01`$ to $`\omega =0.29\pm 0.01`$. The number of backward bonds in the NDOP model remains negligibly small. It is $`B15`$ for $`L=1024`$ if $`\rho =0.40`$ where $`\delta 1.0`$. In that case almost every path has at least one such bond. According to these results we obtain $`lL`$.

Although the roughness exponent increases significantly for $`\rho 0.3`$ and might even come closer to 1 for $`\rho 0.6`$ it stays still smaller than one. Also the length of the shortest path scales linearly with $`L`$ for the values of $`\rho `$ that we could study. Both observations imply a non-fractal shortest path. However, the properties of the shortest path cluster change remarkably when coming closer and closer to $`\rho =0.40`$, where the fractal dimension of its surface becomes $`d_c1.1`$ instead of $`d_c1`$ if $`\rho =0`$. In Fig.6 one can see that for weak correlations the surface of the shortest path cluster is a semicircle whereas for larger correlations it becomes topologically more complicated.

Not only the disorder averaged roughness $`D`$ scales with $`L^\nu `$ but the whole probability distribution $`P_L(D)`$: It is $`P_L(D)=L^\nu p(D/L^\nu )`$ as we show in Fig.7 for $`\rho =0.40`$. For the scaling function $`p(x)`$ we can fit a log-normal distribution given by $`p(x)=(2\pi \sigma ^2)^{1/2}\mathrm{exp}((\mathrm{ln}D/D_0)^2/2\sigma ^2)`$. The parameters $`\sigma `$ and $`D_0`$ do not depend on $`\rho `$ for $`\rho 0`$ (where $`D_00.13`$, $`\sigma 0.58`$) and vary slightly with $`\rho `$ for $`\rho >0`$.

### Three dimensions ($`d=3`$)

The results in $`d=3`$ are qualitatively similar to those of the preceeding section. For uncorrelated randomness we obtain $`\nu =0.62\pm 0.02`$ and $`\omega =0.22\pm 0.01`$ in agreement with ref10 ; ref16 ; ref15 . Both $`\nu `$ and $`\omega `$ increase monotonuously with $`\rho `$, i.e. increasing correlations. While $`\nu `$ does so without any difference between directed and undirected paths, $`\omega `$ differs from this behavior: the stronger the correlations the more estimates for the exponents deviate from each other (Fig.8).

As $`\nu <1`$ the undirected path becomes directed on large scales. We should emphasize that obtaining data in the regime $`\rho >0.25`$ for $`d=3`$ is much more delicate than in $`d=2`$. The local slopes of the energy fluctuations indicate that even for the maximal size $`L=128`$ we are not yet in the asymptotic regime for $`\rho >0.3`$. Therefore we explicitly restrict ourselves to values $`\rho 0.3`$. This result refers to both kinds of lattices and, so is not a consequence of saturation effects. The distinct behavior of both kinds of paths can be demonstrated more obviously by plotting the exponent relation $`2\nu \omega =1`$ (see inset of Fig.8).

Even though the scaling regimes in $`d=3`$ are not very wide our estimates for $`\nu `$ and $`\omega `$ deviate significantly from their values in the uncorrelated case. Significantly means here that their estimated errors (obtained in the same way as in our $`2d`$-study before) are smaller than the deviation from the uncorrelated values. Therefore we can infer that for strong enough correlations ($`\rho 0.3`$) the DOP and NDOP problems constitute new universality classes. The precise determination of the critical value for $`\rho `$, beyond which correlations become relevant, is, however, beyond our numerical precision.

For completeness we mention that even in case of $`L=128`$ and $`\rho =0.40`$ the number of bonds turning backward is negligibly small: $`B1.6`$ compared to $`B3`$ in $`d=2`$. We do not expect a significant change in $`B`$ for larger lattices as this quantity scales by an exponent $`\delta 1.0`$ but the number of paths including such bonds should increase from $`60\%`$ to $`100\%`$.

In contrast to the difficulties above, the fractal dimension $`d_c`$ of the surface can be extracted quite clearly. It is $`d_c2.35\pm 0.10`$ if we choose $`\rho =0.40`$ in contrast to $`d_c2.0`$ for uncorrelated numbers. We illustrate the change of the cluster by cutting the system along the $`xt`$plane. Fig.9 shows such cuts for $`y=H/2`$ where $`a)`$ denotes random costs and $`b)`$ corresponds to strongest correlations.

## 4 Summary

In the present paper we have considered the effect of isotropically correlated bond weights on the scaling behavior of both undirected and directed paths in two and three dimensions. We found that in $`d=2`$ the algebraic correlation of the disorder are relevant for a decay slower than $`1/r`$ (i.e. $`\rho =0`$) and the roughness exponent $`\nu `$ and the energy fluctuation exponent $`\omega `$ increase monotonously for $`\rho 0`$. The Flory estimate (8) for the roughness exponent is fulfilled as a lower bound of our numerical estimates. A precise value for $`\rho _c`$, above which correlations modify the universality class, is hard to estimate numerically, but we observe that it is close to and slightly smaller than zero in $`d=2`$ and $`d=3`$, in agreement with the Flory argument (8).

Moreover, the results in $`d=2`$ indicate that the scaling exponents are independent from the directedness even in case of very strong correlations. In contrast to this, in $`d=3`$ directed and undirected lattices yield different results for $`\rho >\rho _c0`$ indicating that both cases constitute different universality classes. In $`2d`$ the scaling relation $`2\nu \omega `$ appears to remain valid even for stronge correlations (although Galilean invariance is broken), whereas for $`d=3`$ we observe in the directed case significant deviations from it (in the undirected case it appears to remain valid).

Finally we could exclude the possibility that optimal paths become fractal for strong disorder correlations as long as they decay algebraically. As a consequence $`\nu `$ stays smaller than one and overhangs turn out to be irrelevant. This behavior changes if we consider the following disorder correlations: $`(\eta _𝐣\eta _{𝐣+𝐫})^2r^\alpha `$, with $`\alpha >0`$, which increase with distance $`r=|𝐫|`$. These type of correlators are relevant for an effective model of dislocations in a vortex line lattice vinokur . In this case we found that the optimal paths are fractals with a fractal dimension significantly larger than one and depending on the value of $`\alpha `$ ref28 .

###### Acknowledgements.

This work has been supported by the Deutsche Forschungsgemeinschaft (DFG).

## Appendix A Appendix

Here we explain how to initialize the correlation function $`g_\rho (𝐣)`$ and we show that $`N=2^dLH^{d1}`$ correlated random numbers have to be created for our purposes. We need to generate $`2N`$ uncorrelated random numbers and have to perform a Fourier transformation twice. If we consider the simple square lattice of nodes with lattice spacing $`a`$ we can immediately see that we have to distinguish between positions of nodes and positions of bonds. The latter build a simple lattice that is rotated towards the node lattice by $`\pi /4`$ with a lattice spacing enlarged by a factor $`\sqrt{2}`$ (Fig.10). But we have to create correlations according to positions of *bonds*. Let us assume now this lattice to be a two dimensional array. The most comfortable solution of the generation of correlations due to bond positions is to keep the original lattice structure, meaning we switch to a spacing $`a/2`$ and finally focus on those positions within this array that correspond to bonds in our lattice. Then we define the center of that lattice (array) and initialize each position of the array by $`|𝐣|^{2\rho 1}`$, where $`H/2𝐣_i<H/2`$ or $`L/2𝐣_i<L/2`$ is the position vector. This array is called the correlator $`g_\rho (𝐣)`$. There are two problems: by doing so, we generate twice the quantum of correlated random numbers we really need (the information of the black squared positions) and we still have to discuss how to define $`g_\rho (0)`$. As mentioned in section II, we already have to take into account a factor $`2`$ from the longitudinal expansion so that $`N`$ becomes $`N=22^dLH^{d1}`$, where again all sites in the left part of Fig.10 would be occupied. In order to restrict memory usage we compress the lattice along one axes (here the x-axes), whereby we loose the information of all the positions indicated by open symbols. By doing so, we avoid the necessity of defining $`g_\rho (0)`$ corresponding to a virtual node at the origin and, consequently, being not part of the right side of Fig.10. More important, $`N`$ is reduced by a factor 2. Note that the Fourier transform on the compressed lattice also yields the desired correlations.