text:
# Inattainability of Carnot efficiency in the Brownian heat engine

## Abstract

We discuss the reversibility of Brownian heat engine. We perform asymptotic analysis of Kramers equation on Büttiker-Landauer system and show quantitatively that Carnot efficiency is inattainable even in a fully overdamping limit. The inattainability is attributed to the inevitable irreversible heat flow over the temperature boundary.

How efficiently can Brownian heat engine work? This question is important not only for the construction of theory of molecular motors but also for foundation of non-equilibrium statistical physics. Like Carnot cycle, Brownian heat engine can extract work from the difference of the temperature of heat baths, where Brownian working material operates as a transducer of thermal energy into mechanical work. The feature of this engine is: 1) It operates autonomously. 2) It is driven by finite difference of the temperature of heat baths both of which contact with the working material simultaneously. Thus, this engine works because the system is out of equilibrium. Feynman devised what is called Feynman’s ratchet that can rectify the thermal fluctuation for work using the difference of the temperature of two thermal baths. Büttiker and Landauer proposed a simpler type of Brownian motor and pointed out that one could extract work even by the simple heat engine where a Brownian particle is subject to a spatially periodic heat baths in a periodic potential.

One crucial point on the Brownian engines is the efficiency. Feynman claimed that his thermal ratchet can operate reversibly, resulting in Carnot efficiency. Recently, however, some authors claimed that Feynman’s claim was incorrect, while some author supported it: Parrondo and Español discussed that Feynman’s ratchet should not work reversibly since the engine is simultaneously in contact with heat baths at different temperatures. Sekimoto devised so-called ”stochastic energetics” and applied it into Feynman’s ratchet. He showed numerically that the efficiency is much less than that of Carnot. Hondou and Takagi showed that reversible operation of Feynman’s ratchet is impossible using reductio ad absurdum. Magnasco and Stolovitzky studied how the engine generates motion with detailed analysis of its phase space. On the other hand, Sakaguchi claimed that the Feynman’s ratchet could operate reversibly by proposing a ”stochastic boundary condition”. Similar result is also found in ref. (not on Feynman’s ratchet but on Büttiker-Landauer system), on which detailed discussion will be made later. These studies have reminded us that there is difficulty as to the energetic description on Brownian systems, because naive application of conventional energetics formulated in thermodynamic and/or equilibrium system into Brownian system may lead incorrect result.

Operation of Brownian engines is done by the engines themselves and the engines are, therefore, out of  equilibrium. To clarify the non-equilibrium nature of Brownian heat engine and to find how we should apply energetics on it, it is important to make a quantitative analysis of the efficiency without adopting over-simplification, for analysis, that loses the function of a heat engine. Because Feynman’s ratchet is somewhat complex to make a rigorous analysis, it seems suitable to discuss Büttiker-Landauer system that is the simplest system of Brownian motors. Recently, Matsuo and Sasa analyzed the energetics of Büttiker-Landauer system by renormalization method. They claimed the system under quasi-static process approaches Carnot efficiency in overdamping limit. Their analysis was based on a rigorous calculation starting from Kramers equation, and the result is clear except one point: They assumed that the momentum degree of freedom is always in equilibrium with heat bath because the system is overdamping. This assumption is not easy for us to accept because the system is singular at the transition points where the temperature of the heat bath changes suddenly. We conjectured that the essence of the mechanism of Brownian heat engine be concentrated on this singular point and that the nature of these non-equilibrium engines would emerge by the analysis. Thus we will discuss the energetics of Bütikker-Landauer system with paying our attention to the transition points. The result will also give us an insight about how we should apply “stochastic energetics” to the overdamping systems with space-dependent temperature.

Let us consider a one-dimensional Brownian system that Bütikker and Landauer discussed, where working particles operate due to the broken uniformity of the temperature of the heat baths. . While Bütikker and Landauer started their discussion from overdamping equation of the system, we start from more basic standpoint of underdamped description, from which the overdamped equation is obtained by eliminating the momentum variable. The probability density in phase space, $`\rho (p,q)`$, obeys Kramers equation:

$`{\displaystyle \frac{\rho (p,q)}{t}}=\left({\displaystyle \frac{J_q}{q}}+{\displaystyle \frac{J_p}{p}}\right)=K(q){\displaystyle \frac{\rho (p,q)}{p}}{\displaystyle \frac{p}{m}}{\displaystyle \frac{\rho (p,q)}{q}}`$

$$+\frac{\gamma }{m}\frac{}{p}\left[p\rho (p,q)+mk_BT(q)\frac{\rho (p,q)}{p}\right],$$
(1)
where $`K(q)=\frac{U}{q}`$; $`\gamma `$, $`m`$, $`J_q`$, and $`J_p`$ are a friction constant, mass of a particle, and probability current in space and that in momentum, respectively. The potential $`U(q)`$ satisfies, $`U(q)=U_r(q)+gq`$, where $`U_r(q+L)=U_r(q)`$, $`g`$ $`(>0)`$ is a gradient of the global slope (load) and $`L`$ is a period. The temperature has the same spatial period as the potential, $`T(q+L)=T(q)`$. In this Büttiker-Landauer system, there are two heat baths of which the temperature is $`T_h`$ (for hot bath) and $`T_c`$ (for cold bath), respectively. Thus, there are the two transition points in a spatial period where the thermal bath affecting the particle changes. Here we restrict ourselves to the case that $`\frac{T_c}{T_h}=O(1)`$ for simplicity. The system is known to operate as a molecular engine because the particle can move against global gradient of the potential. The globally unidirectional motion was attributed to the difference of the temperatures of the baths, since the hot bath can activate working particle more than the cold bath. Suppose that two working particle climb the potential, where one is in a hot bath and the other is in a cold bath. Then the working particle in a hot bath reaches the top of the potential hill more frequently than that in a cold bath, leading to the global motion in the present system. Thus, one can store work in proportion to the probability current. To make energetic analysis, we consider a “replica” particle, of which the energy is $`E=p^2/2m+U(q)`$. Here, the ensemble average over the replicas corresponds to thermodynamic limit.

It was shown that the efficiency of this engine can have Carnot efficiency if the irreversible heat transfer at the transition points is physically negligible. Any Brownian motor is irreversible when they operate with finite probability current. Thus, the operation in Carnot efficiency, if possible, must be in the “stalled state” , where the probability current in space disappears, $`J_q(q)=0`$ (in an overdamped description) or $`𝑑pJ_q(p,q)=0`$ (in an underdamped description). “Quasi-static” operation requires this “stalled state”. Therefore, we evaluate whether and how the irreversible heat flows at the transition point by solving stationary solution of Kramers equation at the stalled state(Eq.1). For this purpose, we will restrict ourselves to the special region, $`q[l_h,l_c]`$, around the transition point, $`q=0`$, where $`l_h`$ ($`l_c`$) are the width between the transition point and a proper point on hot (cold) bath that satisfies the following inequality:

$$l_{th}l_xL_x(x=h\text{or}c),$$
(2)
where $`L_h`$ ($`L_c`$) is a width of hot (cold) bath ($`L_h+L_c=L`$), and $`l_{th}`$ is a characteristic length scale of the transition region in which the probability density is different from that of thermal equilibrium. The length, $`l_{th}`$, is the product of thermal velocity, $`v_{th}(\sqrt{k_BT/m})`$, and the velocity relaxation time, $`\tau (=m/\gamma )`$: $`l_{th}v_{th}\tau `$. The choice of $`l_x`$ does not alter the following result as long as the inequality, Eq.(2), is satisfied. Although we will discuss only one transition region, the asymptotic behavior does not differ in the other transition region. Hereafter we apply the normalization of the probability density $`\rho (p,q)`$ as

$$_{\mathrm{}}^{\mathrm{}}𝑑p_{l_h}^{l_c}𝑑q\rho (p,q)=1.$$
(3)
Now, we will formulate the irreversible heat transfer from a heat bath to the working particle. R.h.s. of Kramers equation (Eq.1) has two parts. The first and the second terms are a Liouville operator on the probability density $`\rho (p,q)`$ and thus preserve the energy. The last term is what describes the energy transfer between the heat bath and the particle, of which the probability current in momentum space is written: $`J_p^{irr}=\frac{\gamma }{m}\left[p\rho (p,q)+mk_BT(q)\frac{\rho (p,q)}{p}\right]`$. Because the probability current disappears, $`J_p^{irr}=0`$, for the probability density of equilibrium $`\rho (p,q)\mathrm{exp}\{\frac{p^2/2m+U(q)}{k_BT}\}`$, the energy flow through $`J_p^{irr}`$ can be sufficiently described only where $`q[l_h,l_c]`$. Average heat transfer from hot bath to the particle per unit time, $`\frac{dQ_h}{dt}`$, reads

$`{\displaystyle \frac{dQ_h}{dt}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑p{\displaystyle _{l_h}^0}𝑑q{\displaystyle \frac{E}{p}}J_p^{irr}`$

$$=_{\mathrm{}}^{\mathrm{}}𝑑p_{l_h}^0𝑑q\frac{p}{m}\frac{\gamma }{m}\left[p\rho (p,q)+mk_BT(q)\frac{\rho (p,q)}{p}\right].$$
(4)
By integration by part through momentum space, $`p`$, and the property that $`\rho (p,q)`$ exponentially decreases to zero as $`p\pm \mathrm{}`$, we obtain

$`{\displaystyle \frac{dQ_h}{dt}}=2{\displaystyle \frac{\gamma }{m}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑p{\displaystyle _{l_h}^0}𝑑q\left({\displaystyle \frac{p^2}{2m}}{\displaystyle \frac{k_BT(q)}{2}}\right)\rho (p,q)`$

$$2\frac{\gamma }{m}\frac{p^2}{2m}\frac{k_BT_h}{2}_h.$$
(5)
This is the formula of the heat transfer from the hot bath to the particle . When the system is in equilibrium with the heat bath, the heat transfer $`\frac{dQ_h}{dt}`$ disappears, because the theory of equipartition requires $`<p^2/2m>=k_BT/2`$. This also shows that the energy exchange between the replica particle and the thermal bath is dominant only near the thermal transition point, $`q=0`$, where the average kinetic energy, $`p^2/2m`$, deviates from $`k_BT/2`$.

We remark here how the energy flows around the transition point. As we are analyzing the stalled state, the probability density $`\rho (p,q)`$ is stationary. Thus, the energy density, $`\rho _E(q)=_{\mathrm{}}^{\mathrm{}}𝑑p(p^2/2m+U(q))\rho (p,q),`$ is stationary. Because there is no work in the stalled state, the conservation of energy requires that $`\frac{d<Q_h+Q_c>}{dt}=0`$. It shows that the same quantity of the heat from hot bath to the particle flows from the particle to the cold bath: $`\frac{dQ_h}{dt}=\frac{dQ_c}{dt}`$ . It should also be remarked that, in the stalled state, the efficiency vanishes except that the both sides of the last equation vanishes, because the work, the numerator of the efficiency, is absent here. Quasi-static operation is reversible only if the last equation vanishes. Therefore the quantity $`\frac{dQ_h}{dt}`$ sufficiently characterizes the operation at the stalled state and thus we will analyze it in detail. Note that the following equality is simultaneously derived ,

$$\frac{dQ_h}{dt}=_{\mathrm{}}^{\mathrm{}}𝑑p\frac{p^2}{2m}\frac{p}{m}\rho (p,q)|_{q=0}.$$
(6)
The formula confirms that the irreversible heat transfer is carried microscopically as a kinetic energy of the particle at a transition point.

It is known, for example by the kinetic theory of gases, that there is finite heat transfer, $`I`$, in the system where a Brownian particle of finite mass and friction is crossing over the two regions with different temperature, even if the two thermal baths have no direct contact. It implies $`I\frac{dQ_h}{dt}>0`$. The authors of Ref. assumed that the heat transfer should disappear in overdamping limit, $`\frac{m}{\gamma }0`$. However, their assumption is not evident a priori. To reveal the validity of the assumption we have to perform proper energetic analysis on Kramers equation which includes the degree of momentum, $`p`$, instead of overdamped Fokker-Planck equation which lacks the degree.

Hereafter, we will consider the asymptotic behavior of the heat transfer, $`I`$, in the limit of overdamping process ($`\gamma +\mathrm{}`$ and/or $`m0`$). To find out the asymptotic behavior, it is convenient to use reference heat transfer $`I^{}`$ of unit mass and friction in an arbitrary set of units: $`I^{}I(m=1,\gamma =1)`$. By Eq.(5), the reference heat transfer $`I^{}`$ reads:

$`I^{}`$ $`=`$ $`2{\displaystyle \frac{p^2}{2}}{\displaystyle \frac{k_BT_h}{2}}_h`$ (7)
$`=`$ $`2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑p{\displaystyle _{l_h}^0}𝑑q\left({\displaystyle \frac{p^2}{2}}{\displaystyle \frac{k_BT_h}{2}}\right)\rho ^{}(p,q),`$ (8)
where $`\rho ^{}(p,q)`$ is a probability density in the reference state: $`m=1`$ and $`\gamma =1`$. We call a probability density $`\rho `$ and heat transfer of arbitrary mass and friction in a proper set of units as the generic probability density and the generic heat transfer. Note that the following result is not altered if we have a different reference state. The choice of the values $`m=1`$ and $`\gamma =1`$ for reference state is only for simplicity. In the reference state, the characteristic length of the transition region, $`l_{th}^{}`$, that the probability density in momentum, $`p`$, is out of equilibrium is: $`l_{th}^{}=v_{th}^{}\tau ^{}=\sqrt{k_BT}`$.

To evaluate the generic heat transfer (Eq.4) in terms of the reference heat transfer (Eq.4), we will find out the relation between the generic probability density with arbitrary mass and friction $`\rho (p,q)`$ and the reference one $`\rho ^{}(p,q)`$. The potential term $`K\frac{\rho }{p}`$ of Kramers equation, Eq.(1), is negligible when one discusses the asymptotic behavior of the overdamping. With stationary condition, $`\frac{}{t}=0`$, on Eq.(1), we obtain the simple equation that describes stationary flow in phase space around the boundary, $`q=0`$:

$$\frac{p}{m}\frac{\rho (p,q)}{q}=\frac{\gamma }{m}\frac{}{p}\left[p\rho (p,q)+mk_BT(q)\frac{\rho (p,q)}{p}\right].$$
(9)
We found here that this equation has a scaling property in mass and friction: The generic probability density $`\rho (p,q)`$ is expressed using the probability density of the reference state $`\rho ^{}(p,q)`$.

$$\rho (p,q)=c\rho ^{}(\frac{p}{\sqrt{m}},\frac{\gamma }{\sqrt{m}}q),$$
(10)
where the constant factor $`c`$ should be determined by normalization (Eq. 3).

As $`q`$ departs from the transition point, $`q=0`$, further than the characteristic length $`\mathrm{}_{th}`$, the probability density approaches that of equilibrium, where $`\rho _h(p,q)=C_h\mathrm{exp}\{\frac{p^2}{2mk_BT_h}\}`$ (for $`ql_{th}`$), and $`\rho _c(p,q)=C_c\mathrm{exp}\{\frac{p^2}{2mk_BT_c}\}`$ (for $`ql_{th}`$). The coefficients, $`C_h`$ and $`C_c`$, are then required to satisfy the condition of the continuity of probability current. Thus, we have

$$C_hT_h^{3/2}=C_cT_c^{3/2},$$
(11)
which is consistent with the condition derived for the overdamping limit. The remaining condition that determine $`C_x`$ is the normalization. Note that the normalization of the probability density, $`\rho `$, is satisfactorily carried out even by neglecting the contribution from the transition region because the characteristic scale of the transition region $`l_{th}`$ is much smaller than the width $`l_x`$: $`l_{th}/l_x1`$ ($`x=h`$ or $`c`$) (Eq.2). Then $`C_h`$ and $`C_c`$ are determined as:

$$C_h=\frac{1}{\sqrt{2\pi mk_BT_h}}\frac{T_c}{T_cl_h+T_hl_c},C_c=\frac{1}{\sqrt{2\pi mk_BT_c}}\frac{T_h}{T_cl_h+T_hl_c}.$$
(12)
With these solutions and Eq.(10), we obtain the relation between the two normalized probability densities $`\rho `$ and $`\rho ^{}`$:

$$\rho (p,q)=\frac{1}{\sqrt{m}}\rho ^{}(\frac{p}{\sqrt{m}},\frac{\gamma }{\sqrt{m}}q).$$
(13)
Note that this equation is valid even within the transition region.

We can now express the heat transfer, $`I`$, in terms of the reference heat transfer, $`I^{}`$. We rewrite $`I`$ as:

$`I`$ $`=`$ $`2{\displaystyle \frac{\gamma }{m}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑p{\displaystyle _{l_h}^0}𝑑q\left({\displaystyle \frac{p^2}{2m}}{\displaystyle \frac{k_BT_h}{2}}\right)\rho (p,q).`$ (14)
By change of variables such that $`p^{}=\frac{p}{\sqrt{m}}`$, $`q^{}=\frac{\gamma }{\sqrt{m}}q`$, we obtain

$`I`$ $`=`$ $`2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑p{\displaystyle _{\frac{\gamma l_h}{\sqrt{m}}}^0}𝑑q\left({\displaystyle \frac{p^2}{2}}{\displaystyle \frac{k_BT_h}{2}}\right)\rho (\sqrt{m}p,{\displaystyle \frac{\sqrt{m}q}{\gamma }}).`$ (15)
This yields using Eq.(13),

$`I=2{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑p{\displaystyle _{\frac{\gamma l_h}{\sqrt{m}}}^0}𝑑q\left({\displaystyle \frac{p^2}{2}}{\displaystyle \frac{k_BT_h}{2}}\right){\displaystyle \frac{1}{\sqrt{m}}}\rho ^{}(p,q).`$ (16)
This integrand is dominant only near the transition point, $`q=0`$, with characteristic length $`l_{th}^{}`$. As we are analyzing the asymptotic behavior such that $`m0`$ and/or $`\gamma \mathrm{}`$, the inequality, $`(l_{th}^{})l_h\gamma l_h/\sqrt{m}`$, is satisfied. Because the contribution from the interval $`q[\gamma l_h/\sqrt{m},l_h]`$ to the integral is negligible in Eq.(16) compared with that from $`q[l_h,0]`$, the interval of this integral may sufficiently be replaced by $`q[l_h,0]`$. Using Eq.(8), we obtain one of the main result of our paper:

$$I\frac{2}{\sqrt{m}}_{\mathrm{}}^{\mathrm{}}𝑑p_{l_h}^0𝑑q\left(\frac{p^2}{2}\frac{k_BT_h}{2}\right)\rho ^{}(p,q)=\frac{1}{\sqrt{m}}I^{}.$$
(17)
Because the characteristic length of the transition region vanishes in the overdamping limit, the scaling property is exact asymptotically.

From this result, we can learn that the irreversible heat transfer at the transition point does not decrease when one takes overdamping limit, which is contrast to the claim in Ref.: One way to take this limit is to increase the friction constant, $`\gamma `$: Then, the heat transfer does not decrease, because the heat transfer $`I`$ does not depend on $`\gamma `$. The other way is to decrease the mass, $`m`$: Then, the heat transfer does not decrease neither, moreover the heat transfer increases in the power of $`1/\sqrt{m}`$. The result justifies the intuitive estimation by Derényi and Astumian. The heat flow is a result of broken symmetry of probability density in momentum at the transition point, because the heat transfer disappears if the probability density is symmetric in phase space, as found by Eq.(6). Since an overdamped equation has no degree of freedom to describe the irreversible flow caused by the discontinuity of the temperature, the previous literature reached Carnot efficiency.

Up to now, we have discussed how heat transfer between the two heat baths behaves in the overdamping process. We found that the irreversible heat transfer does not decrease in the process. One finds, however, that the possible work out of the system may also vary according to the overdamping limit, because the probability current may vary due to the change of the parameters, $`\gamma `$ and $`m`$. Thus, it is not yet obvious whether non-vanishing heat transfer, $`I`$, itself reveals that the system cannot attain Carnot efficiency in any condition including non-stalled state. Thus, in addition to the irreversible heat transfer discussed above, we will estimate work and work-induced heat transfer in an overdamping process.

We will return to the original Kramers equation (Eq.1) for Büttiker-Landauer system. We have analyzed this equation with retaining both degrees of freedom, $`p`$ and $`q`$. However, we do not have to consider a momentum degree of freedom when we discuss work out of the Brownian system, because the work is the function only of the displacement of the position. Thus, we start the evaluation of the work by the overdamped Fokker-Planck equation for the probability density $`P(q)`$ of the system:

$$\frac{P(q)}{t}=\frac{}{q}J(q)=\frac{1}{\gamma }\frac{}{q}\left[\frac{U(q)}{q}+\frac{(k_BT(q))}{q}\right]P(q),$$
(18)
where the periodic boundary condition is applied: $`P(0)=P(L)`$ and $`\frac{dP}{dq}|_{q=0}=\frac{dP}{dq}|_{q=L}`$. Explicit mass dependence on the displacement of the system disappears in the overdamping limit. In stationary state, $`\frac{P(q)}{t}=0`$, the probability current, $`J(q)`$, is independent of $`q`$. The probability current $`J`$ reads

$$J=\frac{1}{\gamma }\left[\frac{U(q)}{q}+\frac{(k_BT(q))}{q}\right]P(q).$$
(19)
The equation for $`P(q)`$ reads

$$\frac{}{q}\left[\frac{U(q)}{q}+\frac{(k_BT(q))}{q}\right]P(q)=0.$$
(20)
This equation shows that the change of the friction constant, $`\gamma `$, does not alter the probability density, $`P(q)`$. Thus, with Eq.(19), the probability current, $`J`$, scales as: $`J\gamma ^1.`$ For a fixed load potential, the work per unit time $`\frac{dW}{dt}`$ is proportional to its probability current. Thus the sole operation, $`\gamma \mathrm{}`$, does not lead the system to Carnot efficiency, because the induced work $`(J)`$ decreases while the irreversible heat (Eq. (17)) does not decrease.

To find out the mass dependence on the work, we consider the working particle obeying the Stokes’ law with its radius $`r_B`$, where the mass and the friction are specified by one parameter, $`r_B`$: $`\gamma r_B`$ and $`mr_B^3`$. Thus, we have $`\frac{dW}{dt}J\gamma ^1r_B^1`$. The irreversible heat transfer $`\frac{dQ_{irr}}{dt}`$ that is independent of work is just the heat transfer, $`I`$ (Eq.(17)). Thus we have $`\frac{dQ_{irr}}{dt}m^{1/2}r_B^{3/2}.`$ The work-induced heat transfer $`Q_W`$ that is proportional to the work $`W`$ is proportional to the probability current $`J`$. Thus, we obtain $`\frac{dQ_W}{dt}Jr_B^1`$. The three components determines the efficiency. Thus, we have

$$\eta =\frac{\frac{dW}{dt}}{\frac{dQ_W}{dt}+\frac{dQ_{irr}}{dt}}=\frac{c_1r_B^1}{c_2r_B^1+c_3r_B^{3/2}}=\frac{c_1}{c_2+\frac{c_3}{\sqrt{r_B}}},$$
(21)
where $`c_1`$, $`c_2`$ and $`c_3`$ are constants. The result tells that the efficiency decreases monotonically to zero when one takes overdamping limit $`r_B0`$. This result is not altered even if one includes another transition point in the same period, because the asymptotic behavior of the two are the same.

In this paper, we have analyzed the energetics of a Brownian motor of Büttiker-Landauer type. We showed quantitatively that irreversible heat transfer does not disappear even if one takes overdamping limit ($`\gamma +\mathrm{}`$ and/or $`m0`$). This result is in contrast to the claims by Refs.. The mass dependence on the irreversible heat is consistent to the intuitive estimation by Ref.. We further analyzed the effect of non-vanishing irreversible heat transfer on the efficiency and showed that, even in fully overdamping limit, Carnot efficiency is inattainable for the particle obeying Stokes’ law. It shows that the maximum efficiency of the Brownian motor is not attained in the stalled state. The result revealed that the Brownian heat engine is qualitatively different from heat engines of which the most efficient operation is in quasi-static: Quasi-static process is the worst condition for the Brownian heat engine to work, while is the best for Carnot cycle.

The location of irreversible heat transfer is the transition region characterized by the thermal length $`l_{th}`$. It is certain that the characteristic length $`l_{th}`$ may disappear in the fully overdamping limit. Then, however, the irreversible effect at the transition region cannot be eliminated. From the result we also learn a lesson how to apply energetics to overdamping systems with a space-dependent temperature: We should apply energetics before taking overdamping limit. Otherwise, we might fail in proper evaluation of irreversible heat transfer within the transition region, because energetic interaction between the heat bath and the particle is carried by the momentum exchange between them. When the particle has smaller kinetic energy than that expected by equipartition theorem, the particle receive kinetic energy from the heat bath in average. Thus, if we lose the degree of momentum as in the overdamped equation, we cannot describe this existing physical process properly.

The present system cannot have maximum efficiency at a quasistatic condition. This means that the maximum efficiency is achieved with finite probability current, which is therefore accompanied by irreversible dissipation. Thus, the next challenging question would be, “Is there any principle that determines the optimal efficiency in Brownian heat engines?”

We would like to thank T. Chawanya, T. Mizuguchi, H. Hayakawa, H. Nishimori, M. Matsuo, S. Sasa, T. Shibata, Y. Hayakawa, M. Sano, F. Takagi, Y. Sawada, and T. Tsuzuki for helpful discussion. This work is supported in part by the Japanese Grant-in-Aid for Science Research Fund from the Ministry of Education, Science and Culture (Nos. 09740301, 11156216 and 12740226) and the Inamori Foundation.


format: text

text:
# Cosmic Statistics of Statistics: 𝑁-point Correlations

## 1. Introduction

Astrophysics provides prime examples of random fields, such as the distribution of galaxies and the fluctuations of the Cosmic Microwave Background (CMB). Such random fields can be characterized statistically by a series of well chosen estimators. The most popular ones include counts-in-cells, $`N`$-point correlation functions, as well as statistics derived from them. Indeed, there are mathematical theorems, which state that, under certain conditions, both series describe a random process fully.

Our goal in astrophysics is not simply estimating these statistics, but to constrain underlying theories. This aim necessitates a firm control over the expected statistical errors from a survey. The theory of errors for finite surveys, the “cosmic statistics of statistics”, is therefore of utmost importance for practical applications. Such a theory was formulated in the past for counts in cells statistics. \[9\],\[3\],\[2\]

For the $`N`$-point correlation functions, however, to date only partial results are published, such as calculation of the discreteness effects for the two-point correlation function\[7\], and for $`N`$-point correlation functions for Poisson and multinomial point processes\[10\], full calculation for the two-point function under the hierarchical assumption with edge effects neglected\[1\], and some results in Fourier space with various degree of approximations.

The aim of this proceedings is to present the general variance calculation for $`N`$-point correlation functions with all major contributions included, such as discreteness effects, arising from sampling the underlying random field with a finite number of galaxies, edge effects, due to the geometry of the survey and the corresponding uneven weighting of $`N`$-tuples, and finite volume effects, caused by fluctuations at and above the scale of a survey. The underlying technique of calculation, as well as the fully general results are presented here; specializations such as power spectrum and bispectrum, and approximations, such as weakly non-linear perturbation theory, hierarchical assumptions, will be presented elsewhere in collaboration with S. Colombi, and A.S. Szalay\[11\].

The next section sets up the general mathematical framework for the calculation using computer algebra packages. §3 presents the results for $`N=2`$, and $`N=3`$. The final discussion section outlines how the formulae can be used in practice, as well as describes developments in the immediate future.

## 2. General Framework

Many interesting statistics, such as the $`N`$-point correlation functions and their Fourier analogs, can be formulated as functions over $`N`$ points in a catalog. The covariance of a pair of such estimators will be calculated for a general point process under the assumption that the average density is a priori known. This is the obvious generalization of the Poisson process when arbitrarily high order correlations are present. The number of objects is thus varied corresponding to a grand canonical ensemble in statistical physics. The following calculations lean heavily on the elegant formalism by Ripley\[8\], which can be consulted for details, and are the direct generalization of the framework set up by Szapudi & Szalay.\[10\]

Let $`D`$ be a catalog of data points to be analyzed, and $`R`$ randomly generated over the same area, with averages $`\lambda `$, and $`\rho `$ respectively. The role of $`R`$ is to perform a Monte Carlo integration compensating for edge effects, therefore eventually the limit $`\rho \mathrm{}`$ will be taken. $`\lambda `$ on the other hand is assumed to be externally estimated with arbitrary precision. No other assumption is taken about the point process. In practice, $`\lambda `$ is usually estimated from the same survey, which gives rise to additional correlations, the “integral constraint bias”. This effect will be investigated in more detail elsewhere.

Following\[10\], let us define symbolically an estimator $`D^pR^q`$, with $`p+q=N`$ for a function $`\mathrm{\Phi }`$ symmetric in its arguments

$$D^pR^q=\mathrm{\Phi }(x_1,\mathrm{},x_p,y_1,\mathrm{},y_q),$$
(1)
with $`x_ix_jD,y_iy_jR`$. As an example, the two point correlation function corresponds to $`\mathrm{\Phi }(x,y)=[x,yD,rd(x,y)r+dr]`$, where $`d(x,y)`$ is the distance between the two points, and $`[condition]`$ equals $`1`$ when $`condition`$ holds, $`0`$ otherwise. Ensemble averages can be estimated via factorial moment measures, $`\nu _s`$.\[5\],\[8\] In the Poisson limit $`\nu _s=\lambda ^s\mu _s`$, where $`\mu _s`$ is the $`s`$ dimensional Lebesgue measure, and in the most general case $`\nu _sf(x_1,\mathrm{},x_s)\lambda ^s\mu _s`$. The function $`\lambda ^sf(x_1,\mathrm{},x_s)=F(x_1,\mathrm{},x_s)`$ can be identified as the full, i.e. non-reduced, $`N`$-point correlation function. The connection between these and the reduced $`N`$-point correlation functions is well known\[10\], and can be obtained through either generating functions, or recursions.

The general covariance of a pair of estimators is

$$E(p_1,p_2,N_1,N_2)=\widehat{D}_a^{p_1}\widehat{R}_a^{q_1}\widehat{D}_b^{p_2}\widehat{R}_b^{q_2}=\underset{i,j}{}\left(\genfrac{}{}{0pt}{}{p_1}{i}\right)\left(\genfrac{}{}{0pt}{}{p_2}{i}\right)i!\left(\genfrac{}{}{0pt}{}{q_1}{j}\right)\left(\genfrac{}{}{0pt}{}{q_2}{j}\right)j!S_{i+j}\lambda ^i\rho ^j,$$
(2)
with $`p_1+q_1=N_1,p_2+q_2=N_2`$, $`S`$ will be specified later, $`\widehat{}`$ denotes normalization with $`\lambda ,\rho `$ respectively, i.e. $`(\widehat{D}=D/\lambda ,\widehat{R}=R/\rho `$. The expression simply describes the fact that out of the $`p_1`$ and $`p_2`$ different data points in $`D`$ we have an $`i`$-fold degeneracy, as well as a $`j`$-fold degeneracy in the random points drawn from $`R`$. To simplify the exposition of the calculation, it is convenient to assume from the very beginning the $`\rho \mathrm{}`$, i.e. the random process employed for the Monte Carlo integration of the shape of the survey is arbitrarily dense. This is usually achievable in practice, thus it will not be considered further. The above equation describes the cross-correlation of two estimators even for two different objects as well: e.g., two particular bins of the two- and three-point correlation functions. An interesting special case, $`N_1=1`$ (the average density in the survey) and $`N_22`$, is needed for calculating the “integral constraint” correction.

When the random process is arbitrarily dense only $`j=0`$ survives,

$`E(p_1,p_2,N_1,N_2)=(\widehat{D}_a^{p_1}\widehat{R}_a^{q_1}\widehat{D}_b^{p_2}\widehat{R}_b^{q_2}=`$ (3)
$`{\displaystyle \underset{i,j}{}}\left({\displaystyle \genfrac{}{}{0pt}{}{p_1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{p_2}{i}}\right)i!\lambda ^i\widehat{S}_if(1,2,\mathrm{},p_1,N_1+1,\mathrm{},N_1+p_2i),`$
where $`\widehat{S}`$ is now an operator acting on $`f`$,

$$\widehat{S}_k=\mathrm{\Phi }_a(1\mathrm{}N_1)\mathrm{\Phi }_b(1\mathrm{}i,N_1+1,\mathrm{},N_1+N_2i)\mathrm{}\mu _{2Nk}.$$
(4)
The operator $`\widehat{S}_k`$ is analogous to the phase space integral $`S_k`$\[10\]. The dot emphasizes that the integral can be performed only after $`\widehat{S}_k`$ is acted on $`f`$ which is part of the measure. The phase space has to be calculated in the general case via the full factorial moment measure of which $`f`$ is an integral part. Throughout the paper we use the convention that $`\left(\genfrac{}{}{0pt}{}{k}{l}\right)`$ is nonzero only for $`k0,l0`$, and $`kl`$, and the variables $`x_i`$ are denoted with $`i`$ for simplicity. Here $`\mathrm{\Phi }_a`$ and $`\mathrm{\Phi }_b`$ denote two different functions, for instance corresponding to two radial bins of two estimators of the same or different orders. In the above formula the symmetry of $`\mathrm{\Phi }`$ in its arguments was heavily relied on to achieve the above “standard” representation of the operator.

The dependence of $`S_k`$ on $`a,b`$, and $`N`$ is not noted for convenience\[10\], but they will be assumed throughout the paper. The estimator\[10\] for the generalized $`N`$-point correlation function is $`(\widehat{D}\widehat{R})^N`$, or more precisely,

$$\stackrel{~}{w}_N=\frac{1}{S}\underset{i}{}\left(\genfrac{}{}{0pt}{}{N}{i}\right)()^{Ni}(\frac{D}{\lambda })^i(\frac{R}{\rho })^{Ni},$$
(5)
where $`S=\mathrm{\Phi }\mu _N`$ (without subscript). In this case $`S`$ is a number since it corresponds to the Poisson catalog with its simple factorial moment measure. The average of the estimator yields

$`\stackrel{~}{w}_N={\displaystyle \frac{1}{S}}{\displaystyle \underset{i}{}}\left({\displaystyle \genfrac{}{}{0pt}{}{N}{i}}\right)()^{Ni}{\displaystyle \mathrm{\Phi }(1,\mathrm{},N)f(1,\mathrm{},N)\mu _N}=`$ (6)
$`{\displaystyle \frac{1}{S}}{\displaystyle \mathrm{\Phi }(1,\mathrm{},N)\xi _N(1,\mathrm{},N)\mu _N}.`$
Since the role of $`\mathrm{\Phi }`$ is effectively a window, with a window operator $`\widehat{W}`$ this can be written symbolically as $`\stackrel{~}{w}_N=\widehat{W}\xi _N/\widehat{W}`$. The asymptotic centered covariance between two estimators of the above for a general point process in the limit of $`\rho \mathrm{}`$ can be written according to the previous considerations as

$`\delta \stackrel{~}{w}_{N_1}\delta \stackrel{~}{w}_{N_2}=\stackrel{~}{w}_{a,N_1}\stackrel{~}{w}_{b,N_2}\stackrel{~}{w}_{a,N_1}\stackrel{~}{w}_{b,N_2}=`$ (7)
$`{\displaystyle \frac{1}{S^2}}{\displaystyle \underset{i,j}{}}\left({\displaystyle \genfrac{}{}{0pt}{}{N_1}{i}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{N_2}{j}}\right)()^{i+j}\left[E(i,j,N_1,N_2)\widehat{S}_0f(1,\mathrm{},i)f(N_1+1,\mathrm{},N_1+j)\right].`$
In the above the operator $`S_0`$ acts on the multiple of the two $`f`$’s on the right. The above equation is the main result of this paper. While it is quite cumbersome, it is easily expandable with the help of computer algebra, as demonstrated by the examples of the next section. The special cases rendered will also illustrate how the simplicity of the proposed class of estimators exactly manifests itself by a mass cancellation of terms. Any other estimator would have extra terms in the variance.\[10\]

## 3. The cosmic error on the two-and three-point correlation functions

The above equation was entered into a computer algebra package. The symmetries and simplicity of the estimator give rise to cancellations and possibilities for collecting similar terms. This is the reason why the final result for the two-point correlation function has only three to six terms, while formally it could have up to about a hundred. Alternative estimators, such as $`DD/RR1`$, etc. would not yield significantly less cancellations, therefore error-calculation for them was not attempted; although the same formalism applies.

### 3.1. The two-point function

The covariance for the two-point function (or any quantity estimated from doublets of points, such as the power spectrum) can be expressed after the cancellations and the possible simplifications as

$`\delta \stackrel{~}{w}_2^a\delta \stackrel{~}{w}_2^b`$ $`={\displaystyle \frac{1}{S^2}}\{.{\displaystyle }\mathrm{\Phi }_a(1,2)\mathrm{\Phi }_b(3,4)[\xi _4(1,2,3,4)+2\xi (1,3)\xi (2,4)]+`$
$`{\displaystyle \frac{4}{\lambda }}{\displaystyle \mathrm{\Phi }_a(1,2)\mathrm{\Phi }_b(1,3)\left[\xi (2,3)+\xi _3(1,2,3)\right]}+`$
$`{\displaystyle \frac{2}{\lambda ^2}}{\displaystyle }\mathrm{\Phi }_a(1,2)\mathrm{\Phi }_b(1,2)[1+\xi (1,2)].\}.`$
The above equation is essentially identical to the result of Hamilton\[6\] where he calculates the variance of $`\delta `$, the fluctuation field itself. This is not at all surprising.\[10\] The estimator contains exactly the same terms and coefficients (regardless of the choice of $`\mathrm{\Phi }`$) as $`\delta `$ itself, which strongly suggests that it is (nearly) optimal.

The above formula contains the different contributions to the error\[9\] entirely mixed. Approximate separation of the different terms appears to be fruitless. The only general point to be made is that discreteness effects are absent in the first term, while they are present (mixed with the other two effects) in the $`1/\lambda ^s`$ terms, with $`s>0`$. This is naturally true for the higher order calculations as well.

It is worth to emphasize again, that this formula applies for the generalized $`2`$-point function, including the “traditional” $`2`$-point function, and any of its incarnations, such as the the power spectrum. In the latter case, esthetic or practical reasons might dictate that the errors on the power spectrum are expressed in terms of the power-spectrum, bi-, and tri-spectrum, instead of the two-, three-, and four-point correlation functions. Since the connection is a simple Fourier transform, this trivial exercise is left for the reader. Explicit formulae, aimed at practical applications for power-spectrum will be presented elsewhere.\[11\]

### 3.2. The three-point correlation function

The same method yields (co)variance for higher order estimators as well. We present another example, the generalized three-point correlation function. Its variance, using again the main result of the proceeding, translates into:

$`\delta \stackrel{~}{w}_3^a\delta \stackrel{~}{w}_3^b`$ $`={\displaystyle \frac{1}{S^2}}\{.`$
$`{\displaystyle }\mathrm{\Phi }_a(1,2,3)\mathrm{\Phi }_b(4,5,6)[.\xi (1,2,3,4,5,6)+`$
$`3\xi (1,2)\xi (3,4,5,6)+9\xi (1,4)\xi (2,3,5,6)+`$
$`3\xi (4,5)\xi (1,2,3,6)+9\xi (1,5,6)\xi (2,3,4)+`$
$`9\xi (1,4)\xi (2,3)\xi (5,6)+6\xi (1,4)\xi (2,5)\xi (3,6).]`$
$`{\displaystyle \frac{9}{\lambda }}{\displaystyle }\mathrm{\Phi }_a(1,2,3)\mathrm{\Phi }_b(1,4,5)[.\xi (1,2,3,4,5)+`$
$`\xi (2,3,4,5)+2\xi (1,2)\xi (3,4,5)+`$
$`2\xi (1,4)\xi (2,3,5)+\xi (2,3)\xi (1,4,5)+`$
$`4\xi (2,5)\xi (1,3,4)+\xi (4,5)\xi (1,2,3)+`$
$`\xi (2,3)\xi (4,5)+2\xi (2,4)\xi (3,5).]+`$
$`{\displaystyle \frac{18}{\lambda ^2}}{\displaystyle }\mathrm{\Phi }_a(1,2,3)\mathrm{\Phi }_b(1,2,4)[.\xi (1,2,3,4)+`$
$`2\xi (1,3,4)+\xi (1,2)\xi (3,4)+`$
$`2\xi (1,3)\xi (2,4)+\xi (3,4).]`$
$`{\displaystyle \frac{6}{\lambda ^3}}{\displaystyle }\mathrm{\Phi }_a(1,2,3)\mathrm{\Phi }_b(1,2,3)[.\xi (1,2,3)+3\xi (1,2)+1.].\}.`$
For simplicity, in the above formula the order of $`\xi `$ is notated with the number of arguments only, e.g., $`\xi _3(1,2,3)=\xi (1,2,3)`$. The above equation is is less obviously useful then that of the two-point correlation function. Nevertheless, given a model for the higher order correlation functions, such as weakly non-linear perturbation theory, or any well specified version of the hierarchical assumption, the equation can easily be turned into a practical computer program.\[11\]

The variance of the four-point and higher order correlation functions can be calculated analogously, but it would make no sense to print the results. When needed, the formulae generated by computer algebra can be transformed into Fortran or C-code directly.

## 4. Discussion

The above method, and the explicit formulae given, can be used to evaluate the cosmic error on any statistical measure based on $`N`$-tuples of a distribution. This includes, but is not limited to, $`N`$-point correlation functions, $`N`$-th order cumulants, cumulant correlators, multi-spectra, etc.

The above calculation was performed only for the best $`N`$-point estimator.\[10\] Any other estimator can be calculated analogously, but be warned that the resulting number of terms can be overwhelming due to the insufficient cancellation arising from suboptimal edge correction.

The fact that the average density is assumed to be given as an outside parameter appears to be fairly restrictive. However, maximum likelihood context, which is the most important potential practical application of the results, it is easy to remedy the situation. The proposed estimator\[10\] can be trivially changed by not normalizing with the average density $`\lambda `$. This introduces only a small modification to the final formulae, and a set of estimators, including that of the the average count, contains all information need for constructing likelihood function. Such a procedure yields full statistical description, takes into account fluctuations in the mean, and the fact the average is estimated from the same surveys (“integral constraint”). Practical demonstration of this procedure will be presented elsewhere.\[11\]

The proposed estimator used here is not connected for $`N4`$.\[10\] Therefore the calculations for the higher order connected estimator should be modified for accurate results for the connected moments. This trivial but tedious generalization is left for future research.

The explicit formulae can be specialized for several cases, which will be presented elsewhere.\[11\] The interesting limits include Poisson, Gaussian, weakly-nonlinear, strong correlations, hierarchy, shot noise limited, continuous limit etc. The detailed discussion of these limits, and specializations to particular statistics, such as $`N`$-point correlation functions, multi-spectra, would go beyond the scope of the present exposition. Similarly, the main equation yields cross correlations between different statistics as well, a must for any investigation in the maximum likelihood framework.

Finally, it is worth to note here, that recent advances in algorithms for calculating $`N`$-point correlation functions render these objects more interesting then ever. Fast algorithms\[4\] will make it possible to calculate $`N`$-point functions from very large catalogs, be it artificial or real data, which undoubtedly will culminate in new insights into the subject. The formulae presented in this proceedings will provide the firm theoretical grounding for any such investigation.

Acknowledgments

This work will form the integral part of a project in collaboration with S. Colombi, and A.S. Szalay soon to be published.\[11\]

## 5. An Alternative Technique

An alternative method of calculation is possible, which is instructive and insightful, even if less rigorous mathematically then the above formalism using factorial moment measures. This alternative technique is well suited for obtaining quick results for low order estimators by hand. We demonstrate the calculation for two-point correlation function, higher order results can be obtained analogously, although it quickly becomes prohibitevely tedious.

Let us assume that the survey is divided into $`K`$ bins, each of them with fluctuations $`\delta _i`$, with $`i`$ running from $`1\mathrm{}K`$. For this configuration our estimator can be written as

$$\stackrel{~}{w}=f_{12}\delta _1\delta _2.$$
(10)
The above equation uses a “shorthand” Einstein convention: $`1,2`$ substituting for $`i_1,i_2`$, and repeated indices summed. The pairwise weights $`f_{12}`$ correspond to $`\mathrm{\Phi }`$ in the main body of the paper, and it is assumed that the two indices cannot overlap.

The ensemble average of the above estimator is clearly $`f_{12}\xi _{12}`$. The continuum limit (co)variance between bins $`a`$ and $`b`$ can be calculated by taking the square of the above, and taking the ensemble average.

$$\delta \stackrel{~}{w}^a\delta \stackrel{~}{w}^b=f_{12}^af_{34}^b\left(\delta _1\delta _2\delta _3\delta _4\delta _1\delta _2\delta _3\delta _4\right).$$
(11)
Note that the averages in this formula are not connected moments, which are distinguished by $`_c`$.

The above equation yields only the continuum limit terms. To add Poisson noise contribution to the error, note that it arises from the possible overlaps between the indices (indices between two pairweights can still overlap!). In the spirit of infinitesimal Poisson models, we replace each overlap with a $`1/\lambda `$ term, and express the results in terms of connected moments. There are three possibilities, i) no overlap (continuum limit)

$$f_{12}^af_{34}^b\left(\xi _{1234}+\xi _{13}\xi _{24}+\xi _{14}\xi _{23}\right),$$
(12)
ii) one overlap (4 possibilities)

$$\frac{4}{\lambda }f_{12}^af_{13}^b\left(\xi _{123}+\xi _{23}\right),$$
(13)
iii) two overlaps (2 possibilities)

$$\frac{2}{\lambda ^2}f_{12}^af_{12}^b\left(1+\xi _{12}\right),$$
(14)
In these equations, for the sake of the Einstein convention we used $`\xi (i,j,k,l)\xi _{ijkl}`$. The above substitutions (rigorously true only in the infinitesimal Poisson sampling limit) become increasingly accurate with decreasing cell size. In that limit, adding the above three equations is equivalent to Eq. 3.1..


text:
# Untitled Document

FTUAM–00-18;

QCD Calculations of Heavy Quarkonium States **Talk given at the conference “ Hyperons, Charm and Beauty Hadrons”, Valencia, June 2000. To be published in Nucl. Phys. Suppl.

Typeset with PHys MaTEX

F. J. Ynduráin Departamento de Física Teórica, C-XI, Universidad Autónoma de Madrid, Canto Blanco, 28049-Madrid, Spain. E-mail: fjy@delta.ft.uam.es

Recent results on the QCD analysis of bound states of heavy $`\overline{q}q`$ quarks are reviewed, paying attention to what can be derived from the theory with a reasonable degree of rigour. We report a calculation of $`\overline{b}c`$ bound states; a very precise evaluation of $`b,c`$ quark masses from quarkonium spectrum; the NNLO evaluation of $`\mathrm{{\rm Y}}e^+e^{}`$; and a discussion of power corrections. For the $`b`$ quark pole mass we get, including $`O(m_c^2/m_b^2)`$ and $`O(\alpha _s^5\mathrm{log}\alpha _s)`$ corrections, $`m_b=5.020\pm 0.058GeV`$; and for the $`\overline{MS}`$ mass the result, correct to $`O(\alpha _s^3)`$, $`O(m_c^2/m_b^2)`$, $`\overline{m}_b(\overline{m}_b)=4.286\pm 0.036GeV`$. For the decay $`\mathrm{{\rm Y}}e^+e^{}`$, higher corrections are too large to permit a reliable calculation, but we can predict a toponium width of $`13\pm 1keV`$.

1. INTRODUCTION In the present note we are going to review some aspects of the QCD analysis of heavy quarkonia, $`\overline{t}t`$, $`\overline{c}c`$, $`b\overline{c}`$ and especially $`\overline{b}b`$ states. The validity of ab initio QCD calculations varies from one case to another. For the energy levels of toponium with $`n5`$, $`n`$ being the principal quantum number, and for the energy of the ground state in bottomium, we have a very satisfactory situation, almost comparable to that in positronium calculations. For bottomium hyperfine splitting, and for the energy of the charmonium ground state, the situation is less favourable, although reasonably under control. Less favourable still is what one has for the energy levels of bottomium with $`n=2`$. Here NLO and NNLO perturbative and the leading nonperturbative corrections are comparable to the LO results. The same is true for the decay $`\mathrm{{\rm Y}}e^+e^{}`$, and even for toponium decay, $`Te^+e^{}`$, higher corrections are a bit too large for comfort. For higher states, for the decays of $`n=2`$ bottomium or for the ground state of charmonium, a rigorous QCD evaluation is out of the question. There are methods that have been devised in the literature to deal with this, essentially more or less justified models; we do not discuss them now. The interested reader may find information, and references, in e.g. my Lisbon lectures.

This does not mean that even in the favourable cases all questions are settled, and in particular we discuss here power-like corrections, at short distances, where there is still some controversy. BRRR 2. NONRELATIVISTIC QCD POTENTIAL, AND CORRECTIONS To analyze the lowest states of heavy quarks one proceeds as follows: first, in the nonrelativistic (NR) approximation one finds the potential that governs the dynamics. This potential is evaluated to increasing orders of perturbation theory; at present it is known to one and two loops. One can thus write

$$H=H^{(0)}+H_{\mathrm{log}};$$
$`(2.1\mathrm{a})`$
where $`H^{(0)}`$ may be solved exactly and contains the coulombic part of the interaction. For $`\overline{q}q`$ states,

$$H^{(0)}=2m+\frac{1}{m}\mathrm{\Delta }\frac{C_F\stackrel{~}{\alpha }_s(\mu ^2)}{r},$$
$`(2.1\mathrm{b})`$
$$\begin{array}{cc}\hfill \stackrel{~}{\alpha }_s(\mu ^2)=& \alpha _s(\mu ^2)\left\{1+c^{(1)}\frac{\alpha _s(\mu ^2)}{\pi }+c^{(2)}\frac{\alpha _s^2}{\pi ^2}\right\};\hfill \\ \hfill c^{(1)}=& a_1+\frac{\gamma _\mathrm{E}\beta _0}{2}\hfill \\ \hfill c^{(2)}=& \gamma _\mathrm{E}\left(a_1\beta _0+\frac{\beta _1}{8}\right)+\left(\frac{\pi ^2}{12}+\gamma _\mathrm{E}^2\right)\frac{\beta _0}{4}+a_2.\hfill \end{array}$$

The expression for the $`a_i`$ may be found in e.g. ref. 1, $`m`$ is the pole mass of the quark, and $`H_{\mathrm{log}}`$ is

$$\begin{array}{cc}\hfill H_{\mathrm{log}}=& \frac{C_F\beta _0\alpha _s(\mu ^2)^2}{2\pi }\frac{\mathrm{log}r\mu }{r}+\frac{C_F\beta _0^2\alpha _s^3}{4\pi ^2}\frac{\mathrm{log}^2r\mu }{r}\hfill \\ \hfill +& \frac{C_F\alpha _s^3}{\pi ^2}\left(a_1\beta _0+\frac{\beta _1}{8}+\frac{\gamma _\mathrm{E}\beta _0^2}{2}\right)\frac{\mathrm{log}r\mu }{r}.\hfill \end{array}$$

One solves exactly the coulombic Schrödinger equation,

$$H^{(0)}\mathrm{\Psi }_{nl}=E_{nl}^{(0)}\mathrm{\Psi }_{nl},$$
$`(2.2)`$
treating the terms in $`H_{\mathrm{log}}`$ as perturbations. Relativistic corrections are then added, to be also treated as perturbations. These relativistic corrections include spin-dependent ones, known at tree level since a long time (they are essentially like for positronium) and to one loop from refs. 3,4. Moreover, and because of the nonabelian character of QCD, we have mixed relativistic–one loop corrections for the spin-independent piece, first evaluated in ref. 5 using the method of equivalent potentials of Gupta et al. (see e.g. ref. 4). Concentrating on the spin-independent hamiltonian we then write

$$H_{\mathrm{SI}}^{\mathrm{e}.\mathrm{p}.}=H^{(0)}+H_{\mathrm{log}}+\frac{1}{4m^3}\mathrm{\Delta }^2+\frac{C_F\alpha _s}{m^2r}\mathrm{\Delta }+\frac{C_Fb_1\alpha _s^2}{2mr^2}$$
$`(2.3)`$
and $`b_1=\frac{1}{2}C_F2C_A`$. Here the superindex “e.p.” means that the hamiltonian was obtained with the method of equivalent potentials.

The hamiltonian has been re-evaluated recently by Brambilla et al.; these authors generalize the nonrelativistic (heavy quark) effective theory (NRQCD), matching then the Green’s functions to a potential description (pNRQCD). In this way they find,

$$\begin{array}{cc}\hfill H_{\mathrm{SI}}^{\mathrm{pHQET}}=H^{(0)}+H_{\mathrm{log}}+\frac{1}{4m^3}\mathrm{\Delta }^2& \\ \hfill +\frac{\pi C_F\alpha _s}{m^2}\delta (𝐫)+\frac{C_F\alpha _s}{2r^3m^2}𝐋^2\frac{C_F\alpha _s}{2m^2}\{\frac{1}{r},\mathrm{\Delta }\}& \\ \hfill +\frac{C_FC_A\alpha _s^2}{2mr^2}& \\ \hfill .& \end{array}$$
$`(2.4)`$
This is different from (2.3). The difference is due to the fact that in the derivation of (2.3) one works with the S-matrix, while for (2.4) one uses Green’s functions. The difference between (2.3) and (2.4) vanishes when taking expectation values between coulombic wave functions, say, solutions of (2.2): so they will produce the same energy spectrum, at least to order $`\alpha _s^4`$. At higher orders, (2.3) and (2.4) would produce different results; but also new terms will appear in the two formalisms: their equivalence, or non-equivalence, has not been proved. For numerical results to $`O(\alpha _s^4)`$ for all states see ref. 9; here we present those for the ground state energy. We write

$$E_{10}^{\mathrm{p}.\mathrm{t}.}=2mm\frac{C_F^2\stackrel{~}{\alpha }_s^2}{4}+\delta _{\mathrm{p}.\mathrm{t}.}E_{10}.$$
$`(2.5\mathrm{a})`$
The label “p.t.” in $`E_{nl}^{\mathrm{p}.\mathrm{t}.}`$ indicates that we have as yet only used results deduced from perturbation theory; the full expression would be

$$E_{nl}=E_{nl}^{\mathrm{p}.\mathrm{t}.}+\delta _{\mathrm{NP}}E_{nl}.$$

$`\delta _{\mathrm{NP}}E_{nl}`$ embodies the nonperturbative contributions, to be discussed below. The $`\delta _{\mathrm{p}.\mathrm{t}.}E_{10}`$ is, with $`a=2/mC_F\stackrel{~}{\alpha }_s`$,

$$\begin{array}{cc}\hfill \delta _{\mathrm{p}.\mathrm{t}.}E_{10}=& \frac{10}{8m^3a^4}\frac{3C_F\alpha _s}{m^2a^3}+\frac{C_Fb_1\alpha _s^2}{ma^2}\hfill \\ & \frac{C_Fc_2^{(L)}\alpha _s^3}{\pi ^2a}\left[\mathrm{log}\frac{a\mu }{2}+\psi (2)\right]\hfill \\ & \frac{C_F\beta _0^2\alpha _s^3}{4\pi ^2a}\{\mathrm{log}^2\frac{a\mu }{2}\frac{\gamma _\mathrm{E}}{2}\mathrm{log}\frac{a\mu }{2}\hfill \\ \hfill +& \frac{3+3\gamma _\mathrm{E}^2\pi ^2+6\zeta _3}{12}\},\hfill \end{array}$$
$`(2.5\mathrm{b})`$
$`c_2^{(L)}=a_1\beta _0+\frac{1}{8}\beta _1+\frac{1}{2}\gamma _\mathrm{E}\beta _0^2.`$

The (nominally) leading, $`\mathrm{log}\alpha _s`$ corrections of next order are also known. They include a logarithmic correction to the static potential, and a relativistic one-loop correction. With the full result (as given in ref. 6) one has, for $`\mu =2/a`$,

$$\delta _{[\alpha _s^5\mathrm{log}\alpha _s]}E_{10}=m[C_F+\frac{3}{2}C_A]C_F^4\alpha _s^5(\mathrm{log}\alpha _s)/\pi .$$
$`(2.5\mathrm{c})`$
In the calculations we will include, for the $`b`$ quark case, a correction of order $`m_c^2/m_b^2`$,

$$\delta _{[m_c^2/m_b^2]}E_{10}=2m_b\frac{3T_F\alpha _s}{2\pi }\frac{m_c^2}{m_b^2}.$$
$`(2.5\mathrm{d})`$
We can invert (2.5) to obtain a precise determination of $`c,b`$ quark masses from those of the $`J/\psi ,\mathrm{{\rm Y}}`$ particles. The determination includes nonperturbative effects, to be discussed below. As input parameters we take the recent determination,

$$\begin{array}{cc}\hfill \mathrm{\Lambda }(n_f=4,\text{three loops})=0.283\pm 0.035GeV& \\ \hfill \left[\alpha _s(M_Z^2)0.117\pm 0.024\right],& \end{array}$$

and for the gluon condensate the value $`\alpha _sG^2=0.06\pm 0.02GeV^4.`$

Then comes the matter of the renormalization point, $`\mu `$. As Eqs. (2.5b,c) show, logarithms are avoided until order $`\alpha _s^5`$ if choosing $`\mu =\mu _0=2/a`$, and this will be our central value here. Then, for the $`b`$ quark, we find

$$\begin{array}{cc}\hfill m_b=& 5.020\pm 0.043(\mathrm{\Lambda })0.005(\alpha _sG^2)\hfill \\ & {}_{+0.037}{}^{0.031}(\text{vary}\mu ^2\mathrm{by}\mathrm{\hspace{0.17em}25}\%)\hfill \\ \hfill \pm & 0.006\pm 0.017(\mathrm{other}\mathrm{th}.\mathrm{uncert}.)\hfill \\ \hfill =& 5.020\pm 0.058GeV.\hfill \end{array}$$
$`(2.6\mathrm{a})`$
To obtain this, we have included the correction of Eq. (2.5c), and for the “theoretical error” also that estimated in ref. 13 for the perturbative evaluations. The value of the $`\overline{\text{MS}}`$ mass that corresponds to this is, taking into account $`O(\alpha _s^3)`$ and leading $`O(m_c^2/m_b^2)`$ corrections,

$$\overline{m}_b(\overline{m}_b)=4272\pm 43MeV.$$
$`(2.6\mathrm{b})`$
We present a summary of QCD calculations of the quark masses, with increasing accuracy.

Ref. $`m_b(\mathrm{pole})`$ $`\overline{m}_b(\overline{m}_b^2)`$ $`m_c(\mathrm{pole})`$                 TY $`4971\pm 72`$ $`4401_{35}^{+21}`$$`4_{4_4}^{4^4}`$ $`1585\pm 20^{}`$   PY $`5065\pm 60`$ $`4455_{29}^{+45}`$$`4_{4_4}^{4^4}`$ $`1866_{133}^{+215}`$   Here $`5022\pm 58`$ $`4272\pm 43`$ $``$          $`b`$, $`c`$ quark masses. $`(^{})`$ Systematic errors not included.    TY: Titard and Ynduráin. $`O\left(\alpha _s^3\right)`$ plus $`O\left(\alpha _s^3\right)v`$, $`O\left(v^2\right)`$ for $`m`$; $`O\left(\alpha _s^2\right)`$ for $`\overline{m}`$. Rescaled for $`\mathrm{\Lambda }\left(n_f=4\right)=283MeV`$. PY: Pineda and Ynduráin. Full $`O\left(\alpha _s^4\right)`$ for $`m`$; $`O\left(\alpha _s^2\right)`$ for $`\overline{m}`$. Rescaled for $`\mathrm{\Lambda }\left(n_f=4\right)=283MeV`$. Here: This calculation. $`O\left(\alpha _s^4\right)`$, $`O\left(\alpha _sm_c^2/m_b^2\right)`$ and $`O\left(\alpha _s^5\mathrm{log}\alpha _s\right)`$ for $`m`$; $`O\left(\alpha _s^3\right)`$ and $`O\left(\alpha _s^2m_c^2/m_b^2\right)`$ for $`\overline{m}`$. Values not given for the $`\overline{\text{MS}}`$ $`c`$ quark mass, as the higher order terms are as large as the leading ones.

It is to be noted that pure QCD quarkonium determinations of $`m_b`$ have remained remarkably stable since the earlier evaluations (ref. 5). \[The variation of the published values is due to the variation of the favoured central value of $`\mathrm{\Lambda }`$, from $`200`$ to $`283`$ MeV\]. Also, the results depend very little on the value of $`\mu `$ chosen. Thus, even taking $`\mu =m5GeV`$ one gets

$$m_b=4901\pm 60MeV,$$

reasonably close to (2.6a), where we had elected $`\mu =2/a2.88GeV`$. The estimates of $`\overline{m}_b(\overline{m}_b)`$ are less stable, doubtlessly because of the large size of the $`O(\alpha _s^3)`$ corrections; but it is interesting that this more precise spectroscopic determination of $`\overline{\text{MS}}`$ mass of the $`b`$ quark mass given here agrees with recent determinations, based on sum rules, also accurate to $`O(\alpha _s^3)`$. These give

$$\overline{m}_b(\overline{m}_b)=4260\pm 100MeV.$$

We will discuss $`\overline{\text{MS}}`$ masses further later on.

The hyperfine splitting can likewise be evaluated, getting a prediction for the $`\eta _b`$ mass:

$$M(\mathrm{{\rm Y}})M(\eta _b)=53.3\pm 12.5MeV.$$

The methods discussed for quarkonium can be extended with some work to hybrid states like $`b\overline{c}`$. The advantage here is that, since the masses of $`b,c`$ quarks can be expressed in terms of those of the $`\mathrm{{\rm Y}},J/\psi `$ particles, the theoretical errors are smaller than what one would get for charmonium. A recent $`O(\alpha _s^4)`$ calculation gives

$$M_{10}(b\overline{c})=6323\pm 10\pm 20MeV,$$

the first error being perturbative, the second coming from estimated nonperturbative effects. 3. POWER CORRECTIONS. THE NONPERTURBATIVE VACUUM; RENORMALONS; SATURATION As stated above, a calculation such as that in (2.6) includes leading nonperturbative effects. These are obtained by realizing that the motion of the $`\overline{q}q`$ pair takes place in the physical vacuum, chock full of soft gluons and light quark pairs so that, in particular, we have a nonzero value for the gluon condensate:

$$\mathrm{vac}|:\alpha _sG^2(0):|\mathrm{vac}\alpha _sG^20.$$

The effects of this were first considered by Leutwyler and Voloshin (see also refs.. 18, 5) and amount to a shift for (say) the ground state energy of quarkonium of

$$\delta _{\mathrm{NP}}^{\alpha _sG^2}=m\frac{\pi ϵ_{10}\alpha _s:G^2:}{(mC_F\alpha _s)^4},$$
$`(3.1)`$
$`ϵ_{10}1.5`$. This is zero to all orders of perturbation theory, because so is the gluon condensate. On dimensional grounds one expects $`\alpha _sG^2\mathrm{\Lambda }^4`$, so that (3.1) is of order $`\mathrm{\Lambda }^4/m^3\alpha _s^4`$.

This is not the only power correction that may appear. We may have renormalons. Let us consider a one-gluon exchange diagram, for $`\overline{q}q`$ scattering. If we dress the gluon propagator with loops then the corresponding potential, in momentum space, is

$$\stackrel{~}{V}(k)=\frac{4\pi C_F}{k^2}\frac{4\pi }{\beta _0\mathrm{log}(k^2/\mathrm{\Lambda }^2)},$$
$`(3.2)`$
and we have substituted the one-loop expression for $`\alpha _s(k^2)`$. (3.2) is undefined for soft gluons, with $`k^2\mathrm{\Lambda }^2`$. As follows from the general theory of singular functions, the ambiguity is of the form $`c\delta (k^2\mathrm{\Lambda }^2)`$: upon Fourier transformation this produces an ambiguity in the $`x`$-space potential of $`\delta V(r)=c[\mathrm{sin}\mathrm{\Lambda }r]/r`$. At short distances we expand this in powers of $`r`$ and find

$$\delta _{\mathrm{renorm}.}V(r)C_0+C_1r^2+\mathrm{}.$$
$`(3.3)`$
The same result may be obtained with the more traditional method of Borel transforms. This coincides with the short distance behaviour of the nonperturbative potential as determined by Dosch and Simonov.

The situation just described applies for states $`\overline{q}q`$ at short distances; but not so short that zero frequency gluons cannot separate the pair. If this last is the case, soft gluons do not resolve the $`\overline{q}q`$ pair and only see a dipole. The generated renormalon may then be seen (ref. 20) to correspond to the contribution of the gluon condensate in the Leutwyler-Voloshin mechanism.

What happens to the Aglietti and Ligeti renormalon? In fact, we have other long-distance power corrections. It is clear that the pole mass is defined purely in perturbation theory, and indeed one can check that a renormalon ambiguity appears already at one loop. Likewise, the coulombic potential is defined so that it vanishes at infinity: but, for confined quarks, “infinity” is equivalent to the confinement radius, $`R1/\mathrm{\Lambda }`$. Actually, mass and potential renormalons cancel for the constant and quadratic terms in (3.2), in the very short distance regime; this has been verified by a detailed calculation in ref. 22, thus substantiating the intuitive arguments of ref. 20, already mentioned.

figure 1.  “Filled in” gluon loop.

Renormalons and the mechanism of Leutwyler–Voloshin are not the only possible sources of power corrections; we next consider saturation. We note that the ambiguities we have found are associated with small momenta or, equivalently, long distances. However, the singularities are clearly spurious. Indeed, not only the theory should be well defined but, because of confinement, long distances are never attained: the theory possesses an internal infrared cut-off of the order of the confinement radius, $`R\mathrm{\Lambda }^1`$. To try and implement it we consider the gluon propagator. To one loop it gets a correction involving the vacuum polarization tensor. Neglecting quarks this is, in $`x`$-space, given by an expression like

$$\begin{array}{cc}\hfill \mathrm{\Pi }_{\alpha \beta }^{aa^{}}(x,0)& g^2f_{abc}f_{a^{}de}0|\mathrm{d}^4y_1\mathrm{d}^4y_2\hfill \\ \hfill \mathrm{T}B_b^\alpha (y_1)_\mu & B_{c\alpha }(y_1)B_d^\beta (y_2)_\nu B_{e\beta }(y_2)|0\hfill \\ \hfill +& \mathrm{}.\hfill \end{array}$$

We can take into account the long distance interactions by introducing a string between the field products at finite distances. In matrix notation for the gluonic fields, $`^\mu =t_aB_a^\mu `$, this is incorporated by replacing

$$\begin{array}{cc}\hfill ^\alpha (y_1)^\beta (y_2)& \\ \hfill ^\alpha (y_1)\mathrm{P}\left(\mathrm{exp}\mathrm{i}_{y_2}^{y_1}dz^\mu _\mu (z)\right)^\beta (y_2).& \end{array}$$

The process may be described as “filling the loop” (see Fig. 1) by introducing all exchanges between the gluonic lines there. If we furthermore replace the perturbative vacuum $`|0`$ by the nonperturbative one $`|\mathrm{vac}`$, then we get a dressed propagator

$$D_{\mathrm{dressed}}^{\mu \nu }(k)=D^{(0)\mu \nu }(k)\frac{4\pi }{\beta _0\mathrm{log}(M^2+k^2)/\mathrm{\Lambda }^2},$$

and $`M^2`$ is related to the gluon condensate at finite distances, $`G(x)G(0)_{\mathrm{vac}}`$. (for more details and references, see e.g. Simonov’s lectures).

This suggests a saturation property of the coupling constant at small momenta (long distances) so that the expression for the running coupling constant should be modified according to

$$\begin{array}{cc}\hfill \alpha _s(k^2)=& \frac{4\pi }{\beta _0\mathrm{log}k^2/\mathrm{\Lambda }^2}\hfill \\ \hfill & \alpha _s^{\mathrm{sat}}(k^2)=\frac{4\pi }{\beta _0\mathrm{log}(k^2+M^2)/\mathrm{\Lambda }^2}.\hfill \end{array}$$
$`(3.4)`$
It is certain that an expression such as (3.4) incorporates, to some extent, long distance properties of the QCD interaction. For example, if we take (3.4) with $`M=\mathrm{\Lambda }`$ in the tree level potential for heavy quarks, this becomes the Richardson potential

$$\begin{array}{cc}\hfill V^{(0)}(𝐤)=& \frac{4\pi C_F\alpha _s(𝐤^2)}{𝐤^2}\hfill \\ \hfill & V^{(0),\mathrm{sat}}(𝐤)=\frac{16\pi ^2C_F}{\beta _0𝐤^2\mathrm{log}(k^2+\mathrm{\Lambda }^2)/\mathrm{\Lambda }^2}.\hfill \end{array}$$

When one has $`𝐤^2\mathrm{\Lambda }^2`$, the short distance coulombic potential is, of course, recovered. For $`𝐤^2\mathrm{\Lambda }^2`$, however,

$$V^{(0),\mathrm{sat}}(𝐤)\underset{𝐤^2\mathrm{\Lambda }^2}{}\frac{16\pi ^2C_F\mathrm{\Lambda }^2}{\beta _0𝐤^4},$$

whose Fourier transform gives

$$V^{(0),\mathrm{sat}}(r)\underset{r\mathrm{\Lambda }^1}{}(\text{constant})\times r,$$

i.e., a linear potential. Indeed, a reasonably accurate description of spin-independent splittings in quarkonia states is obtained with such a potential. However, this long-distance linear potential induced by saturation in the Richardson model is the fourth component of a Lorentz vector, while we know that the Wilson linear potential, as obtained, e. g., in the stochastic vacuum model or in lattice calculations, should be a Lorentz scalar. This is also obvious on phenomenological grounds as this potential has to provide attraction for quark-quark states in e.g. baryons. It thus follows that the linear potential obtained from saturation can be only of limited phenomenological use.

Saturation also gives a linear correction to the short distance coulombic interaction. The possibility of such a correction has been discussed by several people; see, for example refs.. 20, 23. Personally I am not impressed by these arguments. It is clear that the QCD perturbative series is not convergent; in the case of quarkonia, the coefficients are not even analytic as one has terms in $`\mathrm{log}\alpha _s`$ (cf. (2.5c) and Sect. 4 below). The results one finds are then dependent on how one sums and, moreover, it is not guaranteed that results valid for large orders of perturbation theory will hold in the real world, where one knows at most three nontrivial terms. It is the author’s belief that only if the summation method is rooted on solid physics it is likely to represent an improvement; otherwise, estimates of nonperturbative effects become pure guesswork. In this respect, the method of taking into account the nonperturbative nature of the physical vacuum by considering the effect of nonzero values for the correlators stands some chance of being meaningful, as indeed phenomenological calculations seem to indicate. 4. FURTHER DISCUSSION OF THE $`\overline{\text{MS}}`$ MASS: A MATTER OF CONVERGENCE Write, for a heavy quark,

$$\begin{array}{cc}\hfill \overline{m}(\overline{m})& m/\{1+\delta _1+\delta _2+\delta _3+\mathrm{}\},\hfill \\ \hfill \delta _n=& C_n\alpha _s^n,\hfill \end{array}$$
$`(4.1\mathrm{a})`$
with the coefficients $`C_n`$ given in refs. 14 (see also ref. 11). Numerically one has

$$\begin{array}{cc}\hfill \delta _1(b)=& 0.090,\hfill \\ \hfill \delta _2(b)=& 0.045,\hfill \\ \hfill \delta _3(b)=& 0.029;\hfill \end{array}\begin{array}{cc}\hfill \delta _1(c)=& 0.137,\hfill \\ \hfill \delta _2(c)=& 0.108,\hfill \\ \hfill \delta _3(c)=& 0.125.\hfill \end{array}$$
$`(4.2)`$
Thus, the series relating pole and $`\overline{\text{MS}}`$ masses is at the edge of the region of convergence for the $`b`$ quark, and clearly diverges for the third order contribution for the $`c`$ quark.

We can use (4.2) and the formula for $`M(\mathrm{{\rm Y}})`$ in terms of the pole mass to express the former in terms of $`\overline{m}_b`$: One finds, for the perturbative contributions, neglecting $`m_c^2/m_b^2`$, and with $`\mathrm{\Lambda }`$ as before,

$$\begin{array}{cc}\hfill M(\mathrm{{\rm Y}})=& 2\overline{m}_b(\overline{m}_b)\{1+C_F\frac{\alpha _s(\overline{m})}{\pi }\hfill \\ \hfill +& 7.559\left(\frac{\alpha _s(\overline{m})}{\pi }\right)^2+43.502\left(\frac{\alpha _s}{\pi }\right)^3\}.\hfill \end{array}$$
$`(4.3)`$
(We have taken $`n_f=4`$). This does not look particularly convergent, and is certainly not an improvement over the perturbative expression using the pole mass where one has, for the choice $`\alpha _s(\mu =C_Fm_b\alpha _s)`$ and still neglecting the masses of quarks lighter than the $`b`$,

$$\begin{array}{cc}\hfill M(\mathrm{{\rm Y}})=& 2m_b\{12.193\left(\frac{\alpha _s}{\pi }\right)^224.725\left(\frac{\alpha _s}{\pi }\right)^3\hfill \\ \hfill & 458.28\left(\frac{\alpha _s}{\pi }\right)^4+897.93[\mathrm{log}\alpha _s]\left(\frac{\alpha _s}{\pi }\right)^5\}.\hfill \end{array}$$
$`(4.4)`$
To order three, (4.4) is actually better than (4.3); thus, and at least to third order, it is unclear that there is a connection between the rate of convergence and the use of pole or $`\overline{\text{MS}}`$ mass. What is more, at least the size of the known power corrections (namely, those stemming from the gluon condensate) do not favour the expression in terms of the $`\overline{\text{MS}}`$ mass. Indeed, if we evaluate the gluon condensate corrections to the direct formula (4.3) for $`M(\mathrm{{\rm Y}})`$ in terms of $`\overline{m}_b(\overline{m}_b)`$, it turns out that they are larger than than what one would have for the expression in terms of the pole mass, (4.4) ($`80`$ against $`9`$ MeV), because of the definition of the renormalization point. 5. DECAYS OF QUARKONIA Besides energies of bound states, once one has an effective interaction it is possible to evaluate, at least in principle, decay rates such as $`\mathrm{{\rm Y}}e^+e^{}`$. To relative order $`\alpha _s`$ (NLO) this has been done in ref. 5 (see also ref. 9) using the method of variations and the evaluation of ref. 25 for the hard part, $`\delta _{\mathrm{rad}}`$. We give here the exact result of a Rayleigh–Schrödinger perturbative calculation. One finds,

$$\begin{array}{cc}\hfill \mathrm{\Gamma }(\mathrm{{\rm Y}}& e^+e^{})=\mathrm{\Gamma }^{(0)}[1+\delta _{\mathrm{wf}}+\delta _{\mathrm{rad}}+\delta _{\mathrm{NP}}]^2,\hfill \\ \hfill \mathrm{\Gamma }^{(0)}=& 2\left[\frac{Q_b\alpha _{\mathrm{QED}}}{M(V)}\right]^2\left(mC_F\alpha _s(\mu ^2)\right)^3;\hfill \\ \hfill \delta _{\mathrm{wf}}+\delta _{\mathrm{rad}}=& \{\frac{3\beta _0}{4}(\mathrm{log}\frac{a\mu }{2}\gamma _\mathrm{E}\frac{\pi ^2}{9}+\frac{2}{3})\hfill \\ \hfill +& \frac{3c^{(1)}}{2}2C_F\}\frac{\alpha _s}{\pi };\hfill \end{array}$$
$`(5.1)`$
$`c^{(1)}`$ given in (2.1) and the nonperturbative contribution is

$$\delta _{\mathrm{NP}}=\frac{1}{2}\left[\frac{\mathrm{270\hspace{0.17em}459}}{\mathrm{108\hspace{0.17em}800}}+\frac{\mathrm{1\hspace{0.17em}838\hspace{0.17em}781}}{\mathrm{2\hspace{0.17em}890\hspace{0.17em}000}}\right]\frac{\pi \alpha _sG^2}{m^4\stackrel{~}{\alpha }_s^6}.$$

The corrections are very large. Because of this the calculation is unreliable, even for $`\overline{b}b`$, and fails completely for $`\overline{c}c`$. One might hope that this would be arranged by the $`O(\alpha _s^2)`$ corrections (NNLO); but this appears not to be the case. We have, first, “hard” corrections, similar to $`\delta _{\mathrm{rad}}`$; and corrections to the wave function, given in refs. 11,27 where we send for the (rather lengthy) explicit formulas.

Numerically, and with the width in keV,

$`\mathrm{\Gamma }(\mathrm{{\rm Y}}e^+e^{})`$ LO NLO NNLO                 $`\mu =m`$ $`0.41`$ $`1.22`$$`4_{4_4}^{4^4}`$ $`1.13`$   $`\mu =2/a`$ $`0.73`$ $`0.55`$$`4_{4_4}^{4^4}`$ $`0.80`$

The experimental figure is

$$\mathrm{\Gamma }(\mathrm{{\rm Y}}e^+e^{})=1.32\pm 0.04keV.$$

Clearly, the large NLO and NNLO perturbative corrections, both of similar size, and of the leading NP correction, make the theoretical result unstable, as the Table shows clearly.

For the (perhaps measurable) toponium (T) width we get, for $`m_t=175GeV`$ and still in keV,

$`\mathrm{\Gamma }(Te^+e^{})`$ LO NLO NNLO                 $`\mu =m_t`$ $`6.86`$ $`10.53`$$`4_{4_4}^{4^4}`$ $`13.0`$   $`\mu =2/a`$ $`10.24`$ $`10.91`$$`4_{4_4}^{4^4}`$ $`13.5`$

The situation has improved with respect to what we had for bottomium, but there is still a noticeable dependence on the renormalization point, and on the order of perturbation theory considered. We conclude on an estimate of some $`11\mathrm{\hspace{0.17em}14}`$ keV for the width.

ACKNOWLEDGEMENTS The author is grateful to CICYT, Spain, for partial financial support.

I would also like to acknowledge discussions with A. Pineda.

REFERENCES

1. F. J. Ynduráin, “Heavy Quarkonium”, Lectures at the XVII International School of Physics ”QCD: Perturbative or Nonperturbative”, Lisbon, 1999; FTUAM 99-32 (hep-ph/9910399).

2. W. Fischler, Nucl. Phys. B129 (1977) 157; A. Billoire, Phys. Lett. B92 (1980) 343, to one loop. M. Peter, Phys. Rev. Lett. 78 (1997) 602, with corrections in Y. Schröder, Phys. Lett. B447 (1999) 321, to two loops.

3. W. Buchmüller, Y. J. Ng and S.-H. H. Tye, Phys. Rev. D24 (1981) 3003.

4. S. N. Gupta and S. Radford, Phys. Rev. D24 (1981) 2309 and (E) D25 (1982) 3430; S. N. Gupta, S. F. Radford and W. W. Repko, ibid D26 (1982) 3305.

5. S. Titard and F. J. Ynduráin, Phys. Rev. D49 (1994) 6007; Phys. Rev. D51 (1995) 6348.

6. N. Brambilla et al., Phys. Lett. B470 (1999) 215.

7. W. E. Caswell and G. P. Lepage, Phys. Lett. B167 (1986) 437.

8. A. Pineda ad J. Soto, Nucl. Phys. B (Proc. Suppl.) 64 (1998) 428.

9. A. Pineda and F. J. Ynduráin, Phys. Rev. D58 (1998) 094022, and D61, 077505 (2000).

10. N. Brambilla et al., hep-ph/9903355, in press in Phys. Rev. The study of the singularities of the static QCD potential was first carried out by T. Appelquist, M. Dine and I. J. Muzinich, Phys. Rev. D17 (1978) 2074.

11. F. J. Ynduráin, hep-ph/0002237 (2000); FTUAM 00-14, hep-ph/0007333 (Proc. E. Schrödinger Institute Seminar on Confinement, May – June, 20000, to be published).

12. J. Santiago and F. J. Ynduráin, Nucl. Phys. B563 (1999) 45, and references therein

13. W. Lucha and F. F. Schöberl, UWThPh-1999-77 (hep-ph/0001191)

14. R. Coquereaux, Phys. Rev. D 23 (1981) 1365; R. Tarrach, Nucl. Phys. B183 (1981) 384; to two loops, N. Gray et al., Z. Phys. 48 (1990) 673. The $`O\left(\alpha _s^3\right)`$ corrections to the relations between the pole-$`\overline{\text{MS}}`$ masses have been recently evaluated by K. Melnikov and T. van Ritbergen, hep-ph/9912391.

15. M. Beneke and A. Signer, Phys. Lett. B471 (1999) 233.

16. N. Brambilla and A. Vairo, CERN-TH/20000-036 (hep-ph/0002075).

17. M. B. Voloshin, Nucl. Phys. B154 (1979) 365 and Sov. J. Nucl. Phys. 36 (1982) 143; H. Leutwyler, Phys. Lett. B98 (1981) 447.

18.  H. Dosch, Phys. Lett. B190 (1987) 177; Yu. A. Simonov, Nucl. Phys. B307 (1988) 512 and B324 (1989) 56; H. Dosch and Yu. A. Simonov, Phys. Lett. B205 (1988) 339; Yu. A. Simonov, S. Titard and F. J. Ynduráin, Phys. Lett. B354 (1995) 435.

19. U. Aglietti and Z. Ligeti, Phys. Lett. B364 (1995) 75.

20. R. Akhoury and V. I. Zakharov, Phys. Lett. B438 (1998) 165; F. J. Ynduráin, Nucl. Phys. B (Proc. Suppl.) 64 (1998) 433.

21. I. I. Bigi and N. G. Uraltsev, Phys. Lett. B321 (1994) 412; I. I. Bigi et al., Phys. Rev. D50 (1994) 2234.

22. A. Pineda, Ph. D. Thesis, Univ. Barcelona (1998), A. Hoang, M. Smith, T. Stelzer and S. Willenbrock, Phys. Rev. D59 (1999) 114014 and M. Beneke, Phys.Lett. B434 (1998) 115 for the constant term; and N. Brambilla et al., Phys. Rev. D 60 (1999) 091502 and Nucl. Phys. B566 (2000) 275 for the term in $`C_1`$ in (3.3).

23. Yu. A. Simonov, “QCD and Topics in Hadron Physics”, Lectures at the XVII Int’nal School of Physics ”QCD: Perturbative or Nonperturbative”, Lisbon, 1999, hep-ph/9911237.

24. J. L. Richardson, Phys. Lett. B82 (1979) 272. For an application to deep inelastic scattering, see K. Adel, F. Barreiro and F. J. Ynduráin, Nucl. Phys., B495 (1997) 221.

25. R. Barbieri et al., Phys. Lett. 57B (1975) 455; ibid. Nucl. Phys. B154 (1979) 535.

26. A. H. Hoang and T. Teubner, Phys. Rev. D58 (1998) 114023; A. Czarnecki and K. Melnikov, Phys. Rev. Lett. 80 (1998) 2531; M. Beneke, A. Singer and V. A. Smirnov, Phys. Rev. Lett. 80 (1998) 2535.

27. A. A. Penin and A. A. Pivovarov, Nucl. Phys. B549 (1999) 217; K. Melnikov and A. Yelkhovsky, Phys. Rev. D59 (1999) 114009.