text:
# 1 Introduction

## 1 Introduction

One of the main goals of future $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ colliders such as the Next Linear Collider (NLC) and the Future Muon Collider (FMC) will be to measure and to determine the properties of the top quark which was first discovered at the Tevatron  with a mass of $`m=174.3\pm 5\mathrm{GeV}`$ . Although the top quark will be studied at the LHC and the Tevatron (RUN-II) with an expected accuracy for the mass of $`23\mathrm{GeV}`$, the most accurate measurement of the mass with an accuracy of $`0.1\%`$ ($`100200\mathrm{MeV}`$) is expected to be obtained only at the NLC .

Due to the large top quark width, the top-antitop pair cannot hadronize into toponium resonances. The cross section appears therefore to have a smooth line-shape showing only a moderate $`1S`$ peak. In addition the top quark width serves as an infrared cutoff  and as a natural smearing over the energy . As a result, the nonperturbative QCD effects induced by the gluon condensate are small , allowing us to calculate the cross section with high accuracy by using perturbative QCD even in the threshold region. Many theoretical studies at LO and NLO have been done in the past for the total cross section , for the momentum distribution , also accounting for electro-weak corrections , and for the complete NLO correction including non-factorizable corrections .

Recently, the NNLO analysis of the inclusive threshold production cross section has been reported . The results of the NNLO analysis are summarized in a review article . To summarize the results for a standard approach using the pole mass, the NNLO corrections are uncomfortably large, spoiling the possibility for the top quark mass extraction at NLC with good accuracy because the $`1S`$ peak is shifted by about $`0.5\mathrm{GeV}`$ by the NNLO, the last known correction. One of the main reasons for this is the usage of the pole mass in the calculations. It was realized that such type of instability is caused by the fact that the pole mass is a badly defined object within full QCD .

In this paper we suggest a definition of the quark mass alternative to the pole mass. We call it $`\overline{\mathrm{PS}}`$ mass because in the static limit this mass coincides with the potential subtracted (PS) mass, even though the PS mass is defined in a different manner . This is the reason why we include both of them in this paper. The ratio of the small normalization point and the mass introduce relativistic corrections (called recoil corrections) which result in the difference between the two masses at higher orders in $`1/m`$. In contrast to the pole mass, the PS and the $`\overline{\mathrm{PS}}`$ mass are not sensitive to non-perturbative QCD effects. We derive recoil corrections to the relation of the pole mass to the $`\overline{\mathrm{PS}}`$ mass and demonstrate that perturbative predictions for the cross section become much more stable at higher orders of QCD (shifts are below $`0.1\mathrm{GeV}`$) if we use the PS or $`\overline{\mathrm{PS}}`$ mass for the calculations, as it is the case for any of the other threshold masses. This understanding removes one of the obstacles for the accurate top quark mass measurement and it can be expected that the top quark mass will be extracted from a threshold scan at NLC with an accuracy of about $`100200\mathrm{MeV}`$. The necessity of isolating the IR contributions in the mass calculation and the consideration of a short distance mass have been studied intensively . The applications of the PS, LS and $`1S`$ mass have been reported recently by several groups  and reviewed in Ref. . The detailed comparison of our results using the PS mass with results of other groups have been performed in Ref.  so that we refer the reader to this reference for details.

This paper is therefore devoted to two subjects. First, we define the $`\overline{\mathrm{PS}}`$ mass which allows us to calculate recoil corrections near the threshold. Second, we calculate the cross section for $`e^+e^{}t\overline{t}`$ near threshold with NNLO accuracy using three alternative mass concepts, i.e. (1) the pole mass, (2) the static PS mass, and (3) the new $`\overline{\mathrm{PS}}`$ mass. The results of the first approach has already been reported in Ref. , the preliminary results for the top quark pair production near threshold using the static PS mass were presented at the workshop “Physics at Linear Colliders”  by one of the authors. The results on the second and third alternative concepts which we present here are new.

The paper is organized as follows: In Sec. 2 we give the definition of the $`\overline{\mathrm{PS}}`$ mass in terms of the soft part of the self energy and calculate the leading order contribution. In Sec. 3 we calculate the two-loop corrections and present our final result. In Sec. 4 we use the obtained results for the analysis of the top quark pair production. In Sec. 5 we give our conclusions. The Appendix contains explicit expressions for the coefficients which occur in the text.

## 2 Mass definitions

The top quark mass is an input parameter of the Standard Model. Although it is widely accepted that the quark masses are generated due to the Higgs mechanism, the value of the mass cannot be calculated from the Standard Model. Instead, quark masses have to be determined from the comparison of theoretical predictions and experimental data.

It is important to stress that there is no unique definition of the quark mass. Because the quark cannot be observed as a free particle like the electron, the quark mass is a purely theoretical notion and depends on the concept adopted for its definition. The best known definitions are the pole mass and the $`\overline{\mathrm{MS}}`$ mass. However, both definitions are not adequate for the analysis of top quark production near threshold. The pole mass should not be used because it has the renormalon ambiguity and cannot be determined more accurately than $`300400\mathrm{MeV}`$  (see also Refs. ). The $`\overline{\mathrm{MS}}`$ mass is an Euclidean mass, defined at high virtuality, and therefore destroys the non-relativistic expansion. Instead, it was recently suggested to use threshold masses like the low scale (LS) mass , the potential subtracted (PS) mass , or one half of the perturbative mass of a fictious $`1^3S_1`$ ground state (called $`1S`$ mass) . In this paper we study the static PS mass suggested in Ref. ,

$`m_{\mathrm{PS}}=m_{\mathrm{pole}}\delta m_{\mathrm{PS}}\text{with}\delta m_{\mathrm{PS}}={\displaystyle \frac{1}{2}}{\displaystyle ^{|\stackrel{}{k}|<\mu _f}}{\displaystyle \frac{d^3k}{(2\pi )^3}}V_C(|\stackrel{}{k}|)`$ (1)
where $`V_C`$ is the quark-antiquark Coulomb potential. In order to understand why this mass definition is more adequate than the pole mass and to see that the pole mass is very sensitive to long distance effects, it is enough to consider the one-loop expression for the self energy diagram. Taking the residue in $`k_0`$, one obtains a soft self energy contribution which comes from momenta $`k`$ with $`|\stackrel{}{k}|<\mathrm{\Lambda }_{\mathrm{QCD}}`$,

$$\delta m=4\pi \alpha _sC_F\underset{|\stackrel{}{k}|<\mathrm{\Lambda }_{\mathrm{QCD}}}{}\frac{d^3k}{(2\pi )^3}\frac{1}{2|\stackrel{}{k}|^2}=\frac{\alpha _sC_F}{\pi }\mathrm{\Lambda }_{\mathrm{QCD}}.$$
(2)
We observe that the pole mass has a non-perturbative uncertainty of order $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ which then penetrates into consequent perturbative QCD calculations. It is easy to realize that the PS mass is free of this ambiguity. Indeed, the term $`\delta m`$ in Eq. (2) cancels in the definition of the PS mass as given in Eq. (1) as well as in the combination $`2m_{\mathrm{pole}}+V(r)`$. The definition in Eq. (1) has been given in Ref.  but has already been discussed implicitly in Ref. . The remarkable step made in Ref.  is to use this definition beyond one-loop order. It has been proven in Ref.  that the cancellation of the infrared QCD contributions to the PS mass in Eq. (1) holds even at higher loop orders.

In this paper we suggest a natural definition of a short distance quark mass near threshold. Our objective is that the definition should be gauge independent and well-defined within quantum field theory so that radiative and relativistic corrections can be calculated in a systematic way. Nevertheless, we stress that both concepts, as well as the other concepts mentioned before, are equivalently applicable to considerations near the quark production threshold.

### 2.1 The $`\overline{\mathrm{PS}}`$ mass

We define

$`m_{\overline{\mathrm{PS}}}=m_{\mathrm{pole}}\delta m_{\overline{\mathrm{PS}}}\text{with}\delta m_{\overline{\mathrm{PS}}}=\mathrm{\Sigma }_{\mathrm{soft}}(p/)|_{p/=m}`$ (3)
where $`\mathrm{\Sigma }_{\mathrm{soft}}`$ is the soft part of the heavy quark self energy which is defined as the part where at least one of the heavy quark propagators is on-shell. A more precise definition is given below. We will show that the static PS mass and the new $`\overline{\mathrm{PS}}`$ mass definition coincide at leading order in $`1/m`$. At the same time we stress that the new definition accounts for recoil corrections of orders $`1/m`$ and $`1/m^2`$. But first we describe the physical picture and discuss the structure of the quark self energy shown in Fig. 1.

The starting point of our considerations is an on-shell quark with mass $`m`$ and momentum $`p`$ (i.e. $`p^2=m^2`$) which we consider to be at rest, $`p=(m,\stackrel{}{0})`$. This quark interacts with a number of gluons, the subdiagram $`S`$ displayed in Fig. 1 describes the interaction between the gluons. In general the quark lines between the interaction points represent virtual quark states. However, if the virtual quark comes very close to the mass shell and the total momentum of the cloud of virtual gluons becomes soft, this situation gives rise to long-distance nonperturbative QCD interactions. The described (virtual) contributions result in the soft part of the self energy, $`\mathrm{\Sigma }_{\mathrm{soft}}`$. For a precise definition we start with a general self energy diagram as shown in Fig. 1,

$$i\mathrm{\Sigma }(p/)=\underset{m=1}{\overset{M}{}}\frac{d^4l_m}{(2\pi )^4}S_{\{\alpha _n\}}^{\{a_n\}}(\{l_m\})(ig_s\gamma ^{\alpha _{N+1}}T_{a_{N+1}})\underset{n=N}{\overset{1}{}}\frac{i}{p/_nm}(ig_s\gamma ^{\alpha _n}T_{a_n})$$
(4)
where the last factor is a non-commutative product with decreasing index $`n`$. The line momenta $`k_n`$ are linear combinations of the gluon loop momenta $`l_m`$, the particular representation is specified by the structure $`S`$. The symbol $`\{l_m\}`$ means the set of all these loop momenta, the same symbol is used for the Lorentz and color indices. In general we have $`MN`$ which means that line momenta can be correlated. The momenta of the virtual quark states are given by $`p_n=p+k_n`$. Taking this as the starting point we define

$`i\mathrm{\Sigma }_{\mathrm{soft}}(p/)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle }{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{d^4l_m}{(2\pi )^4}}S_{\{\alpha _n\}}^{\{a_n\}}(\{l_m\})(ig_s\gamma ^{\alpha _{N+1}}T_{a_{N+1}}){\displaystyle \underset{n=N}{\overset{i+1}{}}}{\displaystyle \frac{i}{p/_nm}}(ig_s\gamma ^{\alpha _n}T_{a_n})\times `$ (5)
$`\times i(p/_i+m)(i\pi \delta (p_i^2m^2))(ig_s\gamma ^{\alpha _i}T_{a_i}){\displaystyle \underset{n=i1}{\overset{1}{}}}{\displaystyle \frac{i}{p/_nm}}(ig_s\gamma ^{\alpha _n}T_{a_n}).`$
This equation is the definition of the soft part of the quark self energy. One can derive this expression from Eq. (4) by using the identity

$$\frac{1}{p^2m^2+iϵ}=i\pi \delta (p^2m^2)+P\left(\frac{1}{p^2m^2}\right)$$
(6)
and the fact that the principal value integral does not give any infrared sensitive contribution. The delta function can be used to remove the integration over the zero component of $`k_i`$. In order to parameterize the softness of the gluon cloud we impose a cutoff on the spatial component, $`|\stackrel{}{k}_i|<\mu _f`$, and indicate this by a label $`\mu _f`$ written at the upper limit of the three-dimensional integral. This cutoff $`\mu _f`$ is also known as factorization scale. So we can rewrite Eq. (5) as

$$\mathrm{\Sigma }_{\mathrm{soft}}(p/,\mu _f)=\frac{1}{2}\underset{i=1}{\overset{N}{}}^{\mu _f}\frac{d^3k_i}{(2\pi )^3}V(\stackrel{}{k}_i,p)$$
(7)
where

$`V(\stackrel{}{k}_i,p)`$ $`:=`$ $`{\displaystyle }{\displaystyle \underset{m=1}{\overset{M1}{}}}{\displaystyle \frac{d^4l_m}{(2\pi )^4}}S_{\{\alpha _n\}}^{\{a_n\}}(\{l_m\})(ig_s\gamma ^{\alpha _{N+1}}T_{a_{N+1}}){\displaystyle \underset{n=N}{\overset{i+1}{}}}{\displaystyle \frac{i}{p/_nm}}(ig_s\gamma ^{\alpha _n}T_{a_n})\times `$ (8)
$`\times {\displaystyle \frac{p/_i+m}{2p_i^0}}\left(ig_s\gamma ^{\alpha _i}T_{a_i}\right){\displaystyle \underset{n=i1}{\overset{1}{}}}{\displaystyle \frac{i}{p/_nm}}\left(ig_s\gamma ^{\alpha _n}T_{a_n}\right).`$
The range of the index $`m`$ is reduced by one which indicates that one of the loop momenta is extracted as line momentum of the $`i`$-th line. In the following we deal with the different realizations of this compact expression. As we will see explicitly, the function $`V(\stackrel{}{k},p)`$ occurring as integrand can be seen as quark-antiquark potential where we have summed over the spin of the tensor product of a final state with an initial state. Because the static quark-antiquark potential is used in a similar way in Ref. , we recover the result of Ref.  in the static limit. However, there is no kind of hierarchical order between both concepts because both solve the afore mentioned problems with non-perturbative uncertainties of the order $`O(\mathrm{\Lambda }_{\mathrm{QCD}})`$. The role of non-perturbative effects of the order $`1/m`$ have to be studied elsewhere.

### 2.2 The leading order perturbative contribution

The leading order contribution to the self energy of the quark is given by

$$\mathrm{\Sigma }(p/)=i\frac{d^4k}{(2\pi )^4}(ig_s\gamma _\alpha T_a)\frac{i}{p/+/km}(ig_s\gamma ^\alpha T_a)\frac{i}{k^2}$$
(9)
where Feynman gauge is used for the gluon. The soft contribution thus reads

$`\mathrm{\Sigma }_{\mathrm{soft}}(p/)`$ $`=`$ $`ig_s^2C_F{\displaystyle }{\displaystyle \frac{d^4k}{(2\pi )^4k^2}}\gamma _\alpha (p/+/k+m)\gamma ^\alpha (i\pi \delta ((p+k)^2m^2))=`$ (10)
$`=`$ $`\pi g_s^2C_F{\displaystyle }{\displaystyle \frac{d^4k}{(2\pi )^4k^2}}(2(p/+/k)+4m)\delta ((p+k)^2m^2).`$
The projector $`(1+\gamma ^0)/2`$ which represents the on-shell quark at rest can be placed at the left and at the right of the Dirac structure and will lead to a further simplification. In terms of the zero component of the momentum $`k`$ the Dirac delta function has two zeros $`k_0=k_+`$ and $`k_0=k_{}`$ with

$$k_\pm :=\pm \sqrt{\kappa ^2+m^2}m$$
(11)
where $`\kappa =|\stackrel{}{k}|`$. The delta function is therefore written as

$$\delta \left((p+k)^2m^2\right)=\frac{1}{2\sqrt{\kappa ^2+m^2}}\left(\delta (k_0k_+)+\delta (k_0k_{})\right).$$
(12)
The procedure which is done here is shown in Fig. 2(a–b). The cross indicates that we cut the line at this point by imposing the on-shell condition to the corresponding (virtual) momentum. The diagram then proceeds to a quark-antiquark interaction diagram where we have kept the crosses to indicate the position of the cut line. This line carries the momentum $`p+k`$ while the other two external lines carry the momentum $`p`$. Accordingly we obtain

$$\mathrm{\Sigma }_{\mathrm{soft}}(p/,\mu _f)=\frac{1}{2}^{\mu _f}\frac{d^3k}{(2\pi )^3}V(\stackrel{}{k},p),V(\stackrel{}{k},p)=V_+(\stackrel{}{k},p)+V_{}(\stackrel{}{k},p)$$
(13)
where

$$V_\pm (\stackrel{}{k},p)=g_s^2C_F\frac{(m+k_\pm )+2m}{\sqrt{m^2+\kappa ^2}(k_\pm ^2k^2)}=\frac{g_s^2C_F(\sqrt{m^2+\kappa ^2}2m)}{2m\sqrt{m^2+\kappa ^2}(\sqrt{m^2+\kappa ^2}m)}.$$
(14)
For $`m\mu _f`$ the restriction of the three-dimensional integral allows for an expansion in $`\kappa /m`$. We obtain

$`V_+(\stackrel{}{k},p)`$ $`=`$ $`4\pi \alpha _sC_F\left\{{\displaystyle \frac{1}{\kappa ^2}}{\displaystyle \frac{3}{4m^2}}+O\left({\displaystyle \frac{\kappa ^2}{m^4}}\right)\right\},`$
$`V_{}(\stackrel{}{k},p)`$ $`=`$ $`4\pi \alpha _sC_F\left\{{\displaystyle \frac{3}{4m^2}}+O\left({\displaystyle \frac{\kappa ^2}{m^4}}\right)\right\}`$ (15)
and therefore

$$V(\stackrel{}{k},p)=4\pi \alpha _sC_F\left\{\frac{1}{\kappa ^2}\frac{3}{2m^2}+O\left(\frac{\kappa ^2}{m^4}\right)\right\}.$$
(16)
The first term is the Coulomb potential for a quark-antiquark interaction. The second term can be related to the Breit-Fermi potential of the quark-antiquark interaction  by summing over the spin states of the tensor product of the final quark and the final antiquark state and using the same kinematic constraints. Moreover, one can identify $`V_+(\stackrel{}{k},p)`$ with the scattering potential and $`V_{}(\stackrel{}{k},p)`$ with the annihilation potential.

One can integrate this radial symmetric potential over the space components and end up with an expansion in $`1/m`$. Note that there are no contributions to odd powers of $`1/m`$. But one can also integrate the exact expressions for $`V_\pm (\stackrel{}{k},p)`$ in Eq. (14) up to the factorization scale $`\mu _f`$ and obtain

$$\mathrm{\Sigma }_{\mathrm{soft}}(\mu _f)=\frac{\alpha _sC_F}{2\pi }m\left\{3\mathrm{ln}\left(\frac{\mu _f}{m}+\sqrt{\frac{\mu _f^2}{m^2}+1}\right)\frac{\mu _f}{m}\sqrt{\frac{\mu _f^2}{m^2}+1}\right\}.$$
(17)
The expansion of this expression in small values of $`\mu _f/m`$ results in

$$\mathrm{\Sigma }_{\mathrm{soft}}(\mu _f)=\frac{\alpha _sC_F}{\pi }\mu _f\left\{1\frac{\mu _f^2}{2m^2}\right\}.$$
(18)
The first term reproduces the result given in Ref.  to leading order in $`\alpha _s`$ while the second term is the recoil correction to the static limit in this order of perturbation theory. This second term is related to the Breit-Fermi potential but does not coincide with it.

## 3 Two loop contributions

To take a step beyond the leading order perturbation theory, we consider two-loop diagrams for the heavy quark self energy as shown in Fig. 3.

We calculate them in Coulomb gauge, even though we stress that our final result is gauge invariant. The gluon propagator in Coulomb gauge is given by

$$G_{00}^{ab}(k)=\frac{i\delta ^{ab}}{\stackrel{}{k}^2},G_{ij}^{ab}(k)=\frac{i\delta ^{ab}}{k^2}\left(\delta _{ij}\frac{k_ik_j}{\stackrel{}{k}^2}\right),i,j=1,2,3.$$
(19)
The use of Coulomb gauge splits up the gluon propagators into a Coulomb term (Coulomb gluon) and a transverse term (transverse gluon) where the first one couples to the quark via the time components only. This splitting is shown in Fig. 2(b–d).

### 3.1 The abelian diagrams

We start our consideration with the abelian diagrams shown in Figs. 3(a) and (b). In cutting the quark line in all possible ways we obtain a lot of diagrams. However, we find that the final contribution of these diagrams to the soft part of the self energy is suppressed by $`\mu _f^2/m^2`$. There can be found different arguments for this suppression. First, in applying the classical Ward identity, the QED diagrams shown in Fig. 4 cancel exactly at $`|\stackrel{}{k}|0`$, the remaining contribution is of order $`O(\stackrel{}{k}^2/m^2)`$. Concerning this note that the Ward identity for the interaction vertex of a Coulomb gluon with the quark holds even in non-abelian theories .

A second argument is that the interaction between a transverse gluon and a non-relativistic quark as shown in Fig. 5(a) is suppressed by $`\mu _f/m`$, leading to an overall $`\mu _f^2/m^2`$ suppression. In addition, the box diagrams in Figs. 5(b) and (c) are either suppressed by a factor $`\mu _f^2/m^2`$ or give an iteration of the leading order potential. To summarize, the diagrams in Figs. 3(a) and (b) give contributions only of the order $`g_s^4\mu _f^2/m^2`$.

### 3.2 The vacuum polarization of the gluon

Following these arguments, it turns out that the only abelian diagrams which can give a non-suppressed contribution to the soft part of the quark self energy are the diagrams containing the vacuum polarization of the gluon as shown in Fig. 3(d–f). The simple calculation of these diagrams within the $`\overline{\mathrm{MS}}`$ scheme, accounting only for light fermion loops, gluon loop (and ghost loop if Feynman gauge is used) results after renormalization in

$`\mathrm{\Sigma }_{\mathrm{soft}}^A`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle ^\mu }{\displaystyle \frac{d^3k}{(2\pi )^3}}({\displaystyle \frac{4\pi \alpha _s(\mu )C_F}{|\stackrel{}{k}|^2}})\times `$ (20)
$`\times \left\{1+{\displaystyle \frac{\alpha _s(\mu )}{4\pi }}\left({\displaystyle \frac{31C_A}{9}}{\displaystyle \frac{20T_FN_F}{9}}\left({\displaystyle \frac{11C_A}{3}}{\displaystyle \frac{4T_FN_F}{3}}\right)\mathrm{ln}\left({\displaystyle \frac{|\stackrel{}{k}|^2}{\mu ^2}}\right)\right)\right\}`$
$`=`$ $`{\displaystyle \frac{\alpha _S(\mu )C_F}{\pi }}\mu _f\left\{1+{\displaystyle \frac{\alpha _s(\mu )}{4\pi }}\left(a_1\beta _0\mathrm{ln}\left({\displaystyle \frac{\mu _f^2}{\mu ^2}}\right)\right)\right\}.`$
This result has been anticipated because the expression in the curly brackets of the integrand reproduces the next-to-leading order correction to the QCD Coulomb potential.

### 3.3 The non-abelian diagram

In this subsection we calculate the non-abelian diagram shown in Fig. 3(c). In Coulomb gauge this diagram gives rise to seven two-loop diagrams which are shown in Fig. 6.

Direct calculations show that only the diagram in Fig. 6(b) gives a contribution of order $`g_s^4\mu _f/m`$ while the other diagrams are of order $`g_s^4\mu _f^2/m^2`$ or vanish to this order in the coupling after applying the renormalization procedure (see e.g. Ref. ). The calculation of the diagram in Fig. 6(b) is simple and we show it in detail. The contribution of this diagram to the self energy is given by

$`i\mathrm{\Sigma }^{6b}(p/)`$ $`=`$ $`{\displaystyle }{\displaystyle \frac{d^4k_1}{(2\pi )^4}}{\displaystyle \frac{d^4k_2}{(2\pi )^4}}{\displaystyle \frac{i}{\stackrel{}{k}_1^2}}{\displaystyle \frac{i}{\stackrel{}{k}_2^2}}{\displaystyle \frac{i}{(k_1k_2)^2}}(\delta _{ij}{\displaystyle \frac{(k_1k_2)_i(k_1k_2)_j}{(\stackrel{}{k}_1\stackrel{}{k}_2)^2}})g_sf_{abc}(k_1+k_2)^jg^{00}\times `$ (21)
$`\times {\displaystyle \frac{1+\gamma ^0}{2}}(ig_s\gamma ^0T_a){\displaystyle \frac{i}{p/+/k_2m}}(ig_s\gamma ^iT_b){\displaystyle \frac{i}{p/+/k_1m}}(ig_s\gamma ^0T_c){\displaystyle \frac{1+\gamma ^0}{2}}.`$
Here the Lorentz structure of the three-gluon vertex is reduced to $`k_1+k_2`$. The Dirac structure of the integrand can be simplified to

$$\frac{1+\gamma ^0}{2}\left\{(2m+k_{20})\gamma ^i\stackrel{}{k}_1\stackrel{}{\gamma }\stackrel{}{k}_2\stackrel{}{\gamma }\gamma ^i(2m+k_{10})\right\}.$$
(22)
The rest of the diagram is symmetric with respect to the interchange of the two line momenta $`k_1`$ and $`k_2`$. That provides us with a further simplification of the Dirac structure, the self energy reduces to

$$\mathrm{\Sigma }^{6b}(p/)=g_s^4C_FC_A\frac{d^4k_1}{(2\pi )^4}\frac{d^4k_2}{(2\pi )^4}\frac{(4m+k_{10}+k_{20})\left(\stackrel{}{k}_1^2\stackrel{}{k}_2^2(\stackrel{}{k}_1\stackrel{}{k}_2)^2\right)}{\stackrel{}{k}_1^2\stackrel{}{k}_2^2(\stackrel{}{k}_1\stackrel{}{k}_2)^2(k_1k_2)^2((p+k_1)^2m^2)((p+k_2)^2m^2)}$$
(23)
where $`T_aT_bT_cf_{abc}=iC_FC_A/2`$ has been used. We now employ the replacements

$$\frac{1}{(p+k_1)^2m^2}i\pi \delta ((p+k_1)^2m^2),\frac{1}{(p+k_2)^2m^2}i\pi \delta ((p+k_2)^2m^2)$$
(24)
i.e. we set cuts at the two intermediate quark lines separately to obtain the two parts $`\mathrm{\Sigma }_{\mathrm{soft1}}^{6b}`$ and $`\mathrm{\Sigma }_{\mathrm{soft2}}^{6b}`$ of the soft contribution $`\mathrm{\Sigma }_{\mathrm{soft}}^{6b}`$ of the self energy, as shown in Fig. 7.

Taking the first cut, the delta function removes the integration over $`k_{10}`$. At the same time our definition of the soft contribution imposes a restriction $`|\stackrel{}{k}_1|<\mu _f`$ on the space components of the first line momentum. The delta function

$$\delta \left((p+k_1)^2m^2\right)=\frac{1}{2\sqrt{m^2+\stackrel{}{k}_1^2}}\left(\delta (k_{10}k_{1+})+\delta (k_{10}k_1)\right)$$
(25)
with

$$k_{1\pm }=\pm \sqrt{m^2+\stackrel{}{k}_1^2}m$$
(26)
causes two contributions which are known as scattering and annihilation amplitude (according to $`k_{1+}`$ and $`k_1`$, resp.). The integration over $`k_{20}`$ is done by using the residue theorem. Actually there are only two denominator factors which can contribute to poles of the integrand, namely $`(k_1k_2)^2`$ and $`((p+k_2)^2m^2)`$. In summing up the four contributions from the integration over $`k_{20}`$ and $`k_{10}`$ we end up with a two-fold three-dimensional integral over the space components of the two line momenta where the first integral is restricted by $`|\stackrel{}{k}_1|<\mu _f`$ as mentioned earlier. We now impose the restriction $`|\stackrel{}{k}_1|<\mu _f`$ on the integrand to simplify it and obtain

$$\mathrm{\Sigma }_{\mathrm{soft1}}^{6b}=\frac{g_s^4C_FC_A}{4m}^{\mu _f}\frac{d^3k_1}{(2\pi )^3}\frac{d^3k_2}{(2\pi )^3}\frac{\stackrel{}{k}_1^2\stackrel{}{k}_2^2(\stackrel{}{k}_1\stackrel{}{k}_2)^2}{\stackrel{}{k}_1^2(\stackrel{}{k}_2^2)^2(\stackrel{}{k}_1\stackrel{}{k}_2)^2}.$$
(27)
The integral over the space components of $`k_2`$ can be easily done by executing the angular integration followed by the radial integration. We obtain

$$\frac{d^3k_2}{(2\pi )^3}\frac{\stackrel{}{k}_1^2\stackrel{}{k}_2^2(\stackrel{}{k}_1\stackrel{}{k}_2)^2}{\stackrel{}{k}_1^2(\stackrel{}{k}_2^2)^2(\stackrel{}{k}_1\stackrel{}{k}_2)^2}=\frac{1}{16|\stackrel{}{k}_1|}$$
(28)
and therefore, finally

$$\mathrm{\Sigma }_{\mathrm{soft1}}^{6b}=\frac{\alpha _s^2C_FC_A}{16m}\mu _f^2.$$
(29)
Symmetry considerations show that $`\mathrm{\Sigma }_{\mathrm{soft2}}^{6b}`$ gives exactly the same contribution. As mentioned before, there are no other non-abelian contributions, therefore we obtain

$$\mathrm{\Sigma }_{\mathrm{soft}}^{NA}=\frac{\alpha _s^2C_FC_A}{8m}\mu _f^2.$$
(30)
This result has been anticipated, too, to be minus one half of the non-abelian correction to the QCD Coulomb potential, which is known in the literature (see for example Refs. ),

$`\mathrm{\Sigma }_{\mathrm{soft}}^{NA}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle ^{\mu _f}}{\displaystyle \frac{d^3k}{(2\pi )^3}}\left\{{\displaystyle \frac{\pi ^2\alpha _s^2C_FC_A}{m|\stackrel{}{k}|}}\right\}={\displaystyle \frac{\alpha _s^2C_FC_A}{8m}}\mu _f^2.`$ (31)
This calculation concludes the considerations of the two-loop diagrams shown in Fig. 3.

### 3.4 Our final result

Summarizing all contribution up to NNLO accuracy, we obtain

$$m_{\overline{\mathrm{PS}}}(\mu _f)m=\frac{\alpha _s(\mu )C_F}{\pi }\mu _f\left\{1+C_0^{}\frac{\mu _f}{m}+C_0^{\prime \prime }\frac{\mu _f^2}{m^2}+\frac{\alpha _s(\mu )}{4\pi }\left(C_1+C_1^{}\frac{\mu _f}{m}\right)+C_2\left(\frac{\alpha _s(\mu )}{4\pi }\right)^2\right\}$$
(32)
where $`m`$ is the pole mass, $`\mu `$ is the renormalization scale. This result is considered to be of the order $`O(\alpha _s^2)`$ because we imply that the ratio $`\mu _f/m`$ is typically of the order of $`\alpha _s`$ or smaller. The scale $`\mu _f`$ is the factorization scale, and

$`C_0`$ $`=`$ $`1,C_0^{}=0,C_0^{\prime \prime }={\displaystyle \frac{1}{2}},`$
$`C_1`$ $`=`$ $`a_12\beta _0\mathrm{ln}\left({\displaystyle \frac{\mu _f}{\mu }}\right),C_1^{}=C_A{\displaystyle \frac{\pi ^2}{2}},`$
$`C_2`$ $`=`$ $`a_22(2a_1\beta _0+\beta _1)\left(\mathrm{ln}\left({\displaystyle \frac{\mu _f}{\mu }}\right)1\right)+4\beta _0^2\left(\mathrm{ln}^2\left({\displaystyle \frac{\mu _f}{\mu }}\right)2\mathrm{ln}\left({\displaystyle \frac{\mu _f}{\mu }}\right)+2\right).`$ (33)
The constants $`a_1`$, $`a_2`$, $`\beta _0`$, and $`\beta _1`$ are given in Appendix A. The coefficients $`C_1`$ and $`C_2`$ have been derived in Ref.  by using known corrections to the QCD potential. In this work we have derived the coefficients $`C_0^{}`$, $`C_0^{\prime \prime }`$, and $`C_1^{}`$. Note that our result can be represented in a condensed form as

$$m_{\overline{\mathrm{PS}}}(\mu _f)m=\frac{1}{2}^{\mu _f}\frac{d^3k}{(2\pi )^3}\left(V_C(|\stackrel{}{k}|)+V_R(|\stackrel{}{k}|)+V_{NA}(|\stackrel{}{k}|)\right)$$
(34)
where the first term $`V_C`$ is the static Coulomb potential, $`V_R`$ is the relativistic correction (which is related to Breit-Fermi potential but does not coincide with it), and $`V_{NA}`$ is the non-abelian correction.

## 4 Top quark production near threshold <br>with NNLO accuracy

In this section we consider the cross section of the process $`e^+e^{}t\overline{t}`$ in the near threshold region where the velocity $`v`$ of the top quark is small. It is well-known that the conventional perturbative expansion does not work in the non-relativistic region because of the presence of the Coulomb singularities at small velocities $`v0`$. The terms proportional to $`(\alpha _s/v)^n`$ appear due to the instantaneous Coulomb interaction between the top and the antitop quark. The standard technique for resumming the terms $`(\alpha _s/v)^n`$ consists in using the Schrödinger equation for the Coulomb potential and to find the Green function . The Green function is then related to the cross section by the optical theorem. In order to determine NLO corrections to the cross section we need to know the short-distance correction to the vector current , the NLO correction to the Coulomb potential , and the contribution of the non-factorizable corrections . It was proven that the non-factorizable corrections cancel in the inclusive cross section at NLO and beyond .

At NNLO the situation is more complicated. One of the obstacles for a straightforward calculation are the UV divergences coming from relativistic corrections to the Coulomb potential. This problem can be solved by a proper factorization of the amplitudes and by employing effective theories. We want to sketch the derivation of the inclusive cross section. The inclusive cross section can be obtained from the correlation function of two vector currents $`j_\mu (x)=\overline{t}(x)\gamma _\mu t(x)`$,

$$\mathrm{\Pi }_{\mu \nu }(p^2)=id^4xe^{ipx}0|T\{j_\mu (x),j_\nu (0)\}|0.$$
(35)
The first step is to expand the Lagrangian and the currents

$$j_i=\overline{t}\gamma _it=c_1\psi ^{}\sigma _i\chi \frac{c_2}{6m^2}\psi ^{}\sigma _i(i\stackrel{}{D})^2\chi ,$$
(36)
consistently in $`1/m`$. The useful language for treating the one and two non-relativistic quark(s) is provided by the NRQCD  and the PNRQCD , respectively. After the expansion, the relative cross section reads

$$R=\sigma (e^+e^{}t\overline{t})/\sigma _{pt}=e_Q^2N_c\frac{24\pi }{s}C(r_0)\mathrm{Im}\left[\left(1\frac{\stackrel{}{p}^2}{3m^2}\right)G(r_0,r_0|E+i\mathrm{\Gamma })\right]|_{r_00}$$
(37)
where $`\sigma _{pt}=4\pi \alpha ^2/3s`$, $`e_Q`$ is the electric charge of the top quark, $`N_c`$ is the number of colors, $`\sqrt{s}=2m+E`$ is the total center-of-mass energy of the quark-antiquark system, $`m`$ is the top quark pole mass and $`\mathrm{\Gamma }`$ is the top quark width. The unknown coefficient $`C(r_0)`$ can be fixed by using a direct QCD calculation of the vector vertex at NNLO in the so-called intermediate region  and by using the direct matching procedure suggested in Ref. . The result of such a matching procedure for the coefficient $`C(r_0)`$ is given in the Appendix .

The function $`G(\stackrel{}{r},\stackrel{}{r}^{}|E+i\mathrm{\Gamma })`$ is the non-relativistic Green function. It satisfies the Schrödinger equation

$$(HEi\mathrm{\Gamma })G(\stackrel{}{r},\stackrel{}{r}^{}|E+i\mathrm{\Gamma })=\delta (\stackrel{}{r}\stackrel{}{r}^{}).$$
(38)
The non-relativistic Hamiltonian of the heavy quark-antiquark system reads

$$H=H_0+W(r),H_0=\frac{\stackrel{}{p}^2}{m}+V_C(r),$$
(39)
where $`V_C(r)`$ is the static QCD potential of the heavy quark-antiquark system at NNLO order,

$`V_C(r)`$ $`=`$ $`{\displaystyle \frac{\alpha _s(\mu )C_F}{r}}\{1+{\displaystyle \frac{\alpha _s(\mu )C_F}{4\pi }}(2\beta _0\mathrm{ln}(\mu ^{}r)+a_1)`$ (40)
$`+\left({\displaystyle \frac{\alpha _s(\mu )}{4\pi }}\right)^2(\beta _0^2(4\mathrm{ln}^2(\mu ^{}r)+{\displaystyle \frac{\pi ^2}{3}})+2(\beta _1+2\beta _0a_1)\mathrm{ln}(\mu ^{}r)+a_2)\},`$
$`\mu ^{}=\mu e^{\gamma _E}`$, $`\mu `$ is the renormalization scale, and $`\gamma _E`$ is Euler’s constant. The color factors are given by $`C_F=4/3`$, $`C_A=3`$ and $`T_F=1/2`$. The number of light quark flavors is $`N_F=5`$. The coefficients $`a_1`$ and $`a_2`$ were calculated in  and , respectively, and are listed in the Appendix.

The function $`W(r)`$ in Eq. (39) is the QCD generalization of the QED Breit-Fermi Hamiltonian . We consider here the quark-antiquark production in the $`S`$-wave mode. The Breit-Fermi Hamiltonian for the final state with $`\stackrel{}{L}=0`$, $`\stackrel{}{S}^2=2`$ reads

$$W(r)=\frac{\stackrel{}{p}^4}{4m^3}+\frac{11\pi C_F\alpha _s(\mu )}{3m^2}\delta (\stackrel{}{r})\frac{C_F\alpha _s(\mu )}{2m^2}\{\frac{1}{r},\stackrel{}{p}^2\}\frac{C_AC_F\alpha _s^2(\mu )}{2mr^2}.$$
(41)
It has been demonstrated in Refs.  that the Schrödinger equation can be reduced to the equation

$$(H_1E_1)G_1(r,r^{}|E_1)=\delta ^3(rr^{})$$
(42)
with the energy $`E_1=\overline{E}+\overline{E}^2/4m`$, $`\overline{E}=E+i\mathrm{\Gamma }`$, and with the Hamiltonian

$$H_1=\frac{\stackrel{}{p}^2}{m}+V_C(r)+\frac{3\overline{E}}{2m}V_0(r)\left(\frac{2}{3}+\frac{C_A}{C_F}\right)\frac{V_0^2(r)}{2m},V_0(r)=\frac{\alpha _s(\mu )C_F}{r}.$$
(43)
The final expression for the NNLO cross section is given by

$$R^{\mathrm{NNLO}}(E)=\frac{8\pi }{m^2}\left\{14C_F\frac{\alpha _s(m)}{\pi }+C_2(r_0)C_F\left(\frac{\alpha _s(m)}{\pi }\right)^2\right\}\mathrm{Im}\left[\left(1\frac{5\overline{E}}{6m}\right)G_1(r_0,r_0|E_1)\right]$$
(44)
with $`G_1(r_0,r_0|E_1)`$ being the solution of Eq. (42) at $`r=r^{}=r_0`$ while $`C_2(r_0)`$ is taken from Eqs. (AppendixAppendix). For the numerical solution we use the program derived in Ref.  by one of the authors.

### 4.1 Relations between the masses

Our main result is Eq. (32), it enables us to relate the pole mass to the $`\overline{\mathrm{PS}}`$ mass or the PS mass (the latter by leaving out the terms of higher order in $`1/m`$). Note that we did not include electroweak corrections neither in this mass relation nor in the cross section. To relate the pole mass to the $`\overline{\mathrm{MS}}`$ mass we use the three-loop relation which has been derived in Ref. ,

$`{\displaystyle \frac{m_{\mathrm{pole}}}{\overline{m}(\overline{m})}}`$ $`=`$ $`1+{\displaystyle \frac{4}{3}}\left({\displaystyle \frac{\alpha _s}{\pi }}\right)+\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^2(1.0414N_F+13.4434)`$ (45)
$`+\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^3(0.6527N_F^226.655N_F+190.595)`$
where $`\alpha _s=\alpha _s(\overline{m})`$ is taken at the $`\overline{\mathrm{MS}}`$ mass. More precisely, we fix the $`\overline{\mathrm{MS}}`$ mass to take the values $`\overline{m}(\overline{m})=160\mathrm{GeV}`$, $`165\mathrm{GeV}`$, and $`170\mathrm{GeV}`$ and use Eq. (45) to determine the pole mass at LO, NLO, and NNLO. This pole mass is then used as input parameter $`m`$ in Eq. (32) to determine the PS and $`\overline{\mathrm{PS}}`$ masses at LO, NLO, and NNLO. The obtained values for the PS and the $`\overline{\mathrm{PS}}`$ mass differ only in NNLO. The results of these calculations are collected in Table 1, together with the estimates for “large $`\beta _0`$” accuracy . Note that the same values for the $`\overline{\mathrm{MS}}`$ mass have been used for Tables 2 and 3 in Ref. . The obtained mass values can now be used for an analysis of the Schrödinger equation.

Taking the $`\overline{\mathrm{PS}}`$ mass instead of the PS mass, we observe an improvement of the convergence. The differences for the mass values in going from LO to NLO to NNLO to the “large $`\beta _0`$” estimate for e.g. $`\overline{m}(\overline{m})=165\mathrm{GeV}`$ read $`7.64\mathrm{GeV}`$, $`1.64\mathrm{GeV}`$, $`0.52\mathrm{GeV}`$, and $`0.25\mathrm{GeV}`$ for the pole mass, $`6.72\mathrm{GeV}`$, $`1.21\mathrm{GeV}`$, $`0.29\mathrm{GeV}`$, and $`0.08\mathrm{GeV}`$ for the PS mass and finally $`6.72\mathrm{GeV}`$, $`1.21\mathrm{GeV}`$, $`0.27\mathrm{GeV}`$, and $`0.08\mathrm{GeV}`$ for the $`\overline{\mathrm{PS}}`$ mass.

In Fig. 8 we show the difference between the $`\overline{\mathrm{PS}}`$ and the PS mass (in GeV) as a function of the factorization scale $`\mu _f`$ (solid line) at $`\mu =15\mathrm{GeV}`$. It is interesting to observe that the non-abelian part of the difference between the $`\overline{\mathrm{PS}}`$ and the PS mass (dotted line in Fig. 8) can be as large as $`200\mathrm{MeV}`$. But the recoil correction cancel in part the non-abelian one and therefore the final difference is not more than $`50\mathrm{MeV}`$. The dependence of $`m_{\overline{\mathrm{PS}}}m_{\mathrm{PS}}`$ on the renormalization scale at $`\mu _f=30\mathrm{GeV}`$ is given in Fig. 8 by the dashed curve.

### 4.2 The concept with the pole mass

The top quark cross section at LO, NLO, and NNLO is shown in Fig. 9 as a function of the center-of-mass energy. For the top quark pole mass we use $`m_t=175.05\mathrm{GeV}`$, for the top quark width $`\mathrm{\Gamma }_t=1.43\mathrm{GeV}`$, and for the QCD coupling constant $`\alpha _s(m_Z)=0.119`$ . Different values $`\mu =15\mathrm{GeV}`$, $`30\mathrm{GeV}`$, and $`60\mathrm{GeV}`$ for the renormalization scale are selected because they roughly correspond to the typical spatial momenta for the top quark. For solving the Schrödinger equation we use the program written by one of the authors . Note that we do not take into account an initial photon radiation which would change the shape of the cross section. This can be easily included in the Monte Carlo simulation.

The NNLO curve modifies the line shape by the amount of $`2030\%`$ which is as large as the NLO correction. It also shifts the position of the $`1S`$ peak by approximately $`600\mathrm{MeV}`$. These large shifts of the peak position were expected. As we discussed above (and is well-known in the literature ), the pole mass definition suffers from the renormalon ambiguity. The top quark pole mass cannot be defined better than $`O(\mathrm{\Lambda }_{\mathrm{QCD}})`$. Large NNLO corrections and a large shift of the $`1S`$ resonance can spoil the top quark mass measurement at the NLC.

### 4.3 The concept with the PS mass

In this subsection we discuss the calculation concept using the PS mass. We redefine the pole mass through the PS mass using the relation given in Eq. (32) without the $`1/m`$ contributions and then use the PS mass as an input parameter for our numerical analysis at LO, NLO, and NNLO. In Fig. 10 we show the top quark cross section expressed in terms of $`m_{\mathrm{PS}}(\mu _f)`$ at LO, NLO, and NNLO (like in Fig. 9) as a function of the center-of-mass energy. We take $`m_{\mathrm{PS}}(\mu _f=20\mathrm{GeV})=173.30\mathrm{GeV}`$ which corresponds to the pole mass $`m=175.05\mathrm{GeV}`$ according to Table 1. In looking at Fig. 10 we observe an improvement in the stability of the position of the first peak in comparison to the previous analysis as we go from LO to NLO to NNLO. Actually, for the three values $`\mu =15\mathrm{GeV}`$, $`30\mathrm{GeV}`$, and $`60\mathrm{GeV}`$ we obtain the maxima of the $`1S`$ peak for NNLO at $`s_{\mathrm{max}}=347.32\mathrm{GeV}`$, $`347.41\mathrm{GeV}`$, and $`347.48\mathrm{GeV}`$ while the maximal values are given by $`R_{\mathrm{max}}=1.379`$, $`1.184`$, and $`1.088`$, respectively. This demonstrated that a large variation in the renormalization scale $`\mu `$ gives rise only to a shift of about $`160\mathrm{MeV}`$ for the $`1S`$ peak position at NNLO while the variation for $`R_{\mathrm{max}}`$ is still large.

### 4.4 The concept with the $`\overline{\mathrm{PS}}`$ mass

In this subsection we discuss the calculation concept where we use the $`\overline{\mathrm{PS}}`$ mass. We redefine the pole mass through the $`\overline{\mathrm{PS}}`$ mass by using the relation given in Eq. (32) and then use the $`\overline{\mathrm{PS}}`$ mass as an input parameter for the numerical analysis at LO, NLO, and NNLO. In Fig. 11 we show the top quark cross section expressed in terms of $`m_{\overline{\mathrm{PS}}}(\mu _f)`$ at LO, NLO, and NNLO (like in Fig. 9) as a function of the center-of-mass energy. We take $`m_{\overline{\mathrm{PS}}}(\mu _f=20\mathrm{GeV})=173.28\mathrm{GeV}`$. Again we observe a very good stability of the position of the first peak as we go from LO to NLO to NNLO, similar to the one observed for the PS mass case. Studying the size of the NNLO corrections to the peak positions we conclude that the current theoretical uncertainty of the determination of the PS mass from the $`1S`$ peak position is about $`100200\mathrm{MeV}`$.

## 5 Conclusion and discussions

We introduced the so-called $`\overline{\mathrm{PS}}`$ mass which is an infrared-insensitive quark mass, naturally defined by subtracting the soft part of the quark self energy from the pole mass. We demonstrated a deep relation of this definition with the static quark-antiquark potential. At leading order in $`1/m`$ the $`\overline{\mathrm{PS}}`$ mass coincides with the PS mass, defined in a completely different manner. In contrast to the pole mass, the $`\overline{\mathrm{PS}}`$ mass is not sensitive to the non-perturbative QCD effects. Going beyond the static limit, the small normalization point introduces recoil corrections to the relation of the pole mass with the $`\overline{\mathrm{PS}}`$ mass which have been calculated in this paper. In addition we demonstrated that, if we use the PS or the $`\overline{\mathrm{PS}}`$ mass in the calculations, the perturbative predictions for the cross section become much more stable at higher orders of QCD (shifts are below $`0.1\mathrm{GeV}`$). This understanding removes one of the obstacles for an accurate top mass measurement and one can expect that the top quark mass will be extracted from a threshold scan at NLC with an accuracy of about $`100200\mathrm{MeV}`$.
Acknowledgements: We are grateful to Ratindranath Akhoury, Ed Yao and Martin Beneke for valuable discussions. O.Y. acknowledges support from the US Department of Energy (DOE). S.G. acknowledges a grant given by the DFG, FRG, he also would like to thank the members of the theory group at the Newman Lab for their hospitality during his stay.

## Appendix

The coefficients $`a_1`$ and $`a_2`$ were calculated in Refs.  and , respectively, and are given by

$`a_1`$ $`=`$ $`{\displaystyle \frac{31}{9}}C_A{\displaystyle \frac{20}{9}}T_FN_F,`$
$`a_2`$ $`=`$ $`\left({\displaystyle \frac{4343}{162}}+4\pi ^2{\displaystyle \frac{\pi ^4}{4}}+{\displaystyle \frac{22}{3}}\zeta _3\right)C_A^2\left({\displaystyle \frac{1798}{81}}+{\displaystyle \frac{56}{3}}\zeta _3\right)C_AT_FN_F`$ (A1)
$`+\left({\displaystyle \frac{20}{9}}T_FN_F\right)^2\left({\displaystyle \frac{55}{3}}16\zeta _3\right)C_FT_FN_F.`$
The first two coefficients in the expansion of the QCD $`\beta `$-function are

$`\beta _0={\displaystyle \frac{11}{3}}C_A{\displaystyle \frac{4}{3}}T_FN_F,\beta _1={\displaystyle \frac{34}{3}}C_A^2{\displaystyle \frac{20}{3}}C_AT_FN_F4C_FT_FN_F.`$ (A2)
The short distance coefficient reads

$`C(r)`$ $`=`$ $`14C_F{\displaystyle \frac{\alpha _s(\mu _h)}{\pi }}+C_2(r)C_F\left({\displaystyle \frac{\alpha _s(\mu _h)}{\pi }}\right)^2,`$
$`C_2(r)`$ $`=`$ $`A_1\mathrm{log}(r/a)+A_2\mathrm{log}(m/\mu _h)+A_3`$ (A3)
where the hard renormalization scale is taken to be the pole mass, $`\mu _h=m`$, $`a=e^{2\gamma _E}/2m`$, and

$`A_1`$ $`=`$ $`\pi ^2\left(C_A+{\displaystyle \frac{2C_F}{3}}\right),A_2=2\beta _0,`$
$`A_3`$ $`=`$ $`C_2^AC_F+C_2^{NA}C_A+C_2^LT_FN_F+C_2^HT_F`$ (A4)
where $`N_F=5`$ is the number of the light flavors, and

$`C_2^A`$ $`=`$ $`{\displaystyle \frac{39}{4}}\zeta _3+\pi ^2\left({\displaystyle \frac{4}{3}}\mathrm{ln}2{\displaystyle \frac{35}{18}}\right),`$
$`C_2^{NA}`$ $`=`$ $`{\displaystyle \frac{151}{36}}{\displaystyle \frac{13}{2}}\zeta _3+\pi ^2\left({\displaystyle \frac{179}{72}}{\displaystyle \frac{8}{3}}\mathrm{ln}2\right),`$
$`C_2^L`$ $`=`$ $`{\displaystyle \frac{11}{9}},C_2^H={\displaystyle \frac{44}{9}}{\displaystyle \frac{4}{9}}\pi ^2.`$ (A5)


format: text

text:
# Adsorbate aggregation and relaxation of low-frequency vibrations

## I Introduction

Molecules or atoms adsorbed on surfaces of materials can oscillate about the adsorption bond. If this motion corresponds, in the absence of the adsorption bond, to free translation parallel to the surface, then this mode is called a frustrated translation (FT). It plays an important role in such dynamical processes as reactivity, surface diffusion, and desorption. For many systems, especially with a metallic substrate, the frequency of the FT mode is so low that it readily couples to low-energy excitations in the substrate, which in the case of metals include both phonons and electron-hole pair excitations. These bulk excitations provide decay channels for FT modes, which thus acquire a finite lifetime $`\tau `$ and become resonances. Associated with this lifetime is the resonance width, or damping rate, $`\gamma =1/\tau `$. In many cases the electronic contribution to the damping rate is small compared with the phononic contribution.

We recently developed a general theory of phonon-mediated FT relaxation that can be applied to any adsorbate overlayer pattern. Using this theory we derived simple FT damping rate laws for three important classes of adsorbate patterns: an isolated adsorbate, an ordered overlayer, and a randomly disordered overlayer. Although there are no experimental results to compare to for isolated adsorbates, the agreement with experiment for ordered and disordered overlayers is excellent.

However, these three structural models do not cover all possible experimental conditions. Depending on the interactions of the adsorbates with the substrate and with each other, their distribution on the surface can differ strongly from the periodic or random uncorrelated cases. Under certain conditions the adsorbates may gather into aggregates even at very low coverages. In this paper we consider two possible examples of such adsorbate aggregation: (i) adsorbates accumulating along the edges of surface steps, and (ii) adsorbates forming two-dimensional (2D) islands. We assume that the average coverage is low enough that these aggregates are well separated from each other and can be treated as isolated. For each of these cases we derive a corresponding damping rate formula from our general result.

Another goal of this paper is to facilitate the evaluation of damping rates for a disordered overlayer. The main difficulty here is that one needs to integrate the Green’s function for a semi-infinite elastic medium over the in-plane wavevector at a dense array of frequencies. Regrettably, no closed-form expression for this integral exists, and it must be computed numerically. In this paper we present a good analytical approximation of this integral for elastic media of isotropic and cubic symmetries. We have tested our approximation formula for 88 different materials with cubic symmetry, comparing it against well-converged numerical calculations. The difference between the analytical approximation and the reference numerical values never exceeds 6%, and for a vast majority of the materials is under 4%.

The remainder of this paper is organized as follows. In Section II we briefly review our theory of phonon-mediated FT damping. For details of the derivation, the reader is referred to Ref. . This theory is then applied in Sections III and IV to the cases of adsorbates aggregated along step edges and in islands, respectively. In Section V we present our approximation formula for the $`\stackrel{}{k}`$-space integral of the substrate Green’s function, thus facilitating the evaluation of the damping rate of a disordered overlayer. Finally, in Section VI we summarize and draw conclusions.

## II Review of theory

The phonon contribution to the damping of low-frequency adsorbate vibrations was considered earlier by Persson and Ryberg  (PR) and by Hall, Mills, and Black  (HMB). Persson and Ryberg considered a single, isolated adsorbate coupled to a semi-infinite isotropic elastic substrate. Since direct adsorbate-adsorbate interactions are weak for all but the densest coverages, this case was presumed applicable for most experimental conditions. Their model gives a damping rate law that varies with the vibrational frequency as $`\omega _0^4`$. However when applied to various experimental systems, the PR model predicts damping rates that are much smaller than the measured values. Furthermore, the PR model has no explicit coverage dependence; coverage enters the damping-rate law only implicitly through its effect on the resonance frequency. The applicability of the model is limited by extremely low coverage ($`1\%`$). Hall, Mills, and Black considered a periodic overlayer of rare-gas molecules physisorbed at an isotropic substrate. The damping rate implied by their result is proportional to the coverage and to the square of the resonant frequency. The periodic distribution of the adsorbates over the surface implied by the HMB model generally occurs at high coverage ($`25\%`$). In general, most adsorbate systems have structures that differ from both of these special cases.

In our recent work  we developed a general theory for a completely arbitrary overlayer structure. We considered an array of $`N`$ non-interacting adsorbates of mass $`m`$ located at arbitrary positions $`\{\stackrel{}{R}_a\}`$ at the surface of an anisotropic, semi-infinite elastic substrate of mass density $`\rho `$ and elastic modulus tensor $`C_{ijkl}`$, occupying the half-space $`z<0`$. The adsorbates are coupled to the substrate by harmonic springs of frequency $`\omega _0`$ and are allowed to oscillate along one of the Cartesian axes $`\alpha `$. The adsorbates do not interact with each other directly. We showed that the damping rate for this system is described by

$$\gamma =\frac{m}{\rho }\omega _0^2\frac{d^2k}{(2\pi )^2}\left[\omega _0\mathrm{Im}D_{\alpha \alpha }(\stackrel{}{k},\omega _0)\right]S(\stackrel{}{k})\frac{m}{\rho }\omega _0^2I(\omega _0),$$
(1)
where $`D_{ij}(\stackrel{}{k},\omega )`$ is the in-plane Fourier transform of the substrate elastic Green’s function evaluated at the surface ($`z=z^{}=0`$), and $`S(\stackrel{}{k})`$ is the structure factor defined as

$$S(\stackrel{}{k})=\frac{1}{N}\underset{a,b}{}\mathrm{exp}[i\stackrel{}{k}(\stackrel{}{R}_a\stackrel{}{R}_b)],$$
(2)
where the angle brackets denote the appropriate ensemble average. The integration in Eq. (1) is taken over the bare-surface Brillouin zone. The integral $`I(\omega )`$ in Eq. (1) is a convolution of two quantities, one of which (inside the square brackets) describes all the substrate properties, while $`S(\stackrel{}{k})`$ provides all the information on how the adsorbates are distributed over the surface. Except for a unitless factor of order unity, the quantity in the square brackets is equal to the $`\stackrel{}{k}`$-resolved density of elastic modes at the substrate surface for frequency $`\omega _0`$. Any possible explicit dependence of the damping rate upon the coverage or the structural order comes from the structure factor $`S(\stackrel{}{k})`$.

We previously applied this general result to three different structure models: an isolated adsorbate, a periodic array of adsorbates (appropriate at high coverages), and a random array of adsorbates (appropriate at low coverages). For each of these models, the structure factor has a simple algebraic form. Their substitution into Eq. (1) leads to

$$\gamma =\frac{m}{\rho }\omega _0^4\beta $$
(3)
for an isolated adsorbate,

$$\gamma =\frac{m}{\rho }\omega _0^2\frac{\theta }{A_0c_{}^{(\alpha )}}$$
(4)
for a periodic overlayer, and

$$\gamma =(1\theta )\frac{m}{\rho }\omega _0^4\beta +\frac{m}{\rho }\omega _0^2\frac{\theta }{A_0c_{}^{(\alpha )}}$$
(5)
for a random overlayer, where $`\theta `$ is the coverage, $`A_0`$ is the area per surface site, $`c_{}^{(\alpha )}`$ is the speed of an $`\alpha `$-polarized sound wave propagating perpendicular to the surface, and $`\beta `$ is defined according to

$$\beta =\frac{1}{\omega }\frac{d^2k}{(2\pi )^2}\mathrm{Im}D_{\alpha \alpha }(\stackrel{}{k},\omega )$$
(6)
The integral in Eq. (6) is proportional to the density of $`\alpha `$-polarized states divided by frequency. For a 3D elastic medium the density of $`\alpha `$-polarized states is proportional to $`\omega ^2`$, and thus $`\beta `$ is independent of frequency $`\omega `$. The evaluation of this parameter is discussed in Section V. The damping rates given by Eqs. (3) and (4) are just the anisotropic versions of the PR  and HMB  results, respectively. No experimental results are available for an isolated adsorbate. However, we compared the other two models with experiments on CO on Cu(001) by applying Eq. (4) to the case of an ordered overlayer at 50% coverage, and Eq. (5) to the case of a disordered overlayer at 3% coverage. Both models predict the lifetimes in excellent agreement with the measured values. Equation (5) was independently derived later by Persson et al.  for the case of an isotropic substrate. They also considered the case of an isolated adsorbate of finite size.

## III Adsorbates aligned along step edges

In this Section we study the case of adsorbates located at step edges. This geometry is important since steps are always present at surfaces and they generally attract adsorbates to their edges. While the macroscopic coverage is low for this case, the local coverage is high.

In our analysis we ignore the influence of the steps on the substrate Green’s function. The only effect of steps that we consider is to set the distribution of adsorbates. We assume that the adsorbates form one-dimensional (1D), locally periodic arrays. Starting from Eq. (1) we derive the damping rate for an infinite, 1D periodic chain of adsorbates, and then we extend the derivation for a finite, locally periodic chain.

### A One-dimensional infinite chain

We assume that the adsorbates are arranged in an infinite, 1D, periodic array along the $`x`$ axis. The adsorbate positions are given by $`\stackrel{}{R}_a=(X_a,0)`$, and the nearest-neighbor distance is $`d`$. The 1D coverage of the step edge can be defined as $`\theta _{1D}=a_0/d`$, where $`a_0`$ is the nearest-neighbor distance between the substrate atoms along the step edge. Then the structure factor $`S(\stackrel{}{k})`$ takes the form

$$S(\stackrel{}{k})=\frac{2\pi }{d}\underset{n}{}\delta (k_xG_n),$$
(7)
where $`\{G_n\}`$ is the set of 1D reciprocal vectors defined as $`G_n=2\pi n/d`$. Substitution of Eq. (7) into Eq. (1) leads to

$$\gamma =\frac{m}{\rho }\omega _0^3\frac{1}{d}\frac{dk_y}{2\pi }\underset{n}{}\left[\mathrm{Im}D_{\alpha \alpha }(G_n,k_y;\omega _0)\right].$$
(8)
Since there are no resonant phonon states available in the substrate for in-plane wavevectors $`|G_n|>\omega _0/c_R`$, where $`c_R`$ is the speed of the Rayleigh wave, the imaginary part of the Green’s function in Eq. (8) vanishes for these values of $`G_n`$, and thus these terms do not contribute to the sum. Similar to the case of a two-dimensional (2D) periodic overlayer discussed in detail in our earlier work, it turns out that for all 1D coverages presently accessible to experiment ($`\theta _{1D}>\theta _{1D}^{(c)}10\%`$), only the term with $`G_n=0`$ contributes to the damping rate. Moreover, the integral of this term in Eq. (8) does not depend on frequency $`\omega _0`$, and we can define a constant $`\zeta _\alpha `$ according to

$$\zeta _\alpha =\frac{dk_y}{2\pi }\left[\mathrm{Im}D_{\alpha \alpha }(0,k_y;\omega _0)\right].$$
(9)
It turns out that $`\zeta _x`$ (vibration parallel to the step) differs significantly from $`\zeta _y`$ (vibration perpendicular to the step). However, for substrates with isotropic or cubic symmetry, these values, while different from each other, are independent of the orientation of the step. It is, therefore, appropriate to use the notation $`\zeta _{}`$ for oscillations along the step and $`\zeta _{}`$ for oscillations perpendicular to the step but still in the plane. We can now write the damping rate formula compactly as

$$\gamma _,=\frac{m}{\rho }\omega _0^3\frac{\zeta _,}{d}.$$
(10)
Thus, the oscillations along the step are damped at a different rate than the oscillations perpendicular to the step. Typically, the damping rate is higher for vibrations along the step.

Experimentally measured damping rates can be related to Eq. (10) in two possible ways. If the frequencies for both vibrational modes are close, then the measured damping rate will correspond to an average of $`\gamma _{}`$ and $`\gamma _{}`$. Otherwise, both modes can be resolved experimentally, each exhibiting its own damping rate.

### B One-dimensional finite chain

Here we consider a finite, locally periodic chain of $`N`$ adsorbates along one of the step edges. As before, the $`x`$-axis is along the step, and the nearest-neighbor distance in the chain is $`d`$. The total length of the chain is $`L=Nd`$. The structure factor corresponding to this configuration is given by

$$S(\stackrel{}{k})=N\left[\underset{n}{}\frac{\mathrm{sin}\left[(k_xG_n)L/2\right]}{\left[(k_xG_n)L/2\right]}\right]^2.$$
(11)
Most of the weight of this function is in the peaks centered at the reciprocal vectors $`G_n=2\pi n/d`$. Even at relatively small values of $`N`$ the cross terms in Eq. (11) may be neglected, and the square of the sum can be approximated by a sum of squares

$$S(\stackrel{}{k})=N\underset{n}{}\left[\frac{\mathrm{sin}\left[(k_xG_n)L/2\right]}{\left[(k_xG_n)L/2\right]}\right]^2.$$
(12)
Equations (11) and (12) are compared for a chain of $`N=6`$ adsorbates separated by $`d=5.54`$ Å(twice the nearest-neighbor distance for platinum) in Fig. 1. The solid line corresponds to the structure factor computed according to Eq. (11), while the dashed line represents $`S(\stackrel{}{k})`$ given by Eq. (12). As $`N`$ increases, the main peaks become narrower (the width of the peaks is approximately $`2\pi /L`$), and the difference between the two formulas decreases. In the limit $`N\mathrm{}`$, both expressions take the form of Eq. (7).

Finally, we show that the damping rate for a finite chain of the adsorbates can be computed approximately from Eq. (10) for large enough $`N`$ (where “large enough” may be as small as 6). First, we assume that all the weight in the structure factor comes from the main peaks. Then, repeating the argument we gave for the infinite chain we conclude that if $`\omega _0/c_R<2\pi (11/N)/d`$, only the peak centered at $`G=0`$ contributes to the integral in Eq. (1), and the damping rate becomes

$$\gamma =\frac{m}{\rho }\omega _0^3\frac{dk_x}{2\pi }N\left[\frac{\mathrm{sin}\left(k_xL/2\right)}{(k_xL/2)}\right]^2\frac{dk_y}{2\pi }\left[\mathrm{Im}D_{\alpha \alpha }(\stackrel{}{k},\omega _0)\right].$$
(13)
If the peak is narrow enough, the imaginary part of the Green’s function does not change significantly over its width (the region $`k_x<2\pi /L`$). Thus, for each value of $`k_y`$ the Green’s function can be replaced by its value at $`k_x=0`$. The integral over $`k_x`$ can now be carried out explicitly, and Eq. (10) follows.

## IV Adsorbates forming islands

In this Section we analyze the damping of the FT mode when adsorbates form 2D islands at the substrate surface. As with 1D step-edge aggregation, this case corresponds to low macroscopic coverage, but high local coverage. We assume that each island is a locally periodic array of $`N`$ adsorbates with local coverage $`\theta _{\mathrm{loc}}`$ and area $`A`$. We also assume that the distance between islands is large enough that the islands can be treated as isolated.

Unlike in the previous Section, there is no need here for considering the case of an infinite island, because an infinite island is the same as a periodic overlayer with $`\theta _{\mathrm{loc}}=\theta `$, and its FT damping rate is given by Eq. (4). The problem of a finite island appears to be more complicated than that of a finite 1D chain at a step edge. The reason is that, apart from its local coverage and size, the island is also characterized by its shape. It is convenient to represent the shape of the island by a function $`f(\stackrel{}{r})`$ defined as

$$f(\stackrel{}{r})=\{\begin{array}{cc}1\text{if }\stackrel{}{r}\text{ is inside island}\hfill & \\ 0\text{otherwise}\hfill & \end{array}$$
(14)
It is straightforward to show that the structure factor can be expressed in terms of the Fourier transform $`f(\stackrel{}{k})`$ of this function as

$$S(\stackrel{}{k})=\frac{N}{A^2}\underset{\stackrel{}{G},\stackrel{}{G}^{}}{}f(\stackrel{}{k}\stackrel{}{G})f(\stackrel{}{k}\stackrel{}{G}^{}),$$
(15)
where summation is over reciprocal vectors defined by the local periodicity of the adsorbate distribution inside the island.

We further assume that the island is a convex set with similar dimensions in all directions (i.e. its width is of the same order of magnitude as its length). If the number of adsorbates in the island is not too small (the meaning of which is discussed below), then the function $`f(\stackrel{}{k})`$ will have a strong peak at $`\stackrel{}{k}=0`$, with most of the function’s weight within the peak. Analytical examples of this function include

$$f_{\mathrm{rect}}(\stackrel{}{k})=A\frac{\mathrm{sin}\left[k_xa/2\right]}{\left[k_xa/2\right]}\frac{\mathrm{sin}\left[k_yb/2\right]}{\left[k_yb/2\right]}.$$
(16)
for a rectangular island with dimensions $`a`$ along the $`x`$-axis and $`b`$ along the $`y`$-axis, and

$$f_{\mathrm{circ}}(\stackrel{}{k})=2A\frac{J_1(|\stackrel{}{k}|R)}{|\stackrel{}{k}|R}$$
(17)
for a circular island of radius $`R`$, where $`J_1(x)`$ is the first-order Bessel function.

We now follow the arguments of the 1D case. First, we neglect the cross-terms in Eq. (15) and write it as

$$S(\stackrel{}{k})=\frac{N}{A^2}\underset{\stackrel{}{G}}{}\left[f(\stackrel{}{k}\stackrel{}{G})\right]^2.$$
(18)
Then, we assume that all the weight in the structure factor comes from the main peaks. If the resonant frequency of the FT mode $`\omega _0`$ is small enough that $`\omega _0/c_R<G_1\mathrm{\Delta }`$, where $`G_1`$ is the magnitude of the smallest reciprocal vector, and $`\mathrm{\Delta }`$ is the maximal width of the main peak in $`S(\stackrel{}{k})`$, then only the $`\stackrel{}{G}=0`$ term will contribute to the integral in Eq. (1). Thus this equation now takes the form

$$\gamma =\frac{m}{\rho }\omega _0^2\frac{d^2k}{(2\pi )^2}\left[\omega _0\mathrm{Im}D_{\alpha \alpha }(\stackrel{}{k},\omega _0)\right]\frac{N}{A^2}\left[f(\stackrel{}{k})\right]^2.$$
(19)
If the the peak of $`[f(\stackrel{}{k})]^2`$ is narrow enough, the imaginary part of the Green’s function does not change significantly over its width (i.e., the region $`|\stackrel{}{k}|<\mathrm{\Delta }`$), and it can be replaced by its value at $`\stackrel{}{k}=0`$ given by $`\mathrm{Im}D_{\alpha \alpha }(0,\omega _0)=1/(c_{}^{(\alpha )}\omega _0)`$, and taken out of the integral. The remaining integral can be evaluated with the help of the Parseval’s theorem as

$`{\displaystyle \frac{d^2k}{(2\pi )^2}\frac{N}{A^2}\left[f(\stackrel{}{k})\right]^2}={\displaystyle d^2r\frac{N}{A^2}\left[f(\stackrel{}{r})\right]^2}={\displaystyle \frac{N}{A}}={\displaystyle \frac{\theta _{\mathrm{loc}}}{A_0}},`$

where, as before, $`A_0`$ is the area of the surface unit cell. Substituting this result into Eq. (19) we obtain

$$\gamma =\frac{m}{\rho }\omega _0^2\frac{\theta _{\mathrm{loc}}}{A_0c_{}^{(\alpha )}},$$
(20)
which is exactly the damping rate of a periodic overlayer given by Eq. (4) with the total coverage $`\theta `$ being replaced by the local coverage $`\theta _{\mathrm{loc}}`$.

Finally, we need to determine for which values of $`N`$ the approximation given by Eq. (20) really works. We consider, as an example, an isolated island of hexagonal shape on the Pt(111) surface with local coverage $`\theta _{\mathrm{loc}}=1/4`$. This is about the lowest local coverage observed for islands of small molecular adsorbates. If the approximation given by Eq. (20) holds for this coverage, it will also hold for higher coverages. We have computed the integral $`I(\omega )`$ in Eq. (1) for hexagonal islands of different sizes, and, in Fig. 2, plot $`I(\omega )`$ vs. $`\omega `$ for each island size. The value of the integral for the two limiting cases of an isolated adsorbate and an infinite periodic overlayer are shown by the dashed and solid lines, respectively. The dotted, dash-dotted, and dash-double-dotted lines correspond to the cases of one, two, and three hexagonal shells of adsorbates around the center adsorbate, respectively. The inset to Fig. 2 characterizes the islands by their circumscribing radius $`R`$, and the number of adsorbates, $`N`$, within the island.

The behavior of $`I(\omega )`$ shown in Fig. 2 can readily be explained in terms of resonant phonon wavelengths at frequency $`\omega `$. For small frequencies the phonon wavelengths are much greater than the island size, and the adsorbates in the island behave like a single isolated adsorbate of mass $`Nm`$. Thus, for small frequencies we have $`I(\omega )=N\beta \omega ^2`$ in accord with Eq. (3). Mathematically, this result can be obtained from Eq. (19) by setting $`f(\stackrel{}{k})=f(0)=A`$. Frequencies for which the resonant phonon wavelengths become comparable with the island size correspond to the cross-over region where this quadratic behavior changes. When the phonon wavelengths are smaller than the island size, but still significantly greater than the nearest-neighbor distance between the adsorbates, only phonons with in-plane wavevector component, $`\stackrel{}{k}`$, near zero couple strongly to the FT mode due to the $`\stackrel{}{G}=0`$ peak in the structure factor given by Eq. (18). Terms with $`\stackrel{}{G}0`$ have no corresponding resonant phonons at these frequencies. In this frequency region we expect $`I(\omega )`$ to follow closely its limiting value of a full periodic overlayer. Finally, when $`\omega `$ is large enough that the phonon wavelengths become comparable to the inter-adsorbate nearest-neighbor distance, non-zero $`\stackrel{}{G}`$-vectors in the structure factor become available channels for resonant phonon coupling. However, for realistic local coverages ($`\theta _{\mathrm{loc}}1/4`$) these wavevectors lie at the boundary of the bare-surface Brillouin zone, where continuum models are no longer applicable. For the case $`\theta _{\mathrm{loc}}=1`$, this second cross-over region does not exist, because all non-zero $`\stackrel{}{G}`$-vectors lie outside the surface Brillouin zone.

From Fig. 2 one can determine the range of resonant FT frequencies for which the approximation given by Eq. (20) holds. Even for the island with only one hexagonal shell ($`N=7`$), there is such a frequency range: $`I(\omega )`$ stays close to its limiting value of a periodic overlayer for a range of frequencies between about 3 and 7 meV. For larger islands this range expands to frequencies between 2 and 11 meV.

We can analytically estimate this frequency range for an island of size $`D`$. The low-frequency threshold corresponds to a transverse wave with wavevector magnitude equal to the width, $`\mathrm{\Delta }`$, of a peak in the structure factor, which can be estimated as $`2\pi /D`$. The frequency of such a wave is $`\omega _{\mathrm{low}}=2\pi c_\mathrm{T}/D`$. The high frequency threshold corresponds to a transverse wave with wavevector magnitude equal to $`G_1\mathrm{\Delta }`$, which is approximately $`2(D/d1)\pi c_\mathrm{T}/D`$, where $`d`$ is the inter-adsorbate nearest-neighbor distance, as before. The frequency of this wave is $`\omega _{\mathrm{high}}=(D/d1)\omega _{\mathrm{low}}`$. For a three-shell hexagonal island ($`N=37`$) with $`\theta _{\mathrm{loc}}=1/4`$ on the Pt(111) surface, we have $`d=5.54`$ Å, $`D/d=6`$, and $`c_\mathrm{T}=11`$ Å meV, which leads to the estimates $`\omega _{\mathrm{low}}=2.1`$ meV and $`\omega _{\mathrm{high}}=10.5`$ meV. This frequency range is in excellent agreement with the corresponding plot in Fig. 2. Thus, the example of a hexagonal island on the Pt(111) surface demonstrates that, when most of the adsorbates form widely separated islands, the damping rate can be computed from Eq. (20) for FT resonant frequencies within a broad range.

## V Approximation formula for $`\beta `$ parameter

The main difficulty researchers may encounter in using Eq. (5) to determine the damping rate for low-coverage systems is the necessity to compute the substrate Green’s function and the integral in Eq. (6) numerically. In this Section we present an analytical approximation to the parameter $`\beta `$ defined by Eq. (6). The approximation is developed for the cases of systems with isotropic or cubic symmetry.

We start from the observation that $`\beta `$ must be a function of the substrate parameters of the model, namely, the mass density $`\rho `$ and the elastic constants $`c_{ij}`$. Generally, $`i,j=1,\mathrm{},6`$, but we include only independent constants characteristic of a given symmetry. This function must obey certain global constraints, however. It follows from the equation of motion and the boundary conditions that any quantity depending on the model inputs may depend only on the ratios $`c_{ij}/\rho `$, which we denote as $`k_{ij}`$. Because $`\beta `$ does not depend on the frequency $`\omega `$, one can show from the equations of motion that $`\beta `$, as a function of density-normalized elastic constants $`k_{ij}`$, must satisfy a scaling property given by

$$\beta (vk_{ij})=v^{3/2}\beta (k_{ij}),$$
(21)
where $`v`$ is an arbitrary number. This scaling property reduces the number of independent variables in the function $`\beta (k_{ij})`$ by one. We have found an excellent analytical approximation to this function for the cases of isotropic and cubic substrates.

### A Isotropic substrate

An isotropic medium is characterized by only two density-normalized elastic constants: $`k_{11}`$ and $`k_{44}`$. The transverse and the longitudinal speeds of sound are related to these constants by $`c_\mathrm{T}=\sqrt{k_{44}}`$ and $`c_\mathrm{L}=\sqrt{k_{11}}`$, respectively. Taking into account the scaling property given by Eq. (21), we can express $`\beta ^2`$ as

$$\beta ^2=k_{44}^3g\left(\frac{k_{11}}{k_{44}}\right),$$
(22)
where $`g(x)`$ is an unknown dimensionless function of one variable.

To elucidate the behavior of $`g(x)`$, we have computed $`\beta `$ numerically for a wide range of $`k_{11}/k_{44}`$. The result of these calculations is shown by squares in Fig. 3. From Fig. 3 we obtain that (i) $`g(x)`$ goes to zero at $`x=1`$ as $`(x1)^\alpha `$, where $`1<\alpha <2`$, and (ii) $`g(x)`$ asymptotically approaches a constant as $`x\mathrm{}`$. We find that the function

$$g(x)=P\frac{(x1)^{3/2}}{3/2+(x1)^{3/2}}$$
(23)
almost perfectly approximates the numerical data (solid line in Fig. 3). The parameter $`P`$ is dimensionless and has an optimal value of $`P=9.14\times 10^{11}`$.

### B Cubic substrate

The isotropic approximation is not correct for most single-crystal materials. We therefore extend our approximation formula for $`\beta `$ to describe substrate materials with cubic symmetry, since these are the most commonly used in experiments. Crystals with cubic symmetry have three independent density-normalized elastic constants: $`k_{11}`$, $`k_{12}`$, and $`k_{44}`$. The isotropic limit can be obtained from cubic symmetry by setting $`k_{12}=k_{11}2k_{44}`$. Taking advantage of the scaling property given by Eq. (21), we can write

$$\beta ^2=k_{12}^3h(\frac{k_{11}}{k_{12}},\frac{k_{44}}{k_{12}}),$$
(24)
where $`h(x,y)`$ is an unknown function of two variables. Generally, the form of the function $`h(x,y)`$ may depend on the particular choice of the surface indices. However, numerical calculations of $`\beta `$ for different low-index surfaces of the same material show that $`\beta `$ is essentially independent of the surface orientation with respect to the crystallographic axes. Thus, the same function $`h(x,y)`$ works well for any low-index surface.

First, we study how $`h(x,y)`$ depends on $`x`$ by computing $`\beta `$ numerically as a function of $`k_{11}/k_{12}`$ while keeping $`k_{44}/k_{12}`$ fixed. We find that $`h(x,y)`$ depends linearly on $`x`$ when $`x1`$ and goes to zero at $`x=1`$ as $`(x1)^{3/2}`$. We choose the fitting function to have the form

$$h(x,y)=b(y)\frac{(x1)^{3/2}}{y+(x1)^{1/2}}$$
(25)
which satisfies these asymptotic requirements.

The next step is to find a fitting function for $`b(y)`$. This can be accomplished by requiring that $`\beta `$ for the cubic substrate given by Eq. (24) go to the isotropic limit described by Eq. (22) when $`k_{12}=k_{11}2k_{44}`$. This requirement uniquely specifies that the function $`b(y)`$ has the form

$$b(y)=\frac{P}{2}y^2\left[1+\left(\frac{y}{2}\right)^{1/2}\right]\frac{(y^1+1)^{3/2}}{3/2+(y^1+1)^{3/2}},$$
(26)
where the parameter $`P`$ is the one determined for the isotropic case, namely $`P=9.14\times 10^{11}`$.

To check the quality of the approximation given by Eqs. (24) and (26), we have computed $`\beta `$ according to this approximation and compared it with $`\beta `$ obtained by numerical calculations for 88 different materials with cubic symmetry. The error in the approximated value of $`\beta `$ never exceeds 6% for this set of materials, and for the vast majority of cases it is below 4%. Table I shows this comparison for several important substrate materials. It should be noted that a given error in $`\beta `$ results in a smaller error in the predicted damping rate, $`\gamma `$, due to the $`\beta `$-independent term in Eq. (5).

## VI Conclusion

In this paper we have applied the general result of our continuum elastic model  of phonon-mediated adsorbate vibrational relaxation to derive analytic formulas for the FT damping rate when the adsorbates form aggregates on the surface. We have considered two possible forms of adsorbate aggregation: (i) when the adsorbates are aligned along step edges, and (ii) when the adsorbates form well separated islands. We have derived the expression for the damping rate of an infinite periodic chain of the adsorbates, and then shown that the damping rate of a finite locally periodic chain is described by the same expression for all chains above a small threshold length. For the case of an isolated island we have shown that, under mild restrictions on the island geometry, there is a wide range of possible FT resonant frequencies for which the damping rate is essentially the same as for a periodic overlayer with coverage equal to the island’s local coverage.

Finally, to facilitate the use of our previously derived expression for the FT damping rate at low coverage, we have developed an analytical formula which approximates the parameter $`\beta `$ appearing in the damping law. The formula is valid for both isotropic and cubic materials and achieves remarkable accuracy. We have tested the formula on 88 different cubic materials (all of whose elastic constants were readily available) and have found that our analytic approximation predicts $`\beta `$ with accuracy better than 6%.

## Acknowledgments

Financial support for this research was provided by the National Science Foundation under Grant No. DMR 97-02514, the Air Force Office of Scientific Research, Air Force Materiel Command, USAF, under Grant No. F49620-00-1-0170, and the University of Georgia Research Foundation. SPL acknowledges support of the Donors of The Petroleum Research Fund, administered by the American Chemical Society. AMR acknowledges support of the Alfred P. Sloan Foundation.