text:
# Structural transitions in scale-free networks

## Abstract

Real growing networks like the WWW or personal connection based networks are characterized by a high degree of clustering, in addition to the small-world property and the absence of a characteristic scale. Appropriate modifications of the (Barabási-Albert) preferential attachment network growth capture all these aspects. We present a scaling theory to describe the behavior of the generalized models and the mean field rate equation for the problem. This is solved for a specific case with the result $`C(k)1/k`$ for the clustering of a node of degree $`k`$. Numerical results agree with such a mean-field exponent which also reproduces the clustering of many real networks.

In diverse fields of scientific interest underlying network structures can be recognized, which provide a unifying concept of investigation general . Examples range from biology (metabolic networks metabolic , protein nets in the cell protein ) through sociology (movie actor relationships emergence\_of\_scaling , coauthor networks newman , sexual nets liljeros ) to informatics (Internet faloutsos , WWW WWW ). In all these examples it is easy to identify the constituents of the problem with the nodes of a graph and their relationships with directed or undirected links. During the last few years a great deal of information has accumulated about such structures. Three apparent features seem to characterize them rather robustly: i) a high degree of clustering, i.e., if nodes $`A`$ and $`B`$ are linked to node $`C`$ then there is a good chance that $`A`$ and $`B`$ are also linked; ii) the “Small World” property, i.e., the expected number of links needed to reach from one arbitrarily selected node another one is low; iii) the absence of a characteristic scale, which often appears so that the distribution $`P(k)`$ of the degrees $`k`$ of nodes follows a power law.

It has been noticed from the beginning that clustering in real networks is an essential and an almost ubiquitous feature. It measures the deviation from a structure with vanishing correlations, and it has been used to describe the tendency of networks to form cliques or tightly connected neighborhoods. As an organizing principle, this is most obvious in social networks, where connections are usually created by personal acquaintances, like in the scientific collaboration network. Considerable clustering has also been found in networks of more diverse nature. Prime examples are the WWW, metabolic and protein interaction networks, the actor network, the power grid of the United States, the semantic web of English words hierarchical\_organization , and the backbone of the Internet on both the autonomous system and the router level pastor ; goh . The number of entries in this list is on the rise as new disciplines are being taken under consideration and raw data are made available. A comprehensive examination of a variety of real networks clustering can be found in Ref. hierarchical\_organization . In real networks, as a combination of the properties i) and iii), the clustering coefficient as a function of the degree of the nodes often follows a power law: $`C(k)k^\alpha `$. The value of $`\alpha `$ is in many networks close to 1.

In 1998 Watts and Strogatz created an interesting family of models: introducing a rather low portion of random links between arbitrarily selected pairs of nodes in a regular lattice has the consequence that property ii) gets fulfilled while clustering does not decrease considerably, assuring i) footnote ; watts\_strogatz . However, the distribution of the degrees of nodes shows a characteristic peak instead of the required power law. Barabási and Albert (BA) realized that in the examples mentioned at the beginning an important aspect is that the networks are created by growth. BA proposed preferential attachment (PA) as a growth rule: the new nodes are linked to the old ones with a probability proportional to their the actual degree emergence\_of\_scaling . The structures obtained this way are scale free and have the Small World property. In spite of capturing important aspects of growing networks, the clustering tends rapidly to a constant as a function of the degree $`k`$ and vanishes in the thermodynamic limit.

Recently, attempts have been undertaken to modify the PA network growth models so as to increase clustering. In these models a mechanism, controlled by a new parameter, is introduced to take into account the effect that “friends of friends get friends”. Indeed, it has been possible to create models which have all the three properties i)-iii) tunable\_clustering ; hierarchical\_organization ; growing\_scale\_free .

The aim of this paper is to present a general framework for the transition from a PA graph with zero clustering to still scale free graphs with $`C(k)k^\alpha `$, and to give a corresponding mean-field (MF) and rate-equation theory. As an example we will take the Holme-Kim model tunable\_clustering  (a modified BA one) for which the MF rate equations can be solved exactly, leading to $`\alpha =1`$. This is also shown to describe the simulations very well. At the end, we discuss the assumptions one needs to make, and how these reflect the behavior of clustering in general.

We start from the simplest undirected BA model: a new node $`j`$ with $`m`$ links is added to the system at (discrete) time $`t`$. A link from node $`j`$ to node $`i`$ is drawn with probability $`k_i/k_i`$. It is known that the the average clustering at node $`i`$ is independent of the degree $`k_i`$ growing\_scale\_free :

$`C(k_i)={\displaystyle \frac{m1}{8}}{\displaystyle \frac{(\mathrm{log}N)^2}{N}},`$ (1)
i.e., it is inversely proportional to the number $`N`$ of nodes (with a logarithmic correction) wrongexpr . For the generalization of the BA model with enhanced clustering, we have a parameter $`p`$ representing an imposed tendency to form triangles on the graph. It is chosen such that at $`p=0`$ the original BA model is recovered.

We propose as a scaling ansatz to describe the clustering coefficient $`C`$ as a function of the degree $`k`$, the number of nodes $`N`$ and the parameter $`p`$:

$`C(k,N,p)=N^1f\left({\displaystyle \frac{k}{k^{}(N,p)}}\right),`$ (2)
where $`f(x)`$ is a scaling function with $`f(x)\mathrm{const}.`$ for $`x1`$ and $`f(x)\mathrm{x}^\alpha `$ for $`x1`$ and the behavior in Eq. (1) is already taken into account by fixing the exponent of the prefactor of $`f`$. The characteristic degree $`k^{}`$ is a monotonously increasing function of $`N`$ for fixed $`p`$ and it should decrease as $`p`$ goes to zero. A natural assumption is then:

$`k^{}(N,p)N^\gamma p^\delta .`$ (3)
As for small $`k`$ the clustering $`C`$ in Eq. (2) should go like $`k^\alpha `$ and become independent of $`N`$, we have $`\gamma =1/\alpha `$. The exponent $`\delta \alpha `$ describes, how for $`N\mathrm{}`$ the clustering $`C`$ approaches its limiting value zero as $`p`$ goes to zero. If we accept that in most cases $`\alpha =1`$, there is one exponent to be determined, say $`\delta `$. We now clarify the origin of $`\alpha =1`$ and $`\delta =1`$ for the model employed.

For this purpose we write down the rate equations for the clustering in a general form. We thus need to consider the average rate of change

$$_tn_i=R(k_i,p)\underset{n\mathrm{\Omega }}{}R(k_n,p),$$
(4)
where $`n_i`$ is the number of connected neighbors of site $`i`$, and $`C_i=n_i/(k_i(k_i1)/2)`$. Here $`R`$ is the rate at which $`i`$ gets new links, and we allow, in analogy with the scaling ansatz presented above, the rate to depend on both the degrees of the node in question and the parameter $`p`$. This can be “annealed” or “quenched”, depending on whether the parameter describes stochastic rules (as in the example below) or a fixed property of each node $`i`$. E.g., $`R`$ can simply follow from the preferential attachment rule. $`\mathrm{\Omega }`$ is the set of neighbors of node $`i`$ and the sum caters for the probability that a new node linked to $`i`$ also links to one of the neighbors of $`i`$. This increases $`n_i`$ and enhances clustering. In order to make Eq. (4) more concrete we discuss the triad formation model tunable\_clustering  as an example.

The complications in solving a rate equation like Eq. (4) arise from the correlations that are embedded between the degree of node $`i`$ and the properties of its neighborhood. For the triad formation model, the rules consist of a BA model extended by a triad formation step. Initially, the network contains $`m_0`$ vertices and no edges, and in every time step a new vertex is added with $`m`$ undirected edges. The $`m`$ edges are then one-by-one subsequently linked to $`m`$ different nodes in the network. One performs a preferential attachment step for the first edge as defined in the BA model. With probability $`p`$, the second and further edges edges are joined to a randomly chosen neighbor of the node selected in the previous PA step. Alternatively, with probability $`1p`$, a PA step is performed again.

In the limit when $`p`$ approaches zero, one recovers the original BA model, and by setting $`p`$ to a value between $`0`$ and $`1`$ the average clustering can be adjusted continuously and grows monotonously with an increasing $`p`$. The microscopic mechanisms that increase $`n_i`$ are illustrated in Figure 1 and are (I): the new node connects to node $`i`$ in a PA step, which is potentially followed by several TF steps; (II): the new node connects to one of the neighbors of $`i`$ in a PA step and then $`i`$ conversely gets linked to the new node in one of the subsequent TF steps; (III): the new node connects to node $`i`$ in a PA step and a neighbor of $`i`$ is also selected for connection to the new node in another PA step.

Using these for $`R(k_i,p)`$, the rate equation for $`n_i`$ reads

$`_tn_i`$ $`=`$ $`m_{\mathrm{𝑃𝐴}}{\displaystyle \frac{k_i}{2mt}}m_{\mathrm{𝑇𝐹}}+m_{\mathrm{𝑃𝐴}}{\displaystyle \underset{n\mathrm{\Omega }}{}}{\displaystyle \frac{k_n}{2mt}}{\displaystyle \frac{1}{k_n}}m_{\mathrm{𝑇𝐹}}+`$ (5)
$`+`$ $`m_{\mathrm{𝑃𝐴}}{\displaystyle \frac{k_i}{2mt}}(m_{\mathrm{𝑃𝐴}}1){\displaystyle \underset{n\mathrm{\Omega }}{}}{\displaystyle \frac{k_n}{2mt}}.`$
The first term in the sum gives the increase in $`n_i`$ by mechanism (I). $`m_{\mathrm{𝑃𝐴}}`$ is the number of PA steps attempted per each new node (recall that per time-unit one new node is added). $`k_i/(2mt)`$ is the preferential attachment probability to node $`i`$. $`m_{\mathrm{𝑇𝐹}}`$ is the expected number of triad formation steps that take place on the average after a single PA step. Given this, we have that $`m_{\mathrm{𝑃𝐴}}+m_{\mathrm{𝑃𝐴}}m_{\mathrm{𝑇𝐹}}=m`$.

The second term describes mechanism (II); in this term, the sum goes over all neighbors $`\mathrm{\Omega }`$ of $`i`$, and their degrees are denoted by $`k_n`$. $`1/k_n`$ comes from the fact that the neighboring node where a TF step links is chosen uniformly from the neighbors. We exclude here secondary triangle formation, that takes place if two TF steps from the new node form a triangle with $`i`$ and one of $`i`$’s neighbors. This becomes more relevant for large $`p`$’s. The term for (II) gives the same expression as (I) after simplification.

The last term belongs to (III) and it is the only one that would remain if we considered the simple BA model. It is the product of the probabilities of linking to node $`i`$ and to one of the neighbors of $`i`$, respectively, using only PA steps. The term contains the sum of the degrees of neighboring nodes; this is $`k_i`$ times the average degree of the neighbors. It has been shown that for uncorrelated random BA networks $`k_n=\frac{k}{4}\mathrm{log}t=\frac{m}{2}\mathrm{log}t`$ epidemic\_structured\_sf . In this model the numerical result follows the same scaling not only for $`p1`$ but for $`p`$ general.

Finally, we get $`n_i`$ at the end of the network growth by integrating both sides of Eq. (5). The integral for term (I) or (II) is simply

$`{\displaystyle _1^N}m_{\mathrm{𝑃𝐴}}{\displaystyle \frac{k_i}{2mt}}m_{\mathrm{𝑇𝐹}}𝑑t={\displaystyle \frac{m_{\mathrm{𝑃𝐴}}m_{\mathrm{𝑇𝐹}}}{m}}{\displaystyle _1^N}{\displaystyle \frac{dk_i}{dt}}𝑑t=`$
$`={\displaystyle \frac{m_{\mathrm{𝑃𝐴}}m_{\mathrm{𝑇𝐹}}}{m}}\left[k_i(N)m\right]{\displaystyle \frac{m_{\mathrm{𝑃𝐴}}m_{\mathrm{𝑇𝐹}}}{m}}k_i(N),`$ (6)
where we made use of the fact that $`_tk_i=k_i/(2t)`$ tunable\_clustering . From this, it also follows that $`k_i(t)=m(t/i)^{1/2}`$ and thus integrating (III) gives $`[m_{\mathrm{𝑃𝐴}}(m_{\mathrm{𝑃𝐴}}1)/16m][(\mathrm{log}N)^2/N]k_i^2`$. Combining this with Eq. (6) yields

$`n_i=n_{i,0}+{\displaystyle \frac{2m_{\mathrm{𝑃𝐴}}m_{\mathrm{𝑇𝐹}}}{m}}k_i+{\displaystyle \frac{m_{\mathrm{𝑃𝐴}}(m_{\mathrm{𝑃𝐴}}1)}{16m}}{\displaystyle \frac{(\mathrm{log}N)^2}{N}}k_i^2.`$ (7)
The local clustering coefficient for node $`i`$ becomes

$`C_i(k_i)={\displaystyle \frac{n_i}{k_i(k_i1)/2}}{\displaystyle \frac{4m_{\mathrm{𝑇𝐹}}}{k_i}}+{\displaystyle \frac{m1}{8}}{\displaystyle \frac{(\mathrm{log}N)^2}{N}},`$ (8)
after neglecting $`n_{i,0}`$ and approximating $`m_{\mathrm{𝑃𝐴}}`$ by $`m`$, which is reasonable when the triad formation probability is small. It is not surprising that the constant offset in the expression of $`C_i`$ is for $`p0`$ exactly the constant clustering coefficient of pure BA graphs. The first term, more importantly, can be attributed to the triad formation induced clustering, and shows the $`1/k`$ behavior typical of many real networks and other models hierarchical\_organization ; growing\_scale\_free ; eckmann . $`C_i`$ is composed of a power law and a constant, so perfect power-law behavior follows only when the former one dominates. In the opposite case an effective exponent will be less than $`1`$. Furthermore, since $`n_{i,0}`$ has been neglected, Eq. (8) and the inverse proportionality apply to nodes with $`k_i`$ large enough, only.

For further progress $`m_{\mathrm{𝑇𝐹}}`$, the expected number of links created in the several possible TF steps after a PA step for a particular node, needs to be approximated. Take $`m1`$ edges to be available for successive TF steps (this is an upper limit) and assume node $`i`$ is not saturated yet as far as the connections to the neighbors are concerned. This gives $`m_{\mathrm{𝑇𝐹}}=_{z=1}^{m2}zp^z(1p)+(m1)p^{m1}p`$ for $`p`$ small.

The fact that the local clustering coefficient contains a constant term means that there is a crossover at a certain $`k^{}`$. At this point, a power law turns over to a constant clustering coefficient. $`k^{}`$ can be estimated by taking the two terms in Eq. (8) to be equal:

$`k^{}`$ $``$ $`{\displaystyle \frac{32}{m(\mathrm{log}N)^2}}pN.`$ (9)
Thus we can conclude that the exponents of Eq. (3) are $`\gamma =1/\alpha =1`$ and $`\delta =1`$ for the triad formation model, and from above, $`\alpha =1`$. E.g., in the case of Figure 2, taking $`N=10^6`$ yields $`k^{}400`$.

Simulations of the model consistently confirm the analytical results obtained from the rate equation. In Figure 2 networks of different sizes are shown to undergo such a transition to constant clustering by tuning $`p`$ so that $`k^{}`$ is smaller than the maximum degree in the networks. A similar phenomenon to the transition described above can be observed in the case of the actor network of the IMDB database hierarchical\_organization , where the tail of a decreasing power law becomes constant, although large fluctuations naturally affect this part of the statistics. Figure 3 shows networks well below the transition and thus almost only the power-law part is conceivable.

It is not unusual in the physics of scale free networks that mean-field approaches work well general . This fact is related to the strongly hierarchical nature of the networks grown by preferential attachment and our study demonstrates that this situation remains unaltered even when considering a mechanism which enhances clustering. The agreement between the $`1/k^\alpha `$ dependence with $`\alpha =1`$ obtained in Eq. (8) and that found in real networks indicates that the same “mean-field” mechanisms of clustering are operative. For PA growth with enhanced clustering the simplest interpretation is that for each new link a node $`i`$ gains from a new node introduced to the network, its neigbors (“friends”) have also a constant probability to be linked to the same new one.

It is interesting to ask how robust the mean-field exponent is and what are the limits of the above approach, especially in the light of the recently discovered networks with $`\alpha 1`$ CCR . The rate equations allow to discuss the ways how exponents like such can emerge. Eq. (4) implies that the clustering is crucially dependent on the properties of the nodes in the neighborhood, $`\mathrm{\Omega }`$. If, say, correlations from “assortative” or “disassortative” mixing arise between $`k_i`$ and the average degree $`k_n`$ ($`n\mathrm{\Omega }`$assortative\_mixing , this may either enhance ($`\alpha <1`$) or inhibit ($`\alpha >1`$) clustering from the mean-field result. On the level of models, one can envision changing the $`k`$\- and $`p`$-dependence of the rates. The second possibility is fluctuation effects that limit the validity of the rate-equation theory. It would seem to be of interest to explore both these issues.

In conclusion, we have formulated a scaling assumption and a mean-field theory of the clustering of scale-free networks. A specific example, the triad formation model has been solved and comparisons to the simulations indicate both good agreement and yield the MF-value of the exponent $`\alpha `$. This approach should be amenable to many of the models in the literature, and help in understanding the origins of the statistical properties of clustering, also beyond the $`C(k)`$-function, both in models and in the many real-life examples of networks. We have here considered only growing networks, but obviously the rate equations can be written down also in the case the structural dynamics allows for deleting edges, as well CCR2 ; delete .

Acknowledgements: JK and GS wish to thank for the warm hospitality at HUT and for partial support by OTKA T029985. This work has been supported by the Academy of Finland’s Centre of Excellence Programme and by the Centre for International Mobility (CIMO).


format: text

text:
# 𝐵⁰-𝐵̄⁰ mixing in quenched lattice QCD

## I Introduction

The determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix element $`|V_{td}|`$ plays a crucial role in testing the unitarity relation of the CKM matrix, since the position of the vertex of the unitarity triangle would be essentially identified together with the angle $`\varphi _1`$ of the unitarity triangle. Now that the angle $`\varphi _1`$ has already been measured experimentally by the asymmetric $`B`$ factories Aubert:2001nu ; Abe:2001xe  and its precision is expected to be improved substantially in the near future, the accuracy of $`|V_{td}|`$ really determines the accuracy of the Standard model predictions. Then, other measurements of the CKM matrix elements, such as the determination of $`|V_{ub}|`$ through a measurement of $`bul\nu `$, may be used for a test of the CKM mechanism of the quark flavor mixing and the CP violation in the Standard Model.

The CKM matrix element $`|V_{td}|`$ may be determined using the mass difference $`\mathrm{\Delta }M_d`$ in the neutral $`B`$ meson mixing, as it emerges through a loop diagram mediated by top quark and $`W`$ boson, which is proportional to $`|V_{td}V_{tb}^{}|^2`$. The precision in the current world average ($`\mathrm{\Delta }M_d`$ = 0.489 $`\pm `$ 0.008 ps-1 Hocker:2001jb ) is already as good as 1.6%. The constraint on $`|V_{td}|`$ is, however, limited by the theoretical uncertainty in the calculation of the hadronic matrix element $`\overline{B}^0|𝒪_{L_d}|B^0`$ of the $`\mathrm{\Delta }B=2`$ four-quark operator $`𝒪_{L_d}=\overline{b}\gamma _\mu (1\gamma _5)d\overline{b}\gamma _\mu (1\gamma _5)d`$. It is usually parameterized as $`\frac{8}{3}f_B^2B_BM_B^2`$ using the $`B`$ meson leptonic decay constant $`f_B`$ and the $`B`$ parameter $`B_B`$. In the vacuum saturation approximation, which is valid when both $`b`$ and (anti-)$`d`$ quarks are nonrelativistic, $`B_B`$ is normalized to unity.

The best available theoretical method to calculate $`f_B`$ and $`B_B`$ is the numerical simulation of QCD on the lattice, whose current status is reviewed in Ryan:2001ej ; Kronfeld:2001ss ; Yamada:2002\_rev . For the decay constant $`f_B`$, several groups investigated the systematic errors in the lattice calculation, performing the simulations on several different lattices. It is found that the error associated with the large $`b`$ quark mass is controlled reasonably well if one uses an effective theory for heavy quark such as the non-relativistic QCD (NRQCD) Thacker:1991bm ; Lepage:1992tx  or the Fermilab formalism El-Khadra:1997mp , and the results are insensitive to the lattice spacing Aoki:1998ji ; El-Khadra:1998hq ; Bernard:1998xi ; Ishikawa:2000xu ; Collins:2001ix .

For the $`B`$ parameter $`B_B`$, on the other hand, most lattice calculations relying on the effective theory for heavy quark are limited to the static approximation, in which the $`1/m_b`$ correction is neglected Ewing:1996ih ; Christensen:1997sj ; Gimenez:1999mw , and the study of the systematic uncertainty depending on the lattice spacing has not been made. Recently some of us used the NRQCD action, for the first time, to calculate $`B_B`$ Hashimoto:1999ck , and the $`1/m_b`$ correction was studied at a fixed lattice spacing Hashimoto:2000eh . They also calculated another $`B`$ parameter $`B_S`$ Hashimoto:2000eh ; Hashimoto:2000yh  to parameterize the matrix element of the operator $`𝒪_{S_s}=\overline{b}(1\gamma _5)s\overline{b}(1\gamma _5)s`$, which appears in the heavy quark expansion of the width difference of $`B_s`$ Beneke:1996gn ; Beneke:1999sy ; Beneke:2001cu .

In this paper we extend the previous studies Hashimoto:1999ck ; Hashimoto:2000yh ; Hashimoto:2000eh  to investigate the systematic errors in the calculation of the $`B`$ meson $`B`$ parameters. Using the same NRQCD action as in Hashimoto:1999ck ; Hashimoto:2000eh ; Hashimoto:2000yh  we calculate the $`B`$ parameters at four lattice spacings to estimate the size of systematic errors depending on the lattice spacing. In order to minimize other sources of systematic errors, we use the $`O(a)`$-improved Wilson quark action Sheikholeslami:1985ij  for light quarks with the improvement coefficient $`c_{\mathrm{SW}}`$ calculated at the one-loop level Wohlert:1987rf ; Naik:1993ux ; Luscher:1996vw  and nonperturbatively Luscher:1997ug .

Since NRQCD is an effective theory valid for heavy quark and the action is constructed by an expansion in inverse heavy quark mass, there is a potential source of systematic error due to the truncation of $`1/m_b`$ expansion. Furthermore, in order to match the effective theory to the full theory one has to use perturbation theory, and errors from higher order corrections should also be taken into account. We introduce a method to estimate these systematic errors by treating the neglected higher order terms in different ways. It turned out that the error estimated in this way is quite consistent with a naive order counting assuming typical sizes for the expansion parameters.

Due to the systematic errors discussed above, it is not straightforward to obtain an accuracy better than 10–15% for $`f_B\sqrt{B_B}`$ which is relevant for the determination of $`|V_{td}|`$. Alternatively, one could use the ratio $`\mathrm{\Delta }M_s/\mathrm{\Delta }M_d`$, once the mass difference in the $`B_s\overline{B}_s`$ mixing is measured. (The current experimental bound is $`\mathrm{\Delta }M_s>13.1`$ ps-1 at 95% CL PDG2002 .) It is proportional to $`\xi ^2|V_{ts}/V_{td}|^2`$, where $`\xi `$ is a ratio to describe SU(3) flavor breaking of the hadronic matrix elements given by

$$\xi =\frac{f_{B_s}\sqrt{B_{B_s}}}{f_B\sqrt{B_B}}.$$
(1)
Since the bulk of the systematic errors in the calculations of $`f_B`$ and $`B_B`$ cancels in this ratio, one may achieve much better accuracy, as stressed in Bernard:1998dg . The largest remaining uncertainty comes from the chiral extrapolation of lattice data, which is also discussed in this paper.

The $`B`$ meson $`B`$ parameters have also been calculated using the conventional relativistic actions for heavy quark Becirevic:2001nv ; Lellouch:2001tw . Since the lattice spacing in the present simulations is not small enough compared to the Compton wave length of the $`b`$ quark, one has to extrapolate the results obtained around the charm quark mass to the bottom quark mass, which is a significant source of systematic uncertainty. In fact, the extrapolation with the linear form in $`1/M`$ does not seem to agree with explicit calculations in the static limit Hashimoto:2000bk ; Bernard:2001ki ; Gimenez:2001jj . Therefore one may use the static result to constrain the heavy quark extrapolation in the infinite mass limit Becirevic:2001xt . We present a comparison of our result with these previous calculations.

This paper is organized as follows. In the next section, we summarize some phenomenological formulas for the mass and width difference in the $`B^0\overline{B}^0`$ mixing. The lattice action and operators we employ in this work are defined in Section III, where the method to extract the continuum $`B`$ parameters from lattice matrix elements is also described. Simulation details and results are given in Section IV and V, respectively. In Section VI we present our results for the $`B`$ parameters and their systematic uncertainties are discussed. Using these results we also predict the mass and width difference of $`B_s`$ meson system. Calculation of the SU(3) breaking ratio $`\xi `$ is briefly discussed in Section VII, and our conclusions are given in Section VIII.

Preliminary reports of this work have already been presented in Yamada:2001ym ; Hashimoto:2001zq ; Yamada:2001ik ; Yamada:2001xp .

## II Phenomenological formulas

In this section we summarize the phenomenological formulas which involves the $`B`$ meson $`B`$ parameters. We also present some notations which will be used throughout the paper.

### II.1 Mass difference

In the Standard Model the mass difference in the neutral $`B_q^0\overline{B}_q^0`$ mesons ($`q`$ denotes $`d`$ or $`s`$) is given by

$$\mathrm{\Delta }M_q=|V_{tb}^{}V_{tq}|^2\frac{G_F^2m_W^2}{16\pi ^2M_{B_q}}S_0(x_t)\eta _{2B}\left[\alpha _s(\mu _b)\right]^{\gamma _0/2\beta _0}\left[1+\frac{\alpha _s(\mu _b)}{4\pi }J_5\right]\overline{B}_q^0|𝒪_{L_q}(\mu _b)|B_q^0.$$
(2)
$`S_0(x_t)`$ ($`x_t=m_t^2/m_W^2`$) is the Inami-Lim function Inami:1981fz  and $`\eta _{2B}`$ is the short distance QCD correction Buras:1990fn , whose full expression is found in Buchalla:1996vs .

The four-quark operator $`𝒪_{L_q}(\mu _b)`$ is defined as

$$𝒪_{L_q}(\mu _b)=\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1\gamma _5)q,$$
(3)
which depends on the renormalization scale $`\mu _b`$ if it is defined in the continuum renormalization schemes, such as the naive dimensional regularization (NDR) with the modified minimal subtraction ($`\overline{\mathrm{MS}}`$) scheme. The scale dependence cancels with the prefactor $`\left[\alpha _s(\mu _b)\right]^{\gamma _0/2\beta _0}\left[1+\frac{\alpha _s(\mu _b)}{4\pi }J_5\right]`$ such that the physical mass difference is scale independent. In the NDR-$`\overline{\mathrm{MS}}`$ scheme the anomalous dimensions are written as

$$J_{n_f}=\frac{\gamma _0\beta _1}{2\beta _0^2}\frac{\gamma _1}{2\beta _0},$$
(4)
and

$$\begin{array}{cc}\beta _0=11\frac{2}{3}n_f,& \beta _1=102\frac{38}{3}n_f,\\ \gamma _0=4,& \gamma _1=7+\frac{4}{9}n_f.\end{array}$$
(5)
The renormalization scale $`\mu _b`$ is usually taken at the $`b`$ quark mass $`m_b`$.

The $`B`$ parameter $`B_{B_q}`$ is defined through

$$\overline{B}_q^0|𝒪_{L_q}(\mu _b)|B_q^0=\frac{8}{3}B_{B_q}(\mu _b)f_{B_q}^2M_{B_q}^2,$$
(6)
and the scale-independent $`\widehat{B}_{B_q}`$ is given by

$$\widehat{B}_{B_q}=\left[\alpha _s(\mu _b)\right]^{\gamma _0/2\beta _0}\left[1+\frac{\alpha _s(\mu _b)}{4\pi }J_5\right]B_{B_q}(\mu _b).$$
(7)
The number of flavors $`n_f`$ is 5. To evaluate this expression we use the strong coupling constant $`\alpha _s(\mu _b)`$ defined in the $`\overline{\mathrm{MS}}`$ scheme with $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{(5)}`$ = 225 MeV, which corresponds to $`\alpha _s^{(5)}(4.8\text{GeV})`$ = 0.216.

The bulk of the theoretical uncertainties cancels in the ratio

$$\frac{\mathrm{\Delta }M_s}{\mathrm{\Delta }M_d}=\left|\frac{V_{ts}}{V_{td}}\right|^2\frac{M_{B_s}}{M_{B_d}}\xi ^2,$$
(8)
where $`\xi `$ describes the SU(3) flavor breaking of the matrix element $`\overline{B}_q^0|𝒪_{L_q}(\mu _b)|B_q^0`$ as defined in (1). If one assumes the unitarity relation among the CKM matrix elements $`|V_{ts}||V_{cb}|`$, (8) may be used to determine $`|V_{td}|`$.

### II.2 Width difference

Using the heavy quark expansion, the width difference in the neutral $`B_s^0\overline{B}_s^0`$ mixing is calculated as Beneke:1996gn ; Beneke:1999sy

$$\mathrm{\Delta }\mathrm{\Gamma }_{B_s}=2\frac{1}{2M_{B_s}}\overline{B}_q|\mathrm{Im}id^4xT_{\mathrm{eff}}(x)_{\mathrm{eff}}(0)|B_q,$$
(9)
where $`_{\mathrm{eff}}`$ is the $`\mathrm{\Delta }B`$ = 1 weak transition Hamiltonian. The main contribution comes from the transition $`b\overline{s}c\overline{c}`$ followed by $`c\overline{c}\overline{b}s`$, and other contributions mediated by penguin operators are also considered.

The operator product expansion (OPE) may be used to approximate the transition operator $`id^4xT_{\mathrm{eff}}(x)_{\mathrm{eff}}(0)`$, which gives an $`1/m_b`$ expansion. At the leading order in $`1/m_b`$ the $`\mathrm{\Delta }B`$ = 2 four-quark operator $`𝒪_{L_s}`$ defined in (3) and another operator

$$𝒪_{S_s}=\overline{b}(1\gamma _5)s\overline{b}(1\gamma _5)s,$$
(10)
appear. Then, the following formula

$`\left({\displaystyle \frac{\mathrm{\Delta }\mathrm{\Gamma }}{\mathrm{\Gamma }}}\right)_{B_s}`$ $`=`$ $`{\displaystyle \frac{16\pi ^2B(B_sXe\nu )}{g(z)\stackrel{~}{\eta }_{QCD}}}{\displaystyle \frac{f_{B_s}^2M_{B_s}}{m_b^3}}|V_{cs}|^2`$ (11)
$`\times \left(G(z){\displaystyle \frac{8}{3}}B_{B_s}(m_b)+G_S(z){\displaystyle \frac{5}{3}}{\displaystyle \frac{B_{S_s}(m_b)}{(m_b)^2}}+\sqrt{14z}\delta _{1/m}\right)`$
is obtained at the next-to-leading order Beneke:1999sy , where $`m_b`$ is the pole mass of $`b`$ quark. The width difference is normalized by the total decay width of $`B_s`$ meson $`\mathrm{\Gamma }_{B_s}`$, which is written in terms of the semileptonic decay branching ratio $`B(B_sXe\nu )`$ on the right hand side in order to remove an uncertainty in the value of $`|V_{cb}|`$. The phase space factor $`g(z)=18z+8z^3z^412z^2\mathrm{ln}z`$, where $`z=m_c^2/m_b^2`$, and the short distance QCD correction111 Following the treatment in Beneke:1999sy , the approximate form of Kim:1989ac  is used for $`\stackrel{~}{\eta }_{QCD}`$.

$$\stackrel{~}{\eta }_{QCD}=1\frac{2\alpha _s(m_b)}{3\pi }\left[(\pi ^2\frac{31}{4})(1\sqrt{z})^2+\frac{3}{2}\right]$$
(12)
are known factors, and the functions $`G(z)`$ and $`G_S(z)`$ describe the next-to-leading order QCD corrections Beneke:1999sy  appearing in the calculation of the operator product expansion. The correction term $`\delta _{1/m}`$ denotes the next-to-leading order contribution in the $`1/m_b`$ expansion, which is estimated in Beneke:1996gn  using the factorization approximation.

The $`B`$ parameter $`B_{S_s}`$ is defined through

$$\overline{B}_s^0|𝒪_{S_s}(\mu _b)|B_s^0=\frac{5}{3}f_{B_s}^2\frac{B_{S_s}(\mu _b)}{(\mu _b)^2}M_{B_s}^2,$$
(13)
where

$$(\mu _b)=\frac{\overline{m}_b(\mu _b)+\overline{m}_s(\mu _b)}{M_{B_s}}$$
(14)
is the ratio of matrix elements of heavy-light axial vector current and pseudoscalar density and $`\overline{m}(\mu _b)`$ represents a $`\overline{\mathrm{MS}}`$ quark mass.

In the following analysis, the scale $`\mu _b`$ is set to the pole mass of $`b`$ quark, $`m_b`$ = 4.8 GeV, for which $`G(z)`$=0.03 and $`G_S(z)`$=0.937. With input parameters $`z=0.085`$, $`|V_{cs}|`$ = $`1\frac{\lambda ^2}{2}`$=0.976, $`M_{B_s}`$ = 5.37 GeV, $`B(B_sXe\nu )`$ = 0.107, we obtain

$$\left(\frac{\mathrm{\Delta }\mathrm{\Gamma }}{\mathrm{\Gamma }}\right)_{B_s}=\left(\frac{f_{B_s}}{230\mathrm{MeV}}\right)^2\left[0.007B_{B_s}(m_b)+0.207B_{S_s}(m_b)0.077\right].$$
(15)
For the central value of the decay constant, we choose a recent world average of unquenched lattice calculations $`f_{B_s}`$ = 230(30) MeV Ryan:2001ej ; Yamada:2002\_rev . The uncertainties associated with these input parameters are discussed in Section VI.

## III NRQCD action and operators

In this section we describe the lattice NRQCD action and operators used in our calculations. The perturbative matching of the lattice operators to the continuum ones is summarized.

### III.1 NRQCD action

To treat the heavy quark on a lattice with moderate lattice spacing $`a`$, the idea of the heavy quark effective theory (HQET) Eichten:1990zv ; Georgi:1990um ; Grinstein:1990mj  is useful, as it allows us to describe the heavy quark of mass $`M`$ without introducing large systematic errors scaling as a (positive) power of $`aM`$. In HQET the Lagrangian is organized as an expansion in inverse powers of $`M`$, and the terms beyond some fixed order are truncated. Since the physical expansion parameter is $`\mathrm{\Lambda }_{\mathrm{QCD}}/M`$ with $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ 300–500 MeV, one may typically achieve the accuracy of order few percent for $`B`$ mesons at the next-to-leading order, i.e. including terms of order $`1/M`$.

Though NRQCD was originally introduced in the continuum Caswell:1986ui ; Bodwin:1995jh  and on the lattice Thacker:1991bm ; Lepage:1992tx  to describe the quarkonium systems such as charmonium and bottomonium, for which the expansion parameter is a velocity of heavy quark rather than $`\mathrm{\Lambda }_{\mathrm{QCD}}/M`$, the formulation can also be used for the lattice study of heavy-light mesons as first demonstrated in Hashimoto:1994nd . At the next-to-leading order in $`\mathrm{\Lambda }_{\mathrm{QCD}}/M`$, the Lagrangian in the continuum Euclidean space-time is written as

$$_{\mathrm{NRQCD}}^{cont}=Q^{}\left[D_0+\frac{𝐃^2}{2M}+g\frac{\sigma 𝐁}{2M}\right]Q+\chi ^{}\left[D_0\frac{𝐃^2}{2M}g\frac{\sigma 𝐁}{2M}\right]\chi ,$$
(16)
for heavy quark field $`Q`$ and heavy anti-quark field $`\chi `$. Both are represented by a two-component non-relativistic spinor. The derivatives $`D_0`$ and $`𝐃`$ are temporal and spatial covariant derivatives respectively. The leading order term $`D_0`$ represents a heavy quark as a static color source. The leading correction terms of order $`\mathrm{\Lambda }_{\mathrm{QCD}}/M`$ are the non-relativistic kinetic term $`𝐃^2/2M`$ and the spin-(chromo)magnetic interaction term $`\sigma 𝐁/2M`$, where $`𝐁`$ denotes the chromomagnetic field strength. In the usual HQET approach, only the leading terms are kept in the effective Lagrangian and corrections of order $`\mathrm{\Lambda }_{\mathrm{QCD}}/M`$ are treated as operator insertions. Alternatively, in our lattice calculation we include the correction terms in the Lagrangian (16) and evaluate the matrix elements with the heavy quark propagator including the effect of order $`\mathrm{\Lambda }_{\mathrm{QCD}}/M`$.

On the lattice we use a discretized version of the Lagrangian (16), whose explicit form is written as

$$S_{\mathrm{NRQCD}}=\underset{x,y}{}Q^{}(x)(\delta _{x,y}K_Q(x,y))Q(y)+\underset{x,y}{}\chi ^{}(x)(\delta _{x,y}K_\chi (x.y))\chi (y).$$
(17)
The kernel to describe the time evolution of (anti-)heavy quark is given by

$`K_Q(x,y)`$ $``$ $`\left[\left(1{\displaystyle \frac{aH_0}{2n}}\right)^n\left(1{\displaystyle \frac{a\delta H}{2}}\right)\delta _4^{()}U_4^{}\left(1{\displaystyle \frac{a\delta H}{2}}\right)\left(1{\displaystyle \frac{aH_0}{2n}}\right)^n\right](x,y),`$ (18)
$`K_\chi (x,y)`$ $``$ $`\left[\left(1{\displaystyle \frac{aH_0}{2n}}\right)^n\left(1{\displaystyle \frac{a\delta H}{2}}\right)\delta _4^{(+)}U_4\left(1{\displaystyle \frac{a\delta H}{2}}\right)\left(1{\displaystyle \frac{aH_0}{2n}}\right)^n\right](x,y),`$ (19)
where $`n`$ denotes a stabilization parameter introduced in order to remove the instability arising from unphysical momentum modes in the evolution equation Thacker:1991bm ; Lepage:1992tx . The operator $`\delta _4^{(\pm )}`$ is defined as $`\delta _4^{(\pm )}(x,y)\delta _{x_4\pm 1,y_4}\delta _{𝐱,𝐲}`$, and $`H_0`$ and $`\delta H`$ are lattice Hamiltonians defined by

$`H_0`$ $``$ $`{\displaystyle \frac{𝚫^{(2)}}{2aM_0}},`$ (20)
$`\delta H`$ $``$ $`c_B{\displaystyle \frac{g}{2aM_0}}\sigma 𝐁,`$ (21)
where $`𝚫^{(2)}_{i=1}^3\mathrm{\Delta }_i^{(2)}`$ is the Laplacian defined on the lattice with $`\mathrm{\Delta }_i^{(2)}\mathrm{\Delta }_i^{(+)}\mathrm{\Delta }_i^{()}`$, $`\mathrm{\Delta }_i^{(+)}`$ and $`\mathrm{\Delta }_i^{()}`$ being forward and backward covariant derivatives in the $`i`$-th direction. In (21) the chromomagnetic field operator $`𝐁`$ is the usual clover-leaf type lattice field strength Lepage:1992tx . In these definitions, the lattice operators $`\mathrm{\Delta }^{(2)}`$ and $`𝐁`$ are dimensionless, i.e. appropriate powers of $`a`$ are understood. The space-time indices $`x`$ and $`y`$ are implicit in these expressions. The bare heavy quark mass $`M_0`$ is distinguished from the renormalized one $`M`$.

At the tree level, the lattice action (17) describes the continuum NRQCD (16) in the limit of vanishing lattice spacing $`a`$. (We may identify $`M_0`$ with $`M`$ and take the tree level value $`c_B=1`$.) The leading discretization error for the spatial derivative is of order $`(a\mathrm{\Lambda }_{\mathrm{QCD}})^2\mathrm{\Lambda }_{\mathrm{QCD}}/M`$. Since the temporal derivative is discretized asymmetrically, the leading error appears at order $`aD_0^2`$, whose typical size is estimated as $`a\mathrm{\Lambda }_{\mathrm{QCD}}^3/M^2`$ using the equation of motion. The gauge potential part is automatically improved, as it is exponentiated into the temporal link variable $`U_4`$.

In the presence of radiative corrections, the heavy quark mass $`M_0`$ and the chromomagnetic coupling $`c_B`$ have to be tuned in such a way that the continuum values are reproduced at each value of the strong coupling constant $`\alpha _s`$. Furthermore, the radiative corrections generate many other terms which do not exist in the continuum Lagrangian (16), because NRQCD is not a renormalizable field theory. In general these terms appear with some factor of form $`\alpha _s^k/(aM_0)^m`$ with positive integers $`k`$ and $`m`$ ($`k1`$ and $`m2`$).222There are also the lowest dimension operators $`Q^{}Q`$ and $`\chi ^{}\chi `$, but they only give the energy shift and do not contribute to the dynamics of heavy quark.  Therefore, NRQCD should be considered as an effective theory valid for small $`1/(aM)`$ up to higher order terms in $`1/(aM)`$.

Perturbation theory can be used to calculate the renormalization of the parameters. For example, the one-loop calculations of energy shift and mass renormalization were carried out for lattice NRQCD by Davies and Thacker Davies:1992py  and by Morningstar Morningstar:1993de , and by ourselves Hashimoto:1999ck ; Ishikawa:2000xu ; Hashimoto:2000eh  for the above particular form of the NRQCD action.333We note that the evolution kernels (18) and (19) are slightly different from the definition used, for example, in Morningstar:1993de , where the $`(1aH_0/2n)^n`$ terms appear inside of the $`(1a\delta H/2)`$ terms.  To improve the perturbative expansion we utilize the tadpole improvement Lepage:1993xa , namely, all the gauge links in the action (17) are divided by its mean field value $`u_0`$ determined from the plaquette expectation value $`u_0(\mathrm{Tr}U_P/3)^{1/4}`$, where some counter terms are introduced in perturbative calculations. The one-loop tuning of the coupling constant $`c_B`$ for the spin-(chromo)magnetic interaction term (21) has not yet been performed, so we take the tree level value $`c_B`$ = 1 after making the tadpole improvement.

The relativistic four-component Dirac spinor field $`b`$ is related to the two-component non-relativistic field $`Q`$ and $`\chi `$ appearing in the NRQCD action (17) via the Foldy-Wouthuysen-Tani (FWT) transformation

$$b=\left(1\frac{\gamma 𝚫^{(\pm )}}{2aM_0}\right)h,$$
(22)
where

$$h\left(\begin{array}{c}Q\\ \chi ^{}\end{array}\right).$$
(23)
The symbol $`𝚫^{(\pm )}`$ denotes a symmetric covariant differentiation operator $`\mathrm{\Delta }_i^{(\pm )}\mathrm{\Delta }_i^{(+)}\mathrm{\Delta }_i^{()}`$.

### III.2 Bilinear operators

The heavy-light axial-vector current $`𝒜_\mu =\overline{b}\gamma _\mu \gamma _5q`$ and pseudoscalar density $`𝒫=\overline{b}\gamma _5q`$ appear in the definition of the $`B`$ parameters through the vacuum saturation approximation. We use the calligraphic symbols $`𝒜_\mu `$ and $`𝒫`$ to denote the currents defined in the continuum full theory. Since the pseudoscalar density diverges in the continuum, it is renormalized with the $`\overline{\mathrm{MS}}`$ scheme at a scale $`\mu `$. On the other hand, the axial-vector current does not need renormalization, because it is partially conserved in the continuum full theory.

The corresponding lattice operators are

$`J_\mathrm{\Gamma }^{(0)}`$ $`=`$ $`\overline{b}\mathrm{\Gamma }q,`$ (24)
$`J_\mathrm{\Gamma }^{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{2aM_0}}\overline{b}\left(\gamma \stackrel{}{𝚫}^{(\pm )}\right)\mathrm{\Gamma }q,`$ (25)
$`J_\mathrm{\Gamma }^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2aM_0}}\overline{b}\mathrm{\Gamma }\left(\gamma 𝚫^{(\pm )}\right)q,`$ (26)
where $`J_\mathrm{\Gamma }`$ is $`A_4`$ for $`\mathrm{\Gamma }=\gamma _4\gamma _5`$ or $`P`$ for $`\mathrm{\Gamma }=\gamma _5`$. The light quark field $`q`$ is described by the $`O(a)`$-improved Wilson quark action Sheikholeslami:1985ij . We apply the tadpole improvement Lepage:1993xa  for the light quark field using the critical hopping parameter $`\kappa _c`$ to define the mean link variable $`u_0=1/8\kappa _c`$, so that we normalize the light quark field with a factor $`\sqrt{1\frac{3\kappa }{4\kappa _c}}`$. The heavy quark field $`b`$ is defined in (22).

The one-loop matching between the continuum and lattice operators is written as

$$𝒥_\mathrm{\Gamma }=\left[1+\frac{\alpha _s}{4\pi }\rho _\mathrm{\Gamma }^{(0)}\right]J_\mathrm{\Gamma }^{(0)}+\frac{\alpha _s}{4\pi }\rho _\mathrm{\Gamma }^{(1)}J_\mathrm{\Gamma }^{(1)}+\frac{\alpha _s}{4\pi }\rho _\mathrm{\Gamma }^{(2)}J_\mathrm{\Gamma }^{(2)},$$
(27)
with one-loop coefficients $`\rho _\mathrm{\Gamma }^{(i)}`$. The coefficient $`\rho _\mathrm{\Gamma }^{(0)}`$ is written as

$`\rho _A^{(0)}`$ $`=`$ $`2\mathrm{ln}(a^2M^2)+\zeta _A,`$ (28)
$`\rho _P^{(0)}`$ $`=`$ $`{\displaystyle \frac{9}{2}}\mathrm{ln}(\mu ^2/M^2)+{\displaystyle \frac{3}{2}}\mathrm{ln}(a^2M^2)+\zeta _P,`$ (29)
for $`\mathrm{\Gamma }=\gamma _4\gamma _5`$ and $`\gamma _5`$, respectively. In the static limit the numerical constants are $`\zeta _A=16.55`$ Ishikawa:1999rv ; Golden:1991dx ; Borrelli:1992fy  and $`\zeta _P=11.21`$ Ishikawa:1999rv . For the NRQCD action (17) with a finite heavy quark mass $`M_0`$ the numerical values for $`\zeta _A`$ and $`\zeta _P`$ are available in Table III of Hashimoto:2000eh .444The same quantity was previously calculated in Davies:1993ec ; Morningstar:1998ep , but for a slightly different NRQCD action.

In the static limit, while the second term of (27) vanishes, the third term remains finite and describes the $`O(\alpha _sa\mathrm{\Lambda }_{\mathrm{QCD}})`$ improvement, and its coefficient $`\rho _\mathrm{\Gamma }^{(2)}`$ is $`\rho _A^{(2)}/2aM_0=\rho _P^{(2)}/2aM_0=13.01`$ Ishikawa:1999rv . Away from the static limit, these terms give contributions of the $`O(\alpha _sa\mathrm{\Lambda }_{\mathrm{QCD}})`$ and $`O(\alpha _s\mathrm{\Lambda }_{\mathrm{QCD}}/M)`$, and the one-loop coefficients are calculated only for the axial vector current $`\mathrm{\Gamma }=\gamma _4\gamma _5`$ for our choice of the NRQCD action Ishikawa:2000xu .555Note that a different notation is used in Ishikawa:2000xu . Similar calculation was previously made by Morningstar and Shigemitsu Morningstar:1998yx .

### III.3 $`\mathrm{\Delta }B`$=2 operators

We assume that the continuum four-quark operators $`𝒪_L(\mu )`$ and $`𝒪_S(\mu )`$ are renormalized in the $`\overline{\mathrm{MS}}`$ scheme with totally anti-commuting $`\gamma _5`$. In the renormalization of $`𝒪_S(\mu )`$, the subtraction of evanescent operators is made with the definition given by Eqs. (13)-(15) of Beneke:1999sy .

For matching of these four-quark operators, the following lattice operators are involved at the lowest dimension:

$`O_L`$ $`=`$ $`\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1\gamma _5)q,`$ (30)
$`O_R`$ $`=`$ $`\overline{b}\gamma _\mu (1+\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)q,`$ (31)
$`O_S`$ $`=`$ $`\overline{b}(1\gamma _5)q\overline{b}(1\gamma _5)q,`$ (32)
$`O_N`$ $`=`$ $`2\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)q+4\overline{b}(1\gamma _5)q\overline{b}(1+\gamma _5)q,`$ (33)
$`O_M`$ $`=`$ $`2\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)q4\overline{b}(1\gamma _5)q\overline{b}(1+\gamma _5)q,`$ (34)
$`O_P`$ $`=`$ $`2\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)q+12\overline{b}(1\gamma _5)q\overline{b}(1+\gamma _5)q,`$ (35)
$`O_T`$ $`=`$ $`5\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)q34\overline{b}(1\gamma _5)q\overline{b}(1+\gamma _5)q.`$ (36)
As in the bilinear operator case, one has to introduce higher dimensional operators to remove $`O(a)`$ errors at the one-loop level. For the four-quark operators, however, the necessary one-loop calculations to remove the $`O(\alpha _sa\mathrm{\Lambda }_{\mathrm{QCD}})`$ error are made only in the static limit Ishikawa:1999rv ; Hashimoto:2000yh . The higher dimensional operators appearing in that limit are

$`O_{LD}`$ $`=`$ $`\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1\gamma _5)(a\gamma 𝚫^{(\pm )})q,`$ (37)
$`O_{ND}`$ $`=`$ $`2\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)(a\gamma 𝚫^{(\pm )})q`$ (38)
$`+4\overline{b}(1\gamma _5)q\overline{b}(1+\gamma _5)(a\gamma 𝚫^{(\pm )})q,`$
$`O_{SD}`$ $`=`$ $`\overline{b}(1\gamma _5)q\overline{b}(1\gamma _5)(a\gamma 𝚫^{(\pm )})q,`$ (39)
$`O_{PD}`$ $`=`$ $`2\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)(a\gamma 𝚫^{(\pm )})q`$ (40)
$`+12\overline{b}(1\gamma _5)q\overline{b}(1+\gamma _5)(a\gamma 𝚫^{(\pm )})q.`$
The one-loop matching is written as follows:

$`𝒪_L(\mu )`$ $`=`$ $`O_L+{\displaystyle \frac{\alpha _s}{4\pi }}\rho _{L,L}O_L+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{L,S}O_S`$ (41)
$`+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{L,R}O_R+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{L,N}O_N+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{L,M}O_M`$
$`+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{L,LD}O_{LD}+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{L,ND}O_{ND},`$
$`𝒪_S(\mu )`$ $`=`$ $`O_S+{\displaystyle \frac{\alpha _s}{4\pi }}\rho _{S,S}O_S+{\displaystyle \frac{\alpha _s}{4\pi }}\rho _{S,L}O_L`$ (42)
$`+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{S,R}O_R+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{S,P}O_P+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{S,T}O_T`$
$`+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{S,SD}O_{SD}+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{S,LD}O_{LD}+{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{S,PD}O_{PD},`$
where the coefficients $`\rho _{L,L}`$, $`\rho _{S,S}`$, and $`\rho _{S,L}`$ contain the physical scale $`\mu `$ and $`m_b`$ as follows.

$`\rho _{L,L}`$ $`=`$ $`2\mathrm{ln}\left({\displaystyle \frac{\mu ^2}{M^2}}\right)+4\mathrm{ln}(a^2M^2)+\zeta _{L,L},`$ (43)
$`\rho _{S,S}`$ $`=`$ $`{\displaystyle \frac{16}{3}}\mathrm{ln}\left({\displaystyle \frac{\mu ^2}{M^2}}\right)+{\displaystyle \frac{4}{3}}\mathrm{ln}(a^2M^2)+\zeta _{S,S},`$ (44)
$`\rho _{S,L}`$ $`=`$ $`{\displaystyle \frac{1}{3}}\mathrm{ln}\left({\displaystyle \frac{\mu ^2}{M^2}}\right){\displaystyle \frac{2}{3}}\mathrm{ln}(a^2M^2)+\zeta _{S,L}.`$ (45)
The numerical results for the one-loop coefficients $`\zeta _{L,X}`$ ($`X=L`$, $`S`$, $`R`$, $`N`$, and $`M`$) and $`\zeta _{S,X}`$ ($`X=S`$, $`L`$, $`R`$, $`P`$, and $`T`$) are given in Table VI and VIII of Hashimoto:2000eh . In the static limit, $`\zeta _{L,M}`$ and $`\zeta _{S,T}`$ vanish, and others agree with the previous calculations Flynn:1991qz ; Borrelli:1992fy ; DiPierro:1998ty ; Gimenez:1999mw ; Ishikawa:1999rv .666A numerical error in Borrelli:1992fy  was later corrected in DiPierro:1998ty ; Gimenez:1999mw ; Ishikawa:1999rv .

The last lines in (41) and (42) are added to remove the error of $`O(\alpha _sa\mathrm{\Lambda }_{\mathrm{QCD}})`$, but their coefficients are known only in the static limit. Their values are $`\zeta _{L,LD}=17.20`$, $`\zeta _{L,ND}=9.20`$, $`\zeta _{S,SD}=6.88`$, $`\zeta _{S,LD}=2.58`$ and $`\zeta _{S,PD}=1.15`$ Ishikawa:1999rv ; Hashimoto:2000yh .

### III.4 Truncation of expansions

As in the matching of the NRQCD action discussed in Section III.1, we have to truncate the $`1/M`$ and the perturbative expansions in the matching of the bilinear and four-quark operators. The $`1/M`$ expansion is truncated at $`O(1/M)`$, which is consistent with our choice of the NRQCD action (16), and the perturbative corrections of order $`\alpha _s^2`$ and higher are neglected.

In addition, there are mixed corrections of $`O(\alpha _s\mathrm{\Lambda }_{\mathrm{QCD}}/M)`$. In the matching of the bilinear operators the matching coefficients for the mixed corrections are available and such corrections were actually included in Ishikawa:2000xu  by combining with higher dimensional operators as shown in (27). For the four-quark operators, however, the mixing with higher dimensional operators at the one-loop level has not been calculated yet777Except for the static limit, where the mixing terms describe the correction of order $`\alpha _sa`$ rather than $`\alpha _s/M`$. . Thus, in this paper, the mixed corrections are not considered in both of the bilinear and four-quark operators. This means that, for the bilinear operators, only the first term of (27) is taken, thus the matching becomes multiplicative in this approximation.

At this level of accuracy, it is arbitrary to apply the FWT transformation to a heavy quark field, if the heavy quark field forms an operator appearing in the one-loop corrections. Namely, in (41) and (42), we may replace all the four-quark operators $`O_X`$ except for that in the first term by $`O_X^{}`$, where

$$O_L^{}=\overline{h}\gamma _\mu (1\gamma _5)q\overline{h}\gamma _\mu (1\gamma _5)q,$$
(46)
and so on, and the heavy quark field $`h`$ is not rotated by the FWT transformation (22). Therefore $`O_X^{}`$ differ from $`O_X`$ at $`O(\mathrm{\Lambda }_{\mathrm{QCD}}/M)`$. In the naive order counting analysis both choices are equivalent up to unknown corrections of $`O(\alpha _s\mathrm{\Lambda }_{\mathrm{QCD}}/M)`$.

In the calculation of $`B`$ parameters through the ratios

$`B_B(\mu _b)`$ $`=`$ $`{\displaystyle \frac{\overline{B}|𝒪_L(\mu _b)|B}{\frac{8}{3}\overline{B}|𝒜_\mu |00|𝒜_\mu |B}},`$ (47)
$`B_S(\mu _b)`$ $`=`$ $`{\displaystyle \frac{\overline{B}|𝒪_S(\mu _b)|B}{\frac{5}{3}\overline{B}|𝒫(\mu _b)|00|𝒫(\mu _b)|B}},`$ (48)
the perturbative and $`1/M`$ expansions may be truncated in several different ways. A natural choice to match the ratios (47) and (48) is to write the numerator and denominator as they stand:

$`B_B^{(\mathrm{I})}(\mu _b)`$ $`=`$ $`{\displaystyle \frac{\left[1+\frac{\alpha _s}{4\pi }\rho _{L,L}\right]B_{B,L}^{lat}+_{X=S,R,N,M}\frac{\alpha _s}{4\pi }\zeta _{L,X}B_{B,X}^{lat}}{\left[1+\frac{\alpha _s}{4\pi }\rho _A^{(0)}\right]^2}},`$ (49)
$`B_S^{(\mathrm{I})}(\mu _b)`$ $`=`$ $`{\displaystyle \frac{\left[1+\frac{\alpha _s}{4\pi }\rho _{S,S}\right]B_{S,S}^{lat}+\frac{\alpha _s}{4\pi }\rho _{S,L}B_{S,L}^{lat}+_{X=R,P,T}\frac{\alpha _s}{4\pi }\zeta _{S,X}B_{S,X}^{lat}}{\left[1+\frac{\alpha _s}{4\pi }\rho _P^{(0)}\right]^2}}.`$ (50)
A roman numeral, (I) in this case, as a superscript of $`B_B`$ or $`B_S`$ distinguishes the method to truncate the expansion. $`B_{B,X}^{lat}`$ and $`B_{S,X}^{lat}`$ are $`B`$ parameters defined with the lattice operators as

$`B_{B,X}^{lat}`$ $`=`$ $`{\displaystyle \frac{\overline{B}|O_X|B}{\frac{8}{3}\overline{B}|A_4^{(0)}|00|A_4^{(0)}|B}},`$ (51)
$`B_{S,X}^{lat}`$ $`=`$ $`{\displaystyle \frac{\overline{B}|O_X|B}{\frac{5}{3}\overline{B}|P^{(0)}|00|P^{(0)}|B}},`$ (52)
which are directly measured in the numerical simulation from a ratio of correlation functions as we describe in the next section.

Alternatively, one may linearize the perturbative expansion as

$`B_B^{(\mathrm{II})}(\mu _b)`$ $`=`$ $`\left[1+{\displaystyle \frac{\alpha _s}{4\pi }}(\rho _{L,L}2\rho _A^{(0)})\right]B_{B,L}^{lat}+{\displaystyle \underset{X=S,R,N,M}{}}{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{L,X}B_{B,X}^{lat},`$ (53)
$`B_S^{(\mathrm{II})}(\mu _b)`$ $`=`$ $`\left[1+{\displaystyle \frac{\alpha _s}{4\pi }}(\rho _{S,S}2\rho _P^{(0)})\right]B_{S,S}^{lat}+{\displaystyle \frac{\alpha _s}{4\pi }}\rho _{S,L}B_{S,L}^{lat}+{\displaystyle \underset{X=R,P,T}{}}{\displaystyle \frac{\alpha _s}{4\pi }}\zeta _{S,X}B_{S,X}^{lat}.`$ (54)
Formally they are different from the method I by order $`\alpha _s^2`$, which is not known. We expect, however, that perturbative expansion behaves better for the method II, because the contributions from factorized diagrams to the four-quark operators are the same as those of the corresponding bilinear operators, so that the radiative corrections partly cancel in the combination $`\rho _{L,L}2\rho _A^{(0)}`$ or $`\rho _{S,S}2\rho _P^{(0)}`$.

For each method I or II, we also consider the variation of replacing $`O_X`$ and $`J_\mathrm{\Gamma }^{(0)}`$ by $`O_X^{}`$ and $`J_{}^{}{}_{\mathrm{\Gamma }}{}^{(0)}`$, respectively, as discussed above, and define the methods as I’ and II’, where $`J_{}^{}{}_{\mathrm{\Gamma }}{}^{(0)}`$ are defined similarly to (46). The difference of the method I’ (II’) from I (II) is of order $`\alpha _s\mathrm{\Lambda }_{\mathrm{QCD}}/M`$.

Since the level of accuracy of these four methods is equivalent in the naive order counting argument, they can be used to estimate possible systematic errors due to the truncation of expansions.

## IV Lattice simulations

### IV.1 Simulation sets

We have performed numerical simulations at four $`\beta `$ values. For three of them ($`\beta `$=6.1, 5.9, and 5.7), which we call the simulation set A, the $`O(a)`$-improvement coefficient $`c_{\mathrm{SW}}`$ in the light quark action is determined using the one-loop expression $`c_{\mathrm{SW}}=(1/P^{3/4})[1+0.199\alpha _V(1/a)]`$. The one-loop coefficient is calculated in Wohlert:1987rf ; Naik:1993ux ; Luscher:1996vw , and we apply the tadpole improvement Lepage:1993xa  with the plaquette expectation value to define the mean link variable. For the last simulation ($`\beta `$=6.0), which we denote as the simulation set B, the non-perturbative value is used for $`c_{\mathrm{SW}}`$ Luscher:1997ug . Therefore, as far as the light quark sector is concerned, the discretization error is minimized in the set B, for which the leading error is of $`O(a^2)`$, while the effect of $`O(\alpha _s^2a)`$ is remaining in the set A. For the quantities involving heavy quarks, however, both sets of simulations give the same order of accuracy, since the heavy quark action and operators are not improved at the same level.

Simulation parameters are summarized in Table 1. The parameters for the simulation set A are almost the same as in our previous work for the leptonic decay constant Ishikawa:2000xu , except that the number of statistical ensembles is increased in this work to obtain stable signals for three-point functions. The set B is our new simulation set, which is intended for comparison with our recent unquenched simulations Aoki:2001yi ; Aoki:2002yr , and its $`\beta `$ value, $`\beta `$ = 6.0, is chosen so that the inverse lattice spacing becomes about 2 GeV. In this paper, we present only the quenched results leaving the unquenched calculations for future publications.

For both simulation sets, A and B, we take the standard plaquette gauge action, and the configuration generation and gauge fixing are made as in Ishikawa:2000xu . Four values of the light quark hopping parameter are chosen for each $`\beta `$ as given in Table 1. They correspond to the light quark mass $`m_q`$ covering the range $`m_s/2<m_q<2m_s`$, where $`m_s`$ denotes the physical strange quark mass. The hopping parameter corresponding to the strange quark mass is determined using the $`K`$ or $`\varphi `$ meson masses as input, and will be denoted as $`\kappa _{s1}`$ and $`\kappa _{s2}`$ respectively.

The heavy quark mass in our simulation ranges from $`2m_b/3`$ to $`4m_b`$. The smallest heavy quark mass in the lattice unit $`aM_0`$ is limited around unity due to the possibly large systematic error in the matching calculation as discussed in Section III. The limit in the heaviest side is set by the exponentially growing statistical error Hashimoto:1994nd .

The lattice spacing $`a`$ is determined through the string tension (for the set A) or the rho meson mass (for the set B). For the simulation set A it is confirmed that both determinations are in good agreement (3.5% variation depending on $`\beta `$) Aoki:1998ji . Therefore, in effect the lattice spacing is set using the rho meson mass for both data sets.

We use the simulation set B to obtain our central value and the other to investigate the systematic errors depending on the lattice spacing. The primary reason is that the discretization error is minimized by the non-perturbative $`O(a)`$ improvement for the set B. The set B is also advantageous since we have larger statistics and hence the numerical results are more stable.

### IV.2 Correlation functions

The method to calculate two- and three-point functions mostly follows that of Hashimoto:1999ck . We put a local source at the origin of the lattice and solve for the light quark propagator. The heavy quark and anti-quark propagators are obtained from the same local source by solving the evolution kernels (18) and (19), respectively.

Three point functions are constructed as

$$C_X^{(3)}(t_1,t_2)=\underset{\stackrel{}{x}_1}{}\underset{\stackrel{}{x}_2}{}0|𝒯\left[A_{4}^{S}{}_{}{}^{}(t_1,\stackrel{}{x}_1)O_X(0,\stackrel{}{0})A_{4}^{S}{}_{}{}^{}(t_2,\stackrel{}{x}_2)\right]|0,$$
(55)
where $`O_X`$ is one of the four-quark operators defined in (30)–(36). We take $`t_1>0`$ and $`t_2<0`$ so that a $`\overline{B}`$ meson propagates in the positive direction in time and a $`B`$ meson propagates in the opposite direction. We also measure two point functions

$`C_A^{(2)}(t)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}0|𝒯\left[A_{4}^{S}{}_{}{}^{}(t,\stackrel{}{x})A_4^{(0)}(0,\stackrel{}{0})\right]|0,`$ (56)
$`C_P^{(2)}(t)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}0|𝒯\left[A_{4}^{S}{}_{}{}^{}(t,\stackrel{}{x})P^{(0)}(0,\stackrel{}{0})\right]|0,`$ (57)
for positive and negative values of $`t`$.

A smeared current $`A_4^S`$, defined as

$$A_4^S(t,\stackrel{}{x})=\underset{\stackrel{}{y}}{}\varphi (\stackrel{}{y})\overline{b}(t,\stackrel{}{x}+\stackrel{}{y})\gamma _4\gamma _5q(t,\stackrel{}{x}),$$
(58)
is used to enhance the overlap with the ground state $`B`$ meson. We measure the smearing function $`\varphi (\stackrel{}{x})`$ for each set of heavy and light quark masses with a limited number of gauge configurations before starting the main simulation.

We extract the lattice $`B`$ parameters $`B_{B,X}^{lat}`$ (51) and $`B_{S,X}^{lat}`$ (52) from the following ratios

$`R_{B,X}(t_1,t_2){\displaystyle \frac{C_X^{(3)}(t_1,t_2)}{\frac{8}{3}C_A^{(2)}(t_1)C_A^{(2)}(t_2)}}`$ $`\stackrel{|t_i|1}{}`$ $`B_{B,X}^{lat},`$ (59)
$`R_{S,X}(t_1,t_2){\displaystyle \frac{C_X^{(3)}(t_1,t_2)}{\frac{5}{3}C_P^{(2)}(t_1)C_P^{(2)}(t_2)}}`$ $`\stackrel{|t_i|1}{}`$ $`B_{S,X}^{lat},`$ (60)
for large enough $`|t_i|`$ ($`i=1,2`$).

### IV.3 Meson masses

In order to calculate the heavy-light meson masses precisely, we also calculate two-point functions with the smeared source and local sink,

$`C_A^{(2)LS}(t)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}0|𝒯\left[A_4^{(0)}(t,\stackrel{}{x})A_4^S(0,\stackrel{}{0})\right]|0,`$ (61)
$`C_P^{(2)LS}(t)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}0|𝒯\left[A_4^{(0)}(t,\stackrel{}{x})P^S(0,\stackrel{}{0})\right]|0,`$ (62)
for which the statistical signal is much better than those with the local source and smeared sink. The heavy-light meson mass is, then, obtained by adding the binding energy $`E_{\mathrm{bin}}`$ measured from the two-point functions and bare quark mass $`aM_0`$. Including one-loop corrections we use

$$aM_P=Z_maM_0+E_{\mathrm{bin}}\delta m,$$
(63)
where perturbative corrections $`Z_m`$ and $`\delta m`$ are given as

$`Z_m`$ $`=`$ $`1+\alpha _sB,`$ (64)
$`\delta m`$ $`=`$ $`\alpha _sA,`$ (65)
and $`A`$ and $`B`$ for each bare quark mass are given in Table. I of Ishikawa:2000xu .

## V Simulation results

### V.1 Ratio of correlation functions

We first extract the $`B`$ parameters defined on the lattice $`B_{B,X}^{lat}`$ and $`B_{S,X}^{lat}`$, which are obtained from the asymptotic behavior of the ratios $`R_{B,X}(t_1,t_2)`$ and $`R_{S,X}(t_1,t_2)`$ as (59)–(60). In Figures 14 (top and middle panels in each Figure) we plot these ratios as a function of $`t_1`$ for some fixed values of $`t_2`$. For illustration we show the operators giving leading contributions, i.e. $`B_{B,L}^{lat}`$ and $`B_{S,S}^{lat}`$, for the heavy quark mass closest to the physical $`b`$ quark mass and the lightest quark mass (largest $`\kappa `$ value). We obtain similar plots for other mass parameters, but the statistical signal becomes much noisier for larger heavy quark mass.

The range of $`t`$ ($`t_1`$ and $`t_2`$) to be included in the fit of the ratios has to be chosen such that the effect of excited states is negligible. We identify the plateau seen in the plots of $`R_{B,X}(t_1,t_2)`$ and $`R_{S,X}(t_1,t_2)`$ as the region where the ground state contribution dominates. To be more conservative, we also check that the plateau is reached for the effective mass plot of two-point functions $`C_A^{(2)}(t)`$ and $`C_A^{(2)LS}(t)`$, which are calculated for the same smearing function as used in the calculation of three-point functions. The plots are shown in the bottom panel of Figures 14.

In the fit of the ratios we take a range of $`t`$ as wide as possible in order to avoid possible contamination from the statistical fluctuation Aoki:1996bb . The fit is done for a fixed value of $`t_2=t_{1min}`$ and changing $`t_1`$ in the range $`[t_{1min},t_{1max}]`$. The value of $`[t_{1min},t_{1max}]`$ is listed in Table 1. In order to quantify the possible effect from excited state contamination, we also carried out a fit with larger values of $`t_{1min}`$ (= $`t_2`$). Since the statistical error grows rapidly as $`t_{1min}`$ is taken larger, the maximum change for $`t_{1min}`$ is chosen to keep the statistical error smaller than 8–10%. The effect for $`B_B`$ is found to be 1% or less except for $`\beta `$ = 6.1 where it is at most 3%. For $`B_S`$ it brings a 1.5–3% effect except for $`\beta `$=5.7 where it is negligible. The variations among different choices of the fit range are taken into account in the final results.

### V.2 Chiral extrapolation

The results of $`B_{B,X}^{lat}`$ and $`B_{S,X}^{lat}`$ are insensitive to the light quark mass. An example is shown in Figure 5, where the data at $`\beta `$ = 6.0 are plotted as a function of $`am_q\frac{1}{2}\left(\frac{1}{\kappa }\frac{1}{\kappa _c}\right)`$. The $`B`$ parameters for all relevant operators are plotted: $`X`$ = $`L`$, $`R`$, $`S`$, $`N`$ and $`M`$ for $`B_{B,X}^{lat}`$, and $`X`$ = $`S`$, $`L`$, $`R`$, $`P`$, $`T`$ for $`B_{S,X}^{lat}`$. Each $`B`$ parameter is normalized by its vacuum saturation approximation. We averaged the matrix elements with $`X`$ = $`L`$ and $`R`$, as they should be equal in the infinite number of statistics by parity symmetry. In later sections, the averaged matrix elements are denoted by $`X`$ = $`LR`$.

We extrapolate these $`B`$ parameters to the chiral limit of light quark assuming a linear function in $`am_q`$. In most cases the chiral extrapolation changes the value of $`B`$ parameters from the lightest measured data by about 1% or less. Therefore, the chiral extrapolation is extremely stable and the associated systematic error is negligible. To confirm this observation we also tried a quadratic extrapolation for some parameter sets, for which we find that the results are consistent with the linear extrapolation within the statistical error.

In chiral perturbation theory for heavy-light mesons, the logarithmic dependence such as $`m_q\mathrm{ln}m_q`$ is predicted for $`B_B`$ Booth:1995hx ; Sharpe:1996qp . In the quenched approximation the chiral limit is even divergent as $`\mathrm{ln}m_q`$. The more divergent term $`\mathrm{ln}m_q`$ has a factor $`13g^2`$ as its coefficient, and the $`B^{}B\pi `$ coupling $`g`$ is evaluated in the range 0.2–0.7 using $`D^{}D\pi `$ decay Stewart:1998ke , $`D^{}`$ decay width Anastassov:2002cw , quark models Casalbuoni:1997pg , and quenched lattice calculations deDivitiis:1998kj ; Aoki:2001rd ; Abada:2002xe . It means that this divergent logarithm is relatively unimportant because of its small coefficient $`13g^2`$ = 0.2(7). It is, however, difficult to resolve such logarithmic dependences from the data taken in the range of our light quark masses. In this work, therefore, we do not further consider them, leaving the study of the chiral logarithm including the effect of unquenching for future publications.

Results of the linear extrapolation are summarized in Table 29, where we list the values of $`B_{B,X}^{lat}`$ and $`B_{S,X}^{lat}`$ at each $`\beta `$ and $`aM_0`$. The value of $`\kappa `$ corresponding to the physical $`u`$ or $`d`$ quark mass, which we denote $`\kappa _{ud}`$, is very close to the critical value $`\kappa _c`$. The value of $`\kappa _s`$ corresponding to the strange quark mass depends on the input quantity. We list the results at $`\kappa _{s1}`$, for which the $`K`$ meson mass is used as input, and at $`\kappa _{s2}`$, for which $`\varphi `$ meson mass is used.

### V.3 $`1/M_P`$ dependence

The $`1/M_P`$ dependence of the lattice $`B`$ parameters $`B_{B,X}^{lat}`$ ($`B_{S,X}^{lat}`$) is plotted in Figure 6 (Figure 7). The light quark is extrapolated to the chiral limit. Although the data at different $`\beta `$ values are overlaid, they do not necessarily agree because the operators are not renormalized. Comparison with an adequate definition (the $`B`$ parameter in the continuum renormalization scheme) will be discussed in detail in Section V.6.

We find that the mass dependence is small for $`B_{B,LR}^{lat}`$, while it is significant for others. This behavior can be mostly understood using the vacuum saturation approximation (VSA) Hashimoto:1999ck . In VSA the matrix element in the numerator of the $`B`$ parameter is generally decomposed into $`|0|A_\mu |P|^2`$ and $`|0|P|P|^2`$. For $`B_{B,LR}^{lat}`$, however, it is written by $`|0|A_\mu |P|^2`$ only and no $`|0|P|P|^2`$ term appears by definition, so that $`B_{B,LR}^{(\mathrm{VSA})}`$ = 1 is independent of $`1/M_P`$. For others, the term $`|0|P|P|^2`$ gives a strong mass dependence proportional to $`(M_P/M)^2=(1+\overline{\mathrm{\Lambda }}/M)^2`$, where $`\overline{\mathrm{\Lambda }}`$ represents the binding energy produced by the light degrees of freedom. Comparison of the lattice data with VSA is made in Ref. Hashimoto:1999ck .

### V.4 Renormalized $`B`$ parameters

The $`B`$ parameters for the continuum operators are obtained from the lattice $`B`$ parameters using (49) for $`B_B`$ and (50) for $`B_S`$. We consider the truncation method I in this subsection. The results of other truncations are discussed in the next subsection.

In order to see the effect of $`1/M`$ corrections we consider the quantity

$$\mathrm{\Phi }_{B_B}(\mu _b)\left(\frac{\alpha _s(M_P)}{\alpha _s(M_B)}\right)^{2/\beta _0}B_B(\mu _b)$$
(66)
as a function of $`1/M_P`$. The factor $`(\alpha _s(M_P)/\alpha _s(M_B))^{2/\beta _0}`$ is introduced to cancel the logarithmic dependence on $`M`$ coming from the continuum one-loop integral, so that the heavy quark expansion in $`1/M_P`$ is explicit. Up to two-loop corrections $`\mathrm{\Phi }_{B_B}(\mu _b)`$ is equivalent to $`B_B(\mu _b)`$ obtained with a replacement of $`M`$ in $`\rho _{L,L}`$ (43) and in $`\rho _A^{(0)}`$ (28) by the physical $`b`$ quark mass $`m_b`$, which can be confirmed by expanding the factor $`(\alpha _s(M_P)/\alpha _s(M_B))^{2/\beta _0}`$ in $`\alpha _s(M_P)`$ explicitly. With the replacement the static limit ($`1/M_P0`$) of NRQCD simply becomes the conventional static approximation. Therefore, in the calculation of $`\mathrm{\Phi }_{B_B}(\mu _b)`$ we explicitly set the physical $`b`$ quark mass $`m_b`$ in the matching coefficients (43) and (28). It should be also noted that at the physical $`B`$ meson mass, namely $`M_P`$=$`M_B`$, our definition of $`\mathrm{\Phi }_{B_B}(\mu _b)`$ exactly agrees with the definition (66).

Figure 8 (top panel) shows $`\mathrm{\Phi }_{B_B}(\mu _b)`$ at $`\beta =6.0`$. The light quark mass is extrapolated to the chiral limit, and the renormalization scale $`\mu _b`$ is set to $`m_b`$. In the one-loop matching (49) we use the renormalized coupling $`\alpha _V(q^{})`$ defined through the heavy quark potential Lepage:1993xa . The scale $`q^{}`$ represents the momentum region where the relevant one-loop integral dominates. Since it is not known, we use three typical values $`\pi /a`$, $`2/a`$ and $`1/a`$ and consider their variation as an indication of systematic uncertainty from higher order perturbative corrections. We find that the variation among different coupling constants becomes substantial as one goes to the static limit, while it is relatively unimportant in the physical mass region $`1/M_P`$ 0.2 GeV-1. This is because the one-loop coefficients in the matching (49) grows toward the static limit. For the $`1/M_P`$ dependence of $`\mathrm{\Phi }_{B_B}(\mu _b)`$ we observe a slight positive slope and curvature, but the large systematic uncertainty implies that the mass dependence is insignificant.

We obtain a similar plot for $`\mathrm{\Phi }_{B_S}(\mu _b)`$ in Figure 8 (bottom panel), which is an analog of $`\mathrm{\Phi }_{B_B}(\mu _b)`$ but for $`B_S(\mu _b)`$. The definition of $`\mathrm{\Phi }_{B_S}(\mu _b)`$ with the renormalization group improvement as in (66) is more complicated, since the logarithmic dependence appears in more than one coefficients, so that we have to consider a mixing of operators. In this work, however, we avoid this problem by replacing the heavy quark mass in $`\rho _{S,S}`$ and $`\rho _{S,L}`$ by the physical value as we did for $`\mathrm{\Phi }_{B_B}(\mu _b)`$.

For $`B_S(\mu _b)`$ the one-loop coefficients are relatively small and their dependence on the heavy quark mass is mild. Hence, we obtain a smaller variation due to different scale settings in the coupling constant.

### V.5 Effect of truncation of expansions

As we discussed in Section III.4 there are several method to truncate the perturbative and $`1/M`$ expansions. We consider the following four methods. In the methods I and I’ the perturbative matching is truncated in the numerator and denominator separately as in (49) and (50), while in the methods II and II’ the denominator is linearized as (53) and (54). In the primed methods the heavy quark field without the FWT rotation (22) is used for one-loop correction terms.

In Figure 9 we plot $`\mathrm{\Phi }_{B_B}(\mu _b)`$ (top panel) and $`\mathrm{\Phi }_{B_S}(\mu _b)`$ (bottom panel) for four different truncation methods. As can be seen from the figure, the methods I and I’ (or II and II’) agree in the static limit, since their difference is only in the FWT rotation. On the other hand, the difference between the methods I and II (or I’ and II’) is smaller for lighter heavy quarks because the one-loop correction in the denominator $`\rho _A^{(0)}`$ becomes small.

We consider the variation among different truncations as an indication of systematic uncertainties from higher orders of perturbative and $`1/M`$ expansions. The error estimation is given in the next section.

### V.6 Results at different lattice spacings

Since NRQCD is formulated by an expansion in $`1/M`$ and not a renormalizable field theory, it does not allow a continuum limit; instead it has to be considered as an effective theory valid in the region where $`1/(aM_0)`$ is small enough. The dependence of systematic errors on the lattice spacing $`a`$ is not just a simple power series in $`a`$, but contains its inverse powers. Therefore, the question is how one can find a region of $`a`$ where the discretization error is small while the errors scaling as $`1/a^n`$ ($`n`$ is a positive integer) is under control. Although the order counting argument as discussed in Hashimoto:1999ck ; Hashimoto:2000eh  provides a rough estimate of errors, it is essential to confirm it using actual simulation data.

In Figure 10 we plot $`\mathrm{\Phi }_{B_B}(\mu _b)`$ (top panel) and $`\mathrm{\Phi }_{B_S}(\mu _b)`$ (bottom panel) obtained with the truncation method I for four different lattice spacings. The largest (smallest) inverse lattice spacing is 2.3 GeV at $`\beta =6.1`$ (1.1 GeV at $`\beta =5.7`$). We find that around the physical $`B`$ meson mass ($`1/M_B`$ 0.2 GeV-1) the results agree within order 10% for $`\mathrm{\Phi }_{B_B}(\mu _b)`$ or even better for $`\mathrm{\Phi }_{B_S}(\mu _b)`$. The agreement becomes marginal toward the static limit especially for the nonperturbatively improved lattice ($`\beta =6.0`$), but it is not statistically significant.

Results of physical $`B_B(m_b)`$ (top panel) and $`B_S(m_b)`$ (bottom panel) are plotted in Fig. 11 as a function of the lattice spacing, where the variation due to the different choice of fit range is added to the error bar at each $`\beta `$. The $`a`$ dependence is a mixture of the discretization error scaling as $`a^m`$ and the truncation error containing a form like $`1/a^n`$. In addition, the truncation of perturbative expansion gives a functional dependence like $`1/\mathrm{ln}a`$. It is, therefore, difficult to determine the shape of the $`a`$ dependence, but the data imply that none of these errors is diverging in the region we measured.

## VI Physics results

### VI.1 Analysis of systematic errors

As discussed in the previous section, we have performed the calculation of the $`B`$ parameters with four different truncations of $`1/M`$ and $`\alpha _s`$ expansions. Furthermore, the calculations are made at four different lattice spacings. All of these calculations have different amount of various systematic errors, and thus they allow us to estimate the uncertainty in our final results. In this subsection we first list possible sources of systematic errors and estimate their size using a naive order counting. Then, their results are compared with the actual lattice data.

One of the possible systematic errors arises from the discretization of derivatives, which scales as a power of the lattice spacing $`a`$. Because our actions and operators are $`O(a)`$-improved at tree level, the leading error is of order $`a^2`$ and of order $`\alpha _sa`$. Since we are using an effective theory for heavy quark, the truncation of the $`1/M`$ expansion leads to a systematic error. For our choice of actions and operators the leading contribution is of order $`1/M_b^2`$. Again, since the matching of the $`1/M`$ terms is done at the tree level only, we also expect an error of order $`\alpha _s/M_b`$. The perturbative matching of operators are truncated at the one-loop level, so that there is an uncertainty of order $`\alpha _s^2`$. In Table 10 we list their typical size at each $`\beta `$ value using a naive order counting. Where the scale is needed we assume the typical spatial momentum inside a heavy-light meson to be $`p\mathrm{\Lambda }_{\mathrm{QCD}}`$ 300 MeV. For the strong coupling constant we use a typical value $`\alpha _V(2/a)`$ as listed in Table 1.

The contribution of order $`p^2/M_b^2`$ is not investigated in this paper, as we just neglect the $`1/M^2`$ terms in the action and operators. In Table 10 we assign 2% as the corresponding uncertainty rather than a naive order counting 0.4%, taking the estimate from explicit lattice study in Hashimoto:1999ck .

Since we have removed errors of $`\alpha _s/(aM)^m`$ ($`m`$0) by perturbative matching, the leading contribution which prevents us from the continuum extrapolation with the NRQCD action has the form $`\alpha _s^2/(aM)`$. Although its size in the naive order counting is smaller than the pure two-loop correction $`\alpha _s^2`$, we include it in our error analysis (thus in Table 10), as it gives the leading contribution growing toward the continuum limit.

As mentioned in Section III.4, the results from the different truncation methods (I, II, I’ and II’) are expected to differ from each other by $`O(\alpha _s^2)`$ or $`O(\alpha _sp/M_b)`$. We compare their results in Figure 12 for $`B_B`$ and in Figure 13 for $`B_S`$. The results of the four truncation methods when $`\alpha _s(2/a)`$ is used in the one-loop matching are plotted. In these Figures we also show the size of the systematic errors of order $`O(\alpha _s^2)`$ and $`O(\alpha _sp/M_b)`$ estimated with the naive order counting (first two lines of Table 10 added in quadrature). Although the statistical error of the data points makes the comparison somewhat ambiguous, we conclude that the naive order counting reasonably explains the scatter of the lattice data. The same conclusion is reached when we use different values of $`q^{}`$ in the perturbative matching.

An alternative way to estimate the $`O(\alpha _s^2)`$ error is to see the variation of results with $`q^{}`$=$`\pi /a`$, $`2/a`$ and $`1/a`$. From Figure 8, where the data at $`\beta `$ = 6.0 are plotted, we find that the variation with the choice of scale in the coupling constant is consistent with our order counting ($``$ 3%) for the data points around the $`B`$ meson mass. In the static limit, on the other hand, the one-loop coefficient is uncomfortably large, and the variation among the results with different $`q^{}`$ is much larger than our naive estimate.

Other sources of systematic errors to be tested are $`O(a^2p^2)`$, $`O(ap\alpha _s)`$ and $`O(\alpha _s^2/(aM_b))`$. Although it is difficult to disentangle various $`a`$ dependent systematic errors solely from the data, our results are stable against the change of lattice spacing suggesting that the systematic error is well estimated by the naive order counting.

Finally we also investigated the systematic error from the contamination of excited states by changing the fit range for the ratio $`R_{B,X}(t_1,t_2)`$ and $`R_{S,X}(t_1,t_2)`$. We find that the effect is at most 3%. In particular, at $`\beta `$=6.0 it is found to be 1% for $`B_B`$ and 3% for $`B_S`$. These variations are taken into account in the final results.

### VI.2 $`B`$ parameter results

As we discussed above, the systematic errors estimated with the naive order counting actually describe the differences among different calculations. We, therefore, use the order counting argument to quote our estimate of systematic uncertainties.

We take the central values from the data at $`\beta `$ = 6.0 (set $`B`$) with $`q^{}`$=$`2/a`$, and obtain the following results in the quenched approximation.

$`B_B(m_b)`$ $`=`$ $`0.84(3)(5),`$ (67)
$`B_{B_s}(m_b)`$ $`=`$ $`0.86(2)(5)(0),`$ (68)
$`B_S(m_b)`$ $`=`$ $`0.82(2)(5),`$ (69)
$`B_{S_s}(m_b)`$ $`=`$ $`0.85(1)(5)(_0^{+1}),`$ (70)
where the first and second errors represent statistical and systematic ones respectively. In the systematic error the contamination from excited states is added in quadrature with other sources estimated by the order counting. The third error is from the uncertainty of $`m_s`$ arising from the different input physical quantities, i.e. $`m_K`$ or $`m_\varphi `$.

The corresponding renormalization scale independent $`B`$ parameter is obtained from (7)

$`\widehat{B}_B=1.29(5)(8),`$ (71)
which may be compared with the previous calculations using the relativistic actions for heavy quark, 1.38(11)$`(_9^{+0})`$ Becirevic:2001nv  and 1.40(5)$`(_1^{+6})`$ Lellouch:2001tw . These two results are slightly higher than our result, although Becirevic:2001nv  is still consistent within the large error. A possible reason for the high values in the relativistic approach is in the extrapolation in the heavy quark mass from the charm quark mass region to the bottom. In fact, the combined analysis of the HQET and relativistic heavy quark, in which the interpolation in $`1/M`$ can be made, gives 1.34(6)$`(_6^{+8})`$ Becirevic:2001xt , i.e. closer to our result.

### VI.3 Applications

In this subsection we present a few examples of physics applications of our results. It should be noted, however, that our calculation is still in the quenched approximation and there is no rigorous estimate for the associated uncertainty. For the following applications we assume that the quenching effect is negligible for the $`B`$ parameters as suggested by our preliminary calculations Hashimoto:2001zq ; Yamada:2001ik ; Yamada:2001xp .

For the $`B`$ meson decay constant, on the other hand, the large effect of quenching has been found Ryan:2001ej ; Kronfeld:2001ss ; Yamada:2002\_rev . Furthermore, a large uncertainty due to the presence of chiral logarithm is suggested for $`f_B`$ Yamada:2002\_rev ; Kronfeld:2002ab , while the effect is not too large for $`f_{B_s}`$. In the following analysis we therefore consider the quantities for which only $`f_{B_s}`$ is needed, and use the recent world average of unquenched lattice calculations of $`f_{B_s}`$ = 230(30) MeV Ryan:2001ej ; Yamada:2002\_rev  when needed.

Assuming the three generation unitarity $`|V_{ts}||V_{cb}|`$, we obtain the mass difference in $`B_s^0\overline{B}_s^0`$ mixing using (2) as

$$\mathrm{\Delta }M_s=19.4(5.5)\text{ps}^1,$$
(72)
where the statistical and systematic errors in theoretical and experimental quantities are added in quadrature, but the final error is dominated by the uncertainty of $`f_{B_s}`$. The value is consistent with the current lower bound $`\mathrm{\Delta }M_s>13.1`$ ps-1 at a 95% CL PDG2002 . Tevatron Run II is expected to measure the mass difference very precisely in a few years.

The width difference in the $`B_s`$ meson system could also be measured at Tevatron Run II, if it is large enough. Using (15) and (70), we obtain

$$\left(\frac{\mathrm{\Delta }\mathrm{\Gamma }}{\mathrm{\Gamma }}\right)_{B_s}=0.106\pm 0.020\pm 0.028\pm 0.037\pm 0.024,$$
(73)
where the first through third errors are from $`m_b`$ = 4.8(3) GeV, $`f_{B_s}`$ = 230(30) MeV, and $``$ (or $`\overline{m}_b(m_b)`$ = 4.25(25) GeV and $`\overline{m}_s(m_b)`$ = 0.10(3) GeV), respectively. The last error comes from the uncertainty in the estimation of the $`1/m`$ correction, for which we assign 30%. The error from the $`B`$ parameters is much smaller than the others listed above.

The uncertainty in the calculation of $`(\mathrm{\Delta }\mathrm{\Gamma }/\mathrm{\Gamma })_{B_s}`$ is still very large ($``$ 50% if added in quadrature). In order to improve it one has to calculate the $`1/m`$ corrections reliably, as it largely cancels the leading contribution from $`B_S`$ as seen in (15). Currently, only an upper bound is obtained for this quantity experimentally. Our prediction is consistent with the bound $`(\frac{\mathrm{\Delta }\mathrm{\Gamma }}{\mathrm{\Gamma }})_{B_s}<0.31`$ at a 95 % CL PDG2002 .

## VII SU(3) breaking ratio $`\xi `$

Since we expect that the bulk of systematic uncertainty in the calculation of $`f_B`$ and $`B_B`$ cancels in their SU(3) breaking ratio $`f_{B_s}/f_B`$ and $`B_{B_s}/B_B`$, they could be useful to reduce the errors in the determination of $`|V_{td}|`$ through the relation (8).

In the lattice calculation, the deviation of their ratio from unity is the quantity to be calculated and the errors scale as $`B_{B_s}/B_B1`$ rather than $`B_{B_s}/B_B`$ itself. In the present case, the naive estimate of SU(3) breaking is $`O(m_K^2/\mathrm{\Lambda }_\chi ^2)`$ 25%, where $`\mathrm{\Lambda }_\chi `$ is a scale of the chiral symmetry breaking $``$ 1 GeV, and the order counting of uncertainties for the ratio are starting from this order.

As done in Section VI.1, we compare our order counting with data in Figure 14. As we expected, the variation with different truncations is much smaller than the expected systematic errors (dotted line). The dependence on $`\beta `$ is sizable, but not significant compared to the relatively large statistical error.

Our result is

$$\frac{B_{B_s}}{B_B}=1.020(21)(_{16}^{+15})(_0^{+5}),$$
(74)
where the central value is taken from $`\beta `$ = 6.0 (set B), and the errors are statistical, systematic and the uncertainty of $`m_s`$ in the order given.

For the calculation of $`\xi `$ (1) we have to combine $`B_{B_s}/B_B`$ with $`f_{B_s}/f_B`$, which is called the indirect method. On the other hand, it is also possible to directly obtain $`\xi `$ from a ratio of $`O_{L_s}`$ and $`O_{L_d}`$ (direct method). It has been discussed that one may obtain rather large value of $`\xi `$ if one uses the direct method Bernard:1998dg ; Lellouch:2001tw . Therefore, in the following we check if we could obtain consistent results from both methods using our data.

Figure 15 shows the chiral extrapolation of $`O_L`$, as required in the direct method. The data are obtained at $`\beta `$ = 6.0 (set B) for a heavy quark mass closest to the $`b`$ quark mass. The dashed line is obtained by a linear fit to the data (open circles), while the solid curve represents a fit with a linear plus quadratic term in $`am_q`$. Although the data look consistent with the linear fit, the chiral limit with the quadratic fit is higher by about one standard deviation.

An open diamond at the chiral limit, on the other hand, is obtained through the indirect method, i.e. the decay constant and $`B`$ parameter are separately extrapolated to the chiral limit with a linear fit. Although we have not presented a calculation of the decay constant in this paper, they are done on the same set of gauge configurations at $`\beta `$ = 6.0 and the lattice axial current is renormalized as described in Section III.2. The result is completely consistent with the quadratic fit in the direct method. It implies that in the direct method the chiral extrapolation is more difficult and needs enough statistics to control, since $`O_L`$ is effectively the decay constant squared so that the finite $`am_q`$ correction is amplified888 The similar discussion may be found in Lellouch:2001tw . .

## VIII Conclusions

In this work, we calculate the $`B`$ meson $`B`$ parameters on the lattice in the quenched approximation. The calculation is an extension of our previous works Hashimoto:1999ck ; Hashimoto:2000yh , in which $`B_B`$ and $`B_S`$ were calculated for the first time with the lattice NRQCD action.

In the present work we include a detailed study of systematic uncertainties. Using the lattice simulations at four different $`\beta `$ values with the $`O(a)`$-improved actions, we find that the $`B`$ parameters are essentially insensitive to the discretization error. We also investigate the systematic errors associated with the truncation of heavy quark and perturbative expansions, which are necessary in the effective theory approach such as NRQCD. By comparing four different truncations of these expansions, we are able to confirm that the naive order counting argument of the systematic errors could actually give a reasonable estimate.

In our final results for the $`B`$ parameters the systematic error is $``$ 6%, which is already smaller than that in the equivalent calculations of $`f_B`$ (10–20%), owing to the fact that it is defined as a ratio to the vacuum saturation approximation. Further reduction of systematic errors, if it is necessary, requires higher order calculation of perturbation theory and the $`O(a^2)`$ improvement. Approaching to the continuum limit will not help to reduce the errors in the NRQCD approach.

For a precise extraction of the important CKM element $`|V_{td}|`$ through the SU(3) breaking ratio of the $`B\overline{B}`$ mixing one needs $`\xi ^2`$. It is preferable to take the chiral limit separately for $`f_B`$ and $`B_B`$, as they have milder light quark mass dependence. The SU(3) breaking ratio of $`B_B`$ is obtained with accuracy of order a few percent, since the $`B`$ parameter is extremely insensitive to the light quark mass and the large cancellation of systematic errors is expected.

The largest remaining uncertainty in our calculation is in the quenching approximation, though it is not explicitly discussed in the paper. We are currently performing an unquenched simulation with the same lattice action at similar lattice spacing, which will allow us to directly study the quenching effect.

###### Acknowledgements.

This work is supported by the Supercomputer Project No. 79 (FY2002) of High Energy Accelerator Research Organization (KEK), and also in part by the Grant-in-Aid of the Ministry of Education (Nos. 11640294, 12640253, 12740133, 13135204, 13640259, 13640260, 14046202, 14740173). N.Y. is supported by the JSPS Research Fellowship.