# A lattice NRQCD calculation of the 𝐵⁰-𝐵̄⁰ mixing parameter 𝐵_𝐵

## I Introduction

The constraints on the unitarity triangle of the Cabibbo-Kobayashi-Maskawa (CKM) matrix can provide us with one of the most crucial information on the physics beyond the Standard Model . However, due to large theoretical or experimental uncertainties, the current bound is too loose to test the Standard Model or new physics. The $`B^0`$-$`\overline{B^0}`$ mixing sets a constraint on $`|V_{td}V_{tb}^{}|`$ through the currently available experimental data on the mass difference between two neutral $`B`$ mesons $`\mathrm{\Delta }M_B`$= 0.477$`\pm `$0.017 ps<sup>-1</sup> . The experimental achievement is impressive as the error is already quite small $``$ 4%. Theoretical calculation to relate $`\mathrm{\Delta }M_B`$ to $`|V_{td}V_{tb}^{}|`$, on the other hand, involves a large uncertainty in the $`B`$ meson matrix element $`\overline{B^0}|𝒪_L|B^0`$, which requires a method to calculate the non-perturbative QCD effects.

Lattice QCD is an ideal tool to compute such non-perturbative quantities from first principles. There has been a number of calculations of the $`B`$ meson decay constant $`f_B`$ and the $`B`$-parameter $`B_B`$ ($`B_B`$ describes the matrix element through $`\overline{B^0}|𝒪_L|B^0`$ = $`(8/3)B_Bf_B^2M_B^2`$). The calculation of $`f_B`$ is already matured at least in the quenched approximation . Major systematic errors are removed by introducing the non-relativistic effective actions and by improving the action and currents. Remaining $`a`$ (lattice spacing) dependent systematic error is confirmed to be small, and in some papers the continuum extrapolation is made. A recent review  summarized the value of $`f_B`$ as $`f_B`$= 165$`\pm `$20 MeV within the quenched approximation.

An essential ingredient of these calculations is the use of the non-relativistic effective actions. Since the $`b`$-quark mass in lattice unit $`am_b`$ is large for the typical lattices used for simulations, the relativistic (Wilson-type) actions could suffer from a large discretization error of order $`am_b`$ or $`(am_b)^2`$. The non-relativistic QCD (NRQCD) , on the other hand, is formulated as an expansion in $`𝒑/m_Q`$. In the heavy-light meson system, where the typical spatial momentum scale is $`\mathrm{\Lambda }_{QCD}`$, the error from the truncation of higher order terms is controllable. The calculation of $`f_B`$ is now available to order $`1/m_b^2`$, and it is known that the contribution of $`O(1/m_b)`$ is significant ($`20\%`$) while that of $`O(1/m_b^2)`$ is small ($`3\%`$. In addition, the calculations based on the Fermilab approach for heavy quark  agree with the NRQCD results  including their $`1/m_Q`$ dependence. These results make us confident about the non-relativistic effective action approaches in Lattice QCD. Now that the computation of $`f_B`$ is established, the next goal is to apply the similar technique to the computation of $`B_B`$.

The lattice calculations of the $`B`$-parameter have been done in the infinitely heavy quark mass (static) limit , and the results are in reasonable agreement with each other. There is, however, some indication that the $`1/m_b`$ correction would be non-negligible from the study with relativistic quark actions . Their results show that there is small but non-zero negative slope in the $`1/m_Q`$ dependence of $`B`$-parameter, but it is not conclusive because of the possible systematic uncertainty associated with the relativistic quark action for heavy quark. The purpose of this work is to study the $`1/m_Q`$ dependence of $`B_B`$ by explicitly calculating it with the NRQCD action at several values of $`1/m_Q`$. Our result confirms the previous works : there is a small negative slope in $`B_B`$. In addition, we find that the slope comes from the large $`1/m_Q`$ dependence of $`B_N^{lat}`$ and $`B_S^{lat}`$, which are matrix elements of non-leading operators. For the observed $`1/m_Q`$ dependence of the lattice matrix elements, qualitative explanations are given in Discussion section using the vacuum saturation approximation.

The perturbative matching of the continuum and lattice operators introduces a complication to the analysis. Since the one-loop coefficients for four-quark operators are not known yet for the NRQCD action, we use the coefficients in the static limit in Refs. . The systematic error associated with this approximation is of $`\alpha _s/(am_Q)`$ and expected to be small compared to the $`1/m_Q`$ correction itself. This and other systematic errors are discussed in detail in the Discussion section.

This paper is organized as follows. In the next section, we summarize our matching procedure. The simulation method is described in section III, and our analysis and results are presented in section IV. The results are compared with the vacuum saturation approximation on the lattice in section V, and we estimate the remaining uncertainties in section VI. Finally our conclusion is given in section VII. An early result of this work was presented in Ref. .

## II Perturbative matching

In this section, we give our notations and describe the perturbative matching procedure. The mass difference between two neutral $`B_q`$ mesons is given by

$$\mathrm{\Delta }M_{B_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 )\right]^{6/23}\left[1+\frac{\alpha _s(\mu )}{4\pi }J_5\right]\overline{B_q^0}|𝒪_L(\mu )|B_q^0,$$
(1)   
where $`q=d`$ or $`s`$ and $`S_0(x_t)`$ ($`x_t=m_t^2/m_W^2`$) and $`\eta _{2B}`$ are so-called Inami-Lim function and the short distance QCD correction, respectively. Their explicit forms can be found in Ref. . $`𝒪_L(\mu )`$ is a $`\mathrm{\Delta }B`$=2 operator

$$𝒪_L(\mu )=\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1\gamma _5)q,$$
(2)   
renormalized in the $`\overline{\mathrm{MS}}`$ scheme with the naive dimensional regularization (NDR). $`J_{n_f}`$ is related to the anomalous dimension at the next-to-leading order with $`n_f`$ active flavors as

$$J_{n_f}=\frac{\gamma ^{(0)}\beta _1}{2\beta _0^2}\frac{\gamma ^{(1)}}{2\beta _0},$$
(3)   
where

$$\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}$$
(4)   
$`n_f`$=5 when $`\mu `$ is greater than or equal to the $`b`$ quark mass.

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

$$B_{B_q}(\mu )=\frac{\overline{B_q^0}|𝒪_L(\mu )|B_q^0}{\frac{8}{3}\overline{B_q^0}|A_\mu |00|A_\mu |B_q^0},$$
(5)   
where $`A_\mu `$ denotes the axial-vector current $`\overline{b}\gamma _\mu \gamma _5q`$. The renormalization invariant $`B`$-parameter is defined by

$$\widehat{B}_{B_q}=\left[\alpha _s(\mu )\right]^{\frac{6}{23}}\left[1+\frac{\alpha _s(\mu )}{4\pi }J_5\right]B_{B_q}(\mu ),$$
(6)   
which does not depend on the arbitrary scale $`\mu `$ up to the next-to-leading order. The scale $`\mu `$ is conventionally set at the scale of $`b`$-quark mass $`\mu =m_b`$.

In order to calculate the matrix element $`\overline{B_q^0}|𝒪_L(m_b)|B_q^0`$ on the lattice, we have to connect the operator $`𝒪_L(m_b)`$ defined in the continuum renormalization scheme with its lattice counterpart. The matching coefficients can be obtained by requiring the perturbative quark scattering amplitudes at certain momentum with continuum $`O_L`$ operator and with lattice four fermi operators should give identical results.

At the one-loop level the matching gives the following relation

$`𝒪_L(m_b)`$ $`=`$ $`\left(1+{\displaystyle \frac{\alpha _s}{4\pi }}[4\mathrm{ln}(a^2m_b^2)+D_L14]\right)𝒪_L^{lat}(1/a)`$ (8)   
$`+{\displaystyle \frac{\alpha _s}{4\pi }}D_R𝒪_R^{lat}(1/a)+{\displaystyle \frac{\alpha _s}{4\pi }}D_N𝒪_N^{lat}(1/a)+{\displaystyle \frac{\alpha _s}{4\pi }}D_S𝒪_S^{lat}(1/a),`$   
where $`𝒪_{\{L,R,N,S\}}^{lat}`$ denotes the naive local operators defined on the lattice in which the light quarks are not rotated. Their explicit forms are the following,

$`𝒪_R`$ $`=`$ $`\overline{b}\gamma _\mu (1+\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)q,`$ (9)   
$`𝒪_N`$ $`=`$ $`\overline{b}\gamma _\mu (1\gamma _5)q\overline{b}\gamma _\mu (1+\gamma _5)q+\overline{b}\gamma _\mu (1+\gamma _5)q\overline{b}\gamma _\mu (1\gamma _5)q`$ (11)   
$`+2\overline{b}(1\gamma _5)q\overline{b}(1+\gamma _5)q+2\overline{b}(1+\gamma _5)q\overline{b}(1\gamma _5)q,`$   
$`𝒪_S`$ $`=`$ $`\overline{b}(1\gamma _5)q\overline{b}(1\gamma _5)q.`$ (12)   
Unfortunately, the one-loop coefficients $`D_{\{L,R,N\}}`$ for the NRQCD heavy and $`O(a)`$-improved light quark action  are not known yet. In this work, we use the one-loop coefficients in the static limit ,

$$D_L=21.16,D_R=0.52,D_N=6.16,D_S=8.$$
(13)   
The systematic error associated with this approximation is at most $`\alpha _s/(am_Q)`$, since the NRQCD action’s $`m\mathrm{}`$ limit agrees with the static action. The numerical size of the error is discussed later.

The matching of the axial-vector current appearing in the denominator of Eq. (5) can be done in a similar manner 

$`A_0`$ $`=`$ $`\left(1+{\displaystyle \frac{\alpha _s}{4\pi }}[2\mathrm{ln}(a^2m_b^2)+D_A{\displaystyle \frac{8}{3}}]\right)A_0^{lat}(1/a),`$ (14)   
where $`A_0^{lat}`$ is defined on the lattice, and the matching coefficient $`D_A`$ in the static limit $`D_A=13.89`$.

In calculating the ratio of Eq. (5) a large cancellation of perturbative matching corrections takes place between the numerator and denominator, since the large wave function renormalization coming from the tadpole contribution in the lattice theory is the same. To make this cancellation explicit we consider the matching of a ratio $`B_B(m_b)=𝒪_L(m_b)/(8/3)A_0^2`$ itself,

$`B_B(m_b)`$ $`=`$ $`Z_{L/A^2}(m_b;1/a)B_L^{lat}(1/a)+Z_{R/A^2}(m_b;1/a)B_R^{lat}(1/a)`$ (16)   
$`+Z_{N/A^2}(m_b;1/a)B_N^{lat}(1/a)+Z_{S/A^2}(m_b;1/a)B_S^{lat}(1/a),`$   
where $`B_{\{L,R,N,S\}}^{lat}(1/a)=𝒪_{\{L,R,N,S\}}^{lat}(1/a)/(8/3)A_0^{lat}(1/a)^2`$ and the $`B`$ and $`\overline{B}`$ states are understood for the expectation values $`\mathrm{}`$ as in Eq. (5). Then the coefficients become

$`Z_{L/A^2}(m_b;1/a)`$ $`=`$ $`\left(1+{\displaystyle \frac{\alpha _s}{4\pi }}(D_L2D_A{\displaystyle \frac{26}{3}})\right),`$ (17)   
$`Z_{R/A^2}(m_b;1/a)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}D_R,`$ (18)   
$`Z_{N/A^2}(m_b;1/a)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}D_N,`$ (19)   
$`Z_{S/A^2}(m_b;1/a)`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{4\pi }}D_S.`$ (20)   
The Eqs. (17)-(20) are used in the following analysis to obtain $`B_B(m_b)`$. Numerical values of $`Z_{\{L,R,N,S\}/A^2}(m_b;1/a)`$ are given in Table I for the lattice parameters in our simulation. For the coupling constant $`\alpha _s`$ in Eqs. (17)-(20) we use the V-scheme coupling  with $`q^{}=\pi /a`$ or $`q^{}=1/a`$. At $`\beta `$=5.9 those are $`\alpha _V(\pi /a)`$=0.164 and $`\alpha _V(1/a)`$=0.270. The tadpole improvement  does not make any effect on the ratio of Eq. (16), since the tadpole contribution cancels between the numerator and denominator. The $`b`$-quark mass scale $`m_b`$ is set to 5 GeV as usual.

In the previous works in the static approximation , the leading and the next-to-leading logarithmic corrections are resummed to achieve better control in the running from $`m_b`$ to $`1/a`$. In this paper we use the one-loop formula without the resummation for simplicity. This does not introduce significant error, since the mass scale difference between $`m_b`$ and $`1/a`$ is small and the effect of the resummation is not important. In Appendix A we compare the $`Z`$ factors with and without the resummation.

We determine the heavy-light pseudo-scalar meson mass $`M_P`$ from the binding energy $`E^{\mathrm{bin}}`$ measured in the simulation using a formula

$$aM_P=Z_mam_QE_0+aE^{\mathrm{bin}},$$
(21)   
where $`Z_m`$ and $`E_0`$ are the renormalization constant for the kinetic mass and the energy shift, respectively. Both have been calculated perturbatively by Davies and Thacker  and by Morningstar . Since the precise form of their NRQCD action is slightly different from ours, we performed the perturbative calculations for our action. Our results for the coefficient $`A`$ and $`B`$ in the perturbative expansion

$`E_0`$ $`=`$ $`\alpha _V(q^{})A,`$ (22)   
$`Z_m`$ $`=`$ $`1+\alpha _V(q^{})B,`$ (23)   
are summarized in Table II.

## III Simulation Details

Our task is to compute the ratios $`𝒪_{\{L,R,N,S\}}^{lat}(1/a)/(8/3)A_0^{lat}(1/a)^2`$ using lattice NRQCD with the lattice spacing $`a`$. In this section we describe our simulation method to obtain them.

We performed the numerical simulation on a $`16^3\times 48`$ lattice at $`\beta `$= 5.9, for which the inverse lattice spacing fixed with the string tension is $`1/a`$= 1.64 GeV. In the quenched approximation we use 250 gauge configurations, each separated by 2,000 pseudo-heat bath sweeps. For the light quark we use the $`O(a)`$-improved action  at $`\kappa `$=0.1350, 0.1365. The clover coefficient is set to be $`c_{\mathrm{sw}}=1/u_0^3`$, where $`u_0P_{plaq}^{1/4}`$=0.8734 at $`\beta =5.9`$. The critical $`\kappa `$ value is $`\kappa _c=`$0.1401, and $`\kappa _s`$ corresponding to the strange quark mass determined from the $`K`$ meson mass is $`\kappa _s`$=0.1385.

For the heavy quark and anti-quark we use the lattice NRQCD action  with the tadpole improvement $`U_\mu U_\mu /u_0`$. The precise form of the action is the same as the one we used in the previous work . We use both of the $`O(1/m_Q)`$ and $`O(1/m_Q^2)`$ actions in parallel in order to see the effect of the higher order contributions. The heavy (anti-)quark field in the relativistic four-component spinor form is constructed with the inverse Foldy-Wouthuysen-Tani (FWT) transformation defined at the tadpole improved tree level as in Ref. . The heavy quark masses and the stabilization parameters are $`(am_Q,n)`$=(10.0,2), (5.0,2), (3.0,2), (2.6,2) and (2.1,3). These parameters approximately cover a mass scale between $`4m_b`$ and $`m_b`$.

We label the time axis of our lattice as $`t=[24,23]`$. The heavy quark and anti-quark propagators are created from a local source located at the origin ($`t`$=0 on our lattice) and evolve into opposite temporal directions. The light quark propagator is also solved with the same source location and with a Dirichlet boundary condition at $`t=24`$ and $`t=23`$. The $`B`$ and $`\overline{B}`$ mesons are constructed with local sink operators. Thus we have the four-quark operators at the origin and extract the matrix elements from the following three-point correlation function

$$C_X^{(3)}(t_1,t_2)=\underset{\stackrel{}{x}_1}{}\underset{\stackrel{}{x}_2}{}0|A_{0}^{lat}{}_{}{}^{}(t_1,\stackrel{}{x}_1)O_X^{lat}(0,\stackrel{}{0})A_{0}^{lat}{}_{}{}^{}(t_2,\stackrel{}{x}_2)|0,$$
(24)   
where $`X`$ denotes $`L`$, $`R`$, $`N`$ or $`S`$. Because of a symmetry under parity transformation, $`C_L^{(3)}(t_1,t_2)`$ and $`C_R^{(3)}(t_1,t_2)`$ should exactly coincide in infinitely large statistics. Therefore we explicitly average them before the fitting procedure we describe below.

To obtain the ratios $`B_X^{lat}(1/a)`$ we also define the following two-point functions

$`C^{(2)}(t_1)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}0|A_0^{lat}(t_1,\stackrel{}{x})A_{0}^{lat}{}_{}{}^{}(0,\stackrel{}{0})|0,`$ (25)   
$`C^{(2)}(t_2)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}0|A_0^{lat}(0,\stackrel{}{0})A_{0}^{lat}{}_{}{}^{}(t_2,\stackrel{}{x})|0,`$ (26)   
and consider a ratio

$$\frac{C_X^{(3)}(t_1,t_2)}{\frac{8}{3}C^{(2)}(t_1)C^{(2)}(t_2)}|t_i|1\frac{\overline{P^0}|𝒪_X^{lat}(1/a)|P^0}{\frac{8}{3}\overline{P^0}|A_0^{lat}(1/a)|00|A_0^{lat}(1/a)|P^0}=B_X^{lat}(1/a),$$
(27)   
where $`P^0`$ denotes a heavy-light pseudo-scalar meson. The ground state $`P^0`$ meson is achieved in the large $`|t_i|`$ regime. Although we use the local operator for the sinks at $`t_1`$ and $`t_2`$, the ground state extraction is rather easier for finite $`am_Q`$ than in the static approximation, since the statistical error is much smaller for NRQCD . This is another advantage of introducing the $`1/m_Q`$ correction.

The physical $`B_B(m_b)`$ is obtained by extrapolating and interpolating each $`B_X^{lat}(1/a)`$ to the physical $`B`$ meson with $`\kappa `$ and $`m_Q`$, respectively before combining them as Eq. (16). The final result for $`B_B(m_b)`$ may also be obtained by combining the ratio of correlation functions before a constant fit. Namely we use the relation

$$B_P(m_b;t_1,t_2)=\underset{X=L,R,N,S}{}Z_{X/A^2}(m_b,1/a)\frac{C_X^{(3)}(t_1,t_2)}{\frac{8}{3}C^{(2)}(t_1)C^{(2)}(t_2)}|t_i|1B_P(m_b).$$
(28)   
Since the statistical fluctuation in the individual $`B_X^{lat}(1/a)`$ is correlated, the error is expected to be smaller with this method (We use the jackknife method for error estimation). Following Ref.  we refer to this method as the “combine-then-fit” method, while the usual one as Eq. (27) is called the “fit-then-combine” method in the rest of the paper.

## IV Simulation Results

We describe the simulation results in this section. The results are from the $`O(1/m_Q)`$ action unless we specifically mention.

### A Heavy-light meson mass

The binding energy of the heavy-light meson is obtained from a simultaneous fit of two two-point correlation functions. The numerical results are listed in Table III for each $`am_Q`$ and $`\kappa `$. Extrapolation of the light quark mass to the strange quark mass or to the chiral limit is performed assuming a linear dependence in $`1/\kappa `$.

The meson mass is calculated using the perturbative expression Eq. (21). The results with $`\alpha _V(\pi /a)`$ and with $`\alpha _V(1/a)`$ are also given in Table III.

### B $`B_X^{lat}(1/a)`$ and $`B_P(m_b)`$

Figures 1 and 2 show the $`t_1`$ dependence of $`B_P(m_b;t_1,t_2)`$ in the “combine-then-fit” method. The perturbative matching of the continuum and lattice theory is done with the $`V`$-scheme coupling $`\alpha _V(q^{})`$ at $`q^{}`$=$`\pi /a`$ (left) and $`1/a`$ (right). Their difference represents the effect of $`O(\alpha _s^2)`$. The signal is rather noisier at $`am_Q`$=5.0 (Fig. 1) than at $`am_Q`$=2.6 (Fig. 2) from the same reason as in the static limit. But, still, a reasonably good signal is observed even for large $`am_Q`$. A plateau in the $`t_1`$ dependence is reached around $`t_1`$= 8 $``$ 11 for both $`t_2`$ = $``$10 and $``$15. To be conservative we take $`|t_1|`$ as well as $`|t_2|`$ greater than 10 for the fitting region. All data points $`(t_1,t_2)`$ in 10 $`|t_1|`$ 13 and in 10 $`|t_2|`$ 13 are fitted with constant to obtain the result for $`B_P(m_b)`$. We confirm that except for the heaviest quark the results are stable within about one standard deviation under a change of the fitting region by at most two $`t_i`$ steps in the forward and backward direction. The numerical results are listed in Table IV.

The light quark mass ($`1/\kappa `$) dependence of $`B_P(m_b)`$ is presented in Fig. 3. Since its dependence is quite modest, we assume a linear dependence on $`1/\kappa `$ and extrapolate the results to the strange quark mass and to the chiral limit as shown in the plot. Results of the extrapolation are also listed in Table IV.

### C $`1/M_P`$ dependence

In Fig. 4 we plot $`B_P(m_b)`$ in the chiral limit, namely $`B_{P_d}(m_b)`$, as a function of $`1/M_P`$ in the physical unit. We take $`q^{}`$=$`\pi /a`$ (circles) and $`1/a`$ (squares) for the scale of $`\alpha _V`$. Regardless of the choice of the coupling, we observe a small but non-zero negative slope in $`1/M_P`$, which supports the previous results by Bernard, Blum and Soni  and by Lellouch and Lin  using the relativistic fermions.

To investigate the origin of the observed $`1/M_P`$ dependence, we look into the contributions of the individual operators $`𝒪_{\{L,R,N,S\}}^{lat}`$ through the “fit-then-combine” method with the same fitting region as before. We list the results for each $`B_X^{lat}`$ in Table IV. Figure 5 shows the $`1/M_P`$ dependence of $`B_L^{lat}(1/a)`$(=$`B_R^{lat}(1/a)`$), $`B_N^{lat}(1/a)`$ and $`B_S^{lat}(1/a)`$. While no significant $`1/M_P`$ dependence is observed in $`B_L^{lat}(1/a)`$, $`B_N^{lat}(1/a)`$ and $`B_S^{lat}(1/a)`$ have strong slope. Since their sign is opposite and the sign of the matching factors $`Z_{N/A^2}(m_b;1/a)`$ and $`Z_{S/A^2}(m_b;1/a)`$ (see Table I) is the same, a partial cancellation takes place giving a small negative slope for $`B_{P_d}(m_b)`$.

We also make a comparison of the results of the $`1/m_Q`$ action (circles) with those of the $`1/m_Q^2`$ action (triangles) in Fig. 5. There is a small difference between the two results in $`B_N^{lat}(1/a)`$ and in $`B_S^{lat}(1/a)`$ toward large $`1/M_P`$ ($`1/M_P`$ 0.2 GeV<sup>-1</sup>), which is consistent with our expectation that the difference is an $`O(\mathrm{\Lambda }_{QCD}/m_Q)^2`$ effect.

Previous results in the static approximation by Ewing et al. (diamond) , Gimenéz and Martinelli (triangle)  and Christensen, Draper and McNeile (circle)  are plotted with filled symbols at $`1/M_P=0`$. Although the $`\beta `$ value and the light quark action employed (the $`O(a)`$-improved action is used in Refs.  and unimproved action in Ref. ) are different, all the results are in good agreement with each other. A quadratic extrapolation (dashed line) using our $`1/m_Q^2`$ NRQCD result also does agree with these previous static results.

### D Result for $`B_B(m_b)`$

Combining the data for $`B_X^{lat}(1/a)`$ discussed above, we obtain $`B_{P_d}(m_b)`$ with the “fit-then-combine” method. We confirm that the difference in numerical results from both methods are completely negligible. Figure 6 shows the results of “fit-then-combine” method using $`\alpha _V(1/a)`$ with both actions. The comparison with the static results is also made in this plot, where only the statistical error in each calculation is considered and the same matching procedure as ours are applied. We again observe a consistent result.

Interpolating the above NRQCD results to the physical $`B`$ meson mass, we obtain the physical $`B_{B_d}(m_b)`$

$$B_{B_d}(m_b)=\{\begin{array}{c}0.78(3)(q^{}=\pi /a)\\ 0.72(3)(q^{}=1/a)\end{array}$$
(29)   
for the $`O(1/m_Q)`$ action, and

$$B_{B_d}(m_b)=\{\begin{array}{c}0.78(2)(q^{}=\pi /a)\\ 0.71(3)(q^{}=1/a)\end{array}$$
(30)   
for the $`O(1/m_Q^2)`$ action. The quoted error is statistical only. For the ratio of $`B_{B_s}/B_{B_d}`$, we obtain $`B_{B_s}/B_{B_d}`$ = 1.01(1) for $`q^{}=\pi /a`$ and $`B_{B_s}/B_{B_d}`$ = 1.02(1) for $`q^{}=1/a`$ from both actions.

## V Discussion

The strong $`1/M_P`$ dependence in $`B_X^{lat}(1/a)`$ observed in Fig. 5 can be roughly understood using the vacuum saturation approximation (VSA) on the lattice as explained below. Here it should be noted that the terminology of VSA we use here does not immediately mean $`B_B(m_b)=1`$.

The VSA for $`B_{L,R}^{lat}`$ is unity by construction. This is true even for finite $`1/M_P`$, and its prediction is shown by a straight line in Fig. 5(a). The NRQCD data is located slightly below the line ($``$0.9), but the mass dependence is well reproduced by the VSA.

For $`B_N^{lat}`$ and $`B_S^{lat}`$, we require a little algebra to explain their mass dependence under the VSA. Using the Fierz transformation and inserting the vacuum, we obtain

$`\overline{P^0}|𝒪_N^{lat}|P^0`$ $`=`$ $`{\displaystyle \frac{8}{3}}\overline{P^0}|\overline{b}\gamma _\mu \gamma _5q|00|\overline{b}\gamma _\mu \gamma _5q|P^0{\displaystyle \frac{16}{3}}\overline{P^0}|\overline{b}\gamma _5q|00|\overline{b}\gamma _5q|P^0,`$ (31)   
$`\overline{P^0}|𝒪_S^{lat}|P^0`$ $`=`$ $`{\displaystyle \frac{5}{3}}\overline{P^0}|\overline{b}\gamma _5q|00|\overline{b}\gamma _5q|P^0,`$ (32)   
where $`|P^0`$ denotes a heavy-light pseudo-scalar meson at rest, and $`0|\overline{b}\gamma _\mu \gamma _5q|P^0`$ is related to the pseudo-scalar decay constant

$$\overline{P^0}|A_0^{lat}(1/a)|0=0|A_0^{lat}(1/a)|P^0=f_PM_P.$$
(33)   
Let us now consider a decomposition of the $`b`$-quark field $`\overline{b}`$ into the two-component non-relativistic quark $`Q^{}`$ and anti-quark $`\chi `$ fields. Up to $`O(1/m_Q^2)`$ we have

$`\overline{b}\gamma _5q`$ $`=`$ $`(Q^{}0)\left(1+{\displaystyle \frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}}\right)\gamma _5q(0\chi )\left(1+{\displaystyle \frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}}\right)\gamma _5q,`$ (34)   
$`\overline{b}\gamma _0\gamma _5q`$ $`=`$ $`(Q^{}0)\left(1{\displaystyle \frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}}\right)\gamma _5q+(0\chi )\left(1{\displaystyle \frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}}\right)\gamma _5q,`$ (35)   
and then

$`\overline{P^0}|\overline{b}\gamma _5q|0`$ $`=`$ $`\overline{P^0}|(Q^{}0)\left(1+{\displaystyle \frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}}\right)\gamma _5q|0`$ (36)   
$`=`$ $`\overline{P^0}|\overline{b}\gamma _0\gamma _5q|0+2\overline{P^0}|(Q^{}0){\displaystyle \frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}}\gamma _5q|0,`$ (37)   
$`0|\overline{b}\gamma _5q|P^0`$ $`=`$ $`0|(0\chi )\left(1+{\displaystyle \frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}}\right)\gamma _5q|P^0`$ (38)   
$`=`$ $`0|\overline{b}\gamma _0\gamma _5q|P^0+20|(0\chi ){\displaystyle \frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}}\gamma _5q|P^0.`$ (39)   
By defining $`\delta f_P`$ as

$$\overline{P^0}|(Q^{}0)\frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}\gamma _5q|0=0|(0\chi )\frac{\stackrel{}{\gamma }\stackrel{}{D}}{2m_Q}\gamma _5q|P^0\delta f_PM_P,$$
(40)   
we obtain

$`\overline{P^0}|𝒪_N^{lat}|P^0`$ $`=`$ $`{\displaystyle \frac{8}{3}}f_P^2M_P^2\left(18{\displaystyle \frac{\delta f_P}{f_P}}\right),`$ (41)   
$`\overline{P^0}|𝒪_S^{lat}|P^0`$ $`=`$ $`{\displaystyle \frac{5}{3}}f_P^2M_P^2\left(14{\displaystyle \frac{\delta f_P}{f_P}}\right).`$ (42)   
In our previous work  we denoted $`\delta f_P`$ as $`\delta f_P^{(2)}`$.

Thus the VSA for $`B_N^{lat}`$ and for $`B_S^{lat}`$ read

$`B_N^{lat(\mathrm{VSA})}`$ $`=`$ $`18{\displaystyle \frac{\delta f_P}{f_P}},`$ (43)   
$`B_S^{lat(\mathrm{VSA})}`$ $`=`$ $`{\displaystyle \frac{5}{8}}\left[14{\displaystyle \frac{\delta f_P}{f_P}}\right],`$ (44)   
neglecting the higher order contribution of order $`1/m_Q^2`$. Results for $`\delta f_P/f_P`$ is available at $`\beta `$=5.8 in Ref. . We plot them in Fig. 5(b) and 5(c) by crosses, which show a qualitative agreement with the measured values.

Di Pierro and Sachrajda  pointed out that the value of several $`B`$-parameter-like matrix elements of the $`B`$ meson is explained by the VSA surprisingly well in the static limit. Here we find that the $`1/m_Q`$ dependence of $`B_X^{lat}(1/a)`$ can also be reproduced qualitatively. This result suggests that the vacuum saturation is a reasonable qualitative picture for the heavy-light meson. It does not mean, however, that the VSA works quantitatively for $`B_B(\mu )`$, and careful lattice studies are necessary for precise calculation of the $`B`$-parameters.

## VI Remaining uncertainties and the final result

To estimate the systematic uncertainties in lattice calculations is a difficult task. In our case it is even more true, since we have a simulation result only at a single $`\beta `$ value. However we attempt to do it, giving a dimension counting of missing contributions.

The following sources of systematic errors are possible:

* discretization error: Both of the heavy and light quark actions are $`O(a)`$-improved at tree level, and there is no discretization error of $`O(a\mathrm{\Lambda }_{QCD})`$. The leading error is of $`O(a^2\mathrm{\Lambda }_{QCD}^2)`$ and of $`O(\alpha _sa\mathrm{\Lambda }_{QCD})`$. The second one is from the missing one-loop perturbative correction in the $`O(a)`$-improvement (its matching coefficient has been already obtained in Ref. ). We naively estimate the size of them to be $`O(a^2\mathrm{\Lambda }_{QCD}^2)O(a\mathrm{\Lambda }_{QCD}\alpha _s)`$ 5%, assuming $`1/a`$1.6 GeV, $`\mathrm{\Lambda }_{QCD}`$300 MeV and $`\alpha _s`$0.3.
* perturbative error: The operator matching of the continuum and lattice $`\mathrm{\Delta }B`$=2 operators are done at one-loop level. Thus the $`O(\alpha _s^2)`$ correction is another source of error. In addition, we use the one-loop coefficient for the static lattice action, though our simulation has been done with the NRQCD action. The error in this mismatch is as large as $`O(\alpha _s/(am_Q))`$ and $`O(\alpha _s\mathrm{\Lambda }_{QCD}/m_Q)`$. The size of these contributions is $`O(\alpha _s^2)`$ $``$ $`O(\alpha _s/(am_Q))`$ $``$ 10% and $`O(\alpha _s\mathrm{\Lambda }_{QCD}/m_Q)`$ $``$ 2%.
* relativistic error: Since we have performed a set of simulations with the $`O(1/m_Q)`$ action and the $`O(1/m_Q^2)`$ action, we can estimate the error in the truncation of the non-relativistic expansion. As we have shown, the difference between the results with $`O(1/m_Q)`$ and $`O(1/m_Q^2)`$ is small ($``$ 2%) around the $`B`$ meson mass. Then the higher order ($`O(1/m_Q^3)`$) effect is negligible.
* chiral extrapolation: We have only two light quark $`\kappa `$ values. Then the linear behavior in the chiral extrapolation is nothing but an assumption. Although the light quark mass dependence is small and the assumption is a reasonable one, we conservatively estimate the error from the difference between the data at our lightest $`\kappa `$ value ($`\kappa `$=0.1365) and $`\kappa _c`$. It leads 3% for $`B_{B_d}(m_b)`$.
* quenching error: All results are obtained in the quenched approximation. Study of the sea quark effect is left for future work.

Taking them into account, we obtain the following values as our final results from the quenched lattice,

$`B_{B_d}(m_b)`$ $`=`$ $`0.75(3)(12),`$ (45)   
$`{\displaystyle \frac{B_{B_s}}{B_{B_d}}}`$ $`=`$ $`1.01(1)(3),`$ (46)   
where the first error is statistical and the second a sum of all systematic errors in quadrature. In estimating the error of the ratio $`B_{B_s}/B_{B_d}`$ we consider the error from chiral extrapolation only, assuming that other uncertainties cancel in the ratio. The above result is related to the scale invariant $`B`$-parameter $`\widehat{B}_{B_d}`$ as

$$\widehat{B}_{B_d}=\{\begin{array}{c}\left[\alpha _s(m_b)\right]^{6/23}B_{B_d}(m_b)=1.12(4)(18)\hfill \\ \left[\alpha _s(m_b)\right]^{6/23}\left[1+\frac{\alpha _s(m_b)}{4\pi }J_5\right]B_{B_d}(m_b)=1.15(5)(18)\hfill \end{array},$$
(47)   
using the leading and next-to-leading formula, respectively, where we use $`\mathrm{\Lambda }_{QCD}^{(5)}`$=0.237 GeV and the two-loop $`\beta `$-function.

## VII Conclusion

In this paper we investigate the $`O(\mathrm{\Lambda }_{QCD}/m_Q)`$ and $`O(\mathrm{\Lambda }_{QCD}^2/m_Q^2)`$ effects on the $`B`$-parameter. We find that there is no significant mass dependence in the leading operator contribution $`B_L^{lat}(1/a)`$, while the mixing contributions $`B_N^{lat}(1/a)`$ and $`B_S^{lat}(1/a)`$ have large $`O(\mathrm{\Lambda }_{QCD}/m_Q)`$ corrections. The $`O(\mathrm{\Lambda }_{QCD}^2/m_Q^2)`$ correction for each $`B_X^{lat}(1/a)`$ is, however, reasonably small for the $`B`$ meson as we naively expected. The observed $`1/m_Q`$ dependence is qualitatively understood using the vacuum saturation approximation for the lattice matrix elements.

The lattice NRQCD calculation predicts the small but non-zero negative slope in the mass dependence of $`B_P(m_b)`$ and about 10% reduction from static limit to the physical $`B`$ meson. In the present analysis, we combine lattice simulation for finite heavy quark mass with the mass independent matching coefficients determined in the static limit. The dominant uncertainty is, therefore, arising from the finite mass effects in the perturbative matching coefficients. For more complete understanding of the $`1/m_Q`$ dependence, matching coefficients with the finite heavy quark mass are necessary.

## Acknowledgment

We would like to thank S. Tominaga for useful discussion. Numerical calculations have been done on Paragon XP/S at Institute for Numerical Simulations and Applied Mathematics in Hiroshima University. We are grateful to S. Hioki for allowing us to use his program to generate gauge configurations. T.O. is supported by the Grants-in-Aid of the Ministry of Education (No. 10740125). H.M. would like to thank the JSPS for Young Scientists for a research fellowship.

## A

In this appendix, we compare our perturbative matching by simple one-loop formula with the renormalization group (RG) improved ones used in the previous static calculations . Since the matching procedure for determining the RG improved coefficients is given in Refs.  in detail (see also Refs. ), we just show the results appropriate for our actions and definition of operator.

Considering the matching of a ratio $`B_B(m_b)=𝒪_L(m_b)/(8/3)A_0^2`$ again as in section II, the RG improved versions of $`Z_{X/A^2}(m_b;1/a)`$ are as follows,

$`Z_{L/A^2}(m_b;1/a)`$ $`=`$ $`Z_L^{cont}\left(1+{\displaystyle \frac{\alpha _s}{4\pi }}(D_L2D_A)\right),`$ (A1)   
$`Z_{R/A^2}(m_b;1/a)`$ $`=`$ $`Z_L^{cont}\times {\displaystyle \frac{\alpha _s}{4\pi }}D_R,`$ (A2)   
$`Z_{N/A^2}(m_b;1/a)`$ $`=`$ $`Z_L^{cont}\times {\displaystyle \frac{\alpha _s}{4\pi }}D_N,`$ (A3)   
$`Z_{S/A^2}(m_b;1/a)`$ $`=`$ $`Z_S^{cont},`$ (A4)   
and

$`Z_L^{cont}`$ $`=`$ $`\left(1+{\displaystyle \frac{\alpha _s(m_b)}{4\pi }}({\displaystyle \frac{26}{3}})\right)\left(1+{\displaystyle \frac{\alpha _s(1/a)\alpha _s(m_b)}{4\pi }}(0.043)\right)`$ (A6)   
$`+{\displaystyle \frac{\alpha _s(m_b)}{4\pi }}(8)\left(\left({\displaystyle \frac{\alpha _s(m_b)}{\alpha _s(1/a)}}\right)^{8/25}1\right){\displaystyle \frac{1}{4}},`$   
$`Z_S^{cont}`$ $`=`$ $`{\displaystyle \frac{\alpha _s(m_b)}{4\pi }}(8)\left({\displaystyle \frac{\alpha _s(m_b)}{\alpha _s(1/a)}}\right)^{8/25}.`$ (A7)   
Numerical values of $`Z_{\{L,R,N,S\}/A^2}(m_b;1/a)`$ are given in Table I together with those by the simple one-loop formula, where we use the V-scheme coupling  as $`\alpha _s`$ appearing in Eqs. (A1)-(A4) while the couplings in Eqs. (A6) and (A7) are defined in the continuum $`\overline{\mathrm{MS}}`$ scheme with $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{(4)}`$=344 MeV, which corresponds to $`\mathrm{\Lambda }_{\overline{\mathrm{MS}}}^{(5)}`$= 237 MeV.

Now assuming each $`B_X^{lat}`$ is of $`O(1)`$, the dominant effects of resummation arise from $`Z_L/A^2`$ and $`Z_S/A^2`$. Since, however, its difference is at most 5% level and the effects from $`Z_L/A^2`$ and $`Z_S/A^2`$ are destructive, the total effect amount to less than 3%. To be specific, using our data extrapolated to the static limit (see Table IV) to calculate $`B_B^{stat}(m_b)`$, we obtain the results tabulated in Table V from the two matching procedures. In this case the effect of resummation is almost negligible.