Precision lattice QCD computation of the $B^*B\pi$ coupling

The static $B^{*}B\pi$ coupling, $\hat{g}_\chi$, a low energy constant in the leading order heavy meson chiral Lagrangian, is determined using $N_\mathrm{f} = 2$ lattice QCD. We use CLS ensembles with lattice spacings and pion masses down to $a = 0.05\mathrm{fm}$ and $m_{\pi}=270\mathrm{MeV}$, and perform combined continuum and chiral extrapolations of our results which have a much better accuracy than previous numbers in the literature. As a by-product, we determine the coupling between the first radial excitations in the $B$ and $B^{*}$ channels ($\hat{g}_{22}$). Accounting for all uncertainties, which are dominated by the chiral extrapolation, we obtain $\hat{g}_\chi = 0.492(29)$, while $\hat{g}_{22}$ is somewhat smaller. The comparison to a precise quenched computation suggests that there is little influence by the sea quarks and $\hat{g}_\chi$ will not change much when a dynamical strange quark is included.


Introduction
Low energy QCD is described by an effective theory based on spontaneously broken global SU (N f ) L × SU (N f ) R chiral symmetry, where N f is the number of light quark flavors. At the same time, a low energy expansion of hadrons with a single heavy quark with mass m h exists and is known as Heavy Quark Effective Theory [1][2][3][4]. These chiral and heavy quark expansions may be combined to construct effective theories for the low-energy dynamics of hadrons containing a single heavy quark [5][6][7].
The theory that describes mesons is called Heavy Meson Chiral Perturbation Theory (HMχPT) and contains a single additional leading low energy constant with respect to standard χPT. This additional low energy constant,ĝ χ , describes the coupling of heavy mesons to pseudo-Goldstone bosons in the chiral (m 2 π → 0) and static (m h → ∞) limits. The couplingĝ χ is relevant for the computation of B-physics matrix elements from lattice QCD, exemplified by the ALPHA collaboration HQET program [8][9][10][11][12]; it enters in chiral extrapolations of hadronic parameters needed for heavy flavor phenomenology, such as the B-meson decay constant and the B-meson semi-leptonic decay form factors. This coupling is also referred to as the B * Bπ coupling, where the pseudoscalar and vector static-light mesons are denoted B and B * . Note that the static (m h → ∞) limit is implied.
One way to determineĝ χ is through phenomenological fits to experimental data. A determination from D * → Dπ decays [13] yields a value ofĝ = 0.61 (6). However, this extraction is affected by O(1/m c ) and O(m 2 π ) errors, where especially the first ones are hard to estimate. Unfortunately, the process B * → Bπ is kinematically forbidden, complicating the estimation of the O(1/m h ) errors from experimental results. For a recent review of results, including also quark model and QCD sum rules calculations, see Ref. [14].
In this work we employ a different approach. Using lattice QCD simulations, we calculate a matrix element of the (light-light) axial current in QCD which is equivalent toĝ χ in leading order HMχPT. Namely, we compute (ignoring renormalization and improvement in this introduction) where ψ d (ψ u ) annihilates a down(up) quark and the index k = 1, 2, 3 is not summed over. We use the finite volume normalization of states B 0 (p)|B 0 (p) = B * (p)|B * (p) = 2L 3 = 2V , where L is the linear size of the simulated torus. We work directly in the static limit for the heavy quark, but at finite light quark mass. Thereforeĝ χ is eventually obtained by an extrapolation of our results forĝ to the zero light quark mass (chiral) limit as well as the a → 0 (continuum) limit, where a is the lattice spacing of our simulations. There have been previous determinations ofĝ χ using lattice QCD with N f = 0, 2, 3 dynamical light quark flavors directly in the static limit [15][16][17][18], as well as determinations at the charm [19] and bottom [20] points. However, the lattice spacing dependence of this quantity has not yet been thoroughly investigated. In this work we perform a continuum extrapolation in both the N f = 0 and N f = 2 theories and find small lattice spacing effects for our O(a) improved discretisation.
In Fig. 1 we compare our results (before extrapolations) to recent lattice results from Refs. [16][17][18]. We observe that a new quality is reached, reducing previous uncertainties   [16], Becirevic et al. [17] and Detmold et al. [18]. For Ref. [18], which employs N f = 2 + 1 dynamical flavors, we take the results for a single level of link smearing in the static action.
by an order of magnitude in the region of interest, namely at small pion masses 1 . This is achieved by both improved techniques [21] and good statistics.
Additionally, we quote results for the matrix element of the first radial excitations. To this end we defineĝ mn (y, a) = where m, n = 1 are the ground states of B and B * mesons while m, n > 1 refer to their excitations. Apart from our main object of study,ĝ 11 =ĝ, we quote rough numbers for g 22 . Preliminary results forĝ 11 andĝ 22 have appeared in Ref. [22] (together with results for g 12 ) and a preliminary account of our present work can be found in Ref. [23]. The variable y is proportional to the square of the pion mass and will be defined when we discuss the chiral and continuum limit to arrive atĝ χ ≡ĝ 11 (0, 0) andĝ χ 22 ≡ĝ 22 (0, 0). In Sec. 2 we describe our techniques. In Sec. 3.1 we show results in the quenched approximation, where a continuum limit is taken for bothĝ 11 andĝ 22 at a fixed quark mass m q ≈ m strange . In Sec. 3.2 we discuss the N f = 2 results, in particular the chiral and continuum extrapolations. Finally we conclude in Sec. 4.

Methodology
Here we describe some details of the lattice calculation ofĝ mn , namely the proper definition of the axial current and the technology to obtain precise matrix elements from correlation functions including the estimation of systematic errors due to excited state contributions.
We also detail our stochastic technique utilising translation invariance. The ensembles used in the numerical application are explained in Secs. 3.1 and 3.2.

Discretisation and renormalization
We employ both the HYP1 and HYP2 discretizations of the static quark action [24,25] to mitigate the signal-to-noise problem and provide a further check on discretization effects. Generally, results from these two discretizations are compatible within statistical errors, but they are also strongly correlated. We will thus show both of them in the tables, but only use HYP2 in the detailed analysis.
The light quarks are non-perturbatively O(a) improved Wilson quarks [26][27][28] and the improved and renormalized axial current is where A k has exactly the form given before and P (x) = ψ d (x)γ 5 ψ u (x) is the appropriate pseudoscalar density. For the required values of the bare coupling, the renormalization constant Z A is known non-perturbatively for both N f = 0 [29] and N f = 2 [30,31] while for the improvement coefficient b A we use the expansion in the bare coupling g 2 0 to first order with the one-loop coefficient of Ref. [32]. For our zero momentum transfer matrix element g, the ∂ k P term vanishes identically; c A is not needed.

Matrix elements from the GEVP
The matrix elementsĝ mn are accessible in lattice QCD via three-point correlation functions in Euclidean time, which (from the transfer matrix formalism) have the following representation whereÔ i andÔ k j are suitable interpolating fields for the B 0 and B * + k mesons (respectively), ψ im = 0|Ô i |B, m = 0|Ô k i |B * k , m and E m is the energy of the mth state. In the static limit the B and B * energy levels are the same and furthermore C 3pt ij is independent of k as indicated by our notation.
To isolate the desired matrix elements, we also require the two-point correlation functions Rather than analyzing C 3pt (t 1 , t 2 ) directly, we employ where the position of the current insertion is summed over [23,33,34]. The use of this summed correlation function improves the convergence in t but was proposed in Ref. [33] for different reasons.
In order to extract the desired matrix elements we choose a set of N interpolating operators and form the N × N correlation matrices C 2pt ij (t) and D 3pt ij (t). We then employ solutions of a generalized eigenvalue problem (GEVP) [35][36][37] to accelerate the asymptotic (in t) behavior and enable the extraction ofĝ nn for n > 1.
It has been proven recently that the GEVP may be combined with summed insertions [21,22] to achieve a further reduction in the contribution from excited states. It was demonstrated that the summed insertion is particularly advantageous in the extraction of excited state matrix elements when compared to the ordinary GEVP. For completeness, we review the main points. We begin by solving the following GEVP where t 0 ≥ t/2. It can be shown [21] that where ∆ N,n = E N +1 − E n and (., .) denotes an inner product over the GEVP indices. The important result is that (asymptotically) the corrections fall exponentially in ∆ N,n t. The large energy gap ∆ N,n , which in our application is above 1 GeV, is a virtue of the GEVP and the factor t is due to the summed insertion. Without summation, t would be replaced by min(t 1 , t 2 ), see eq. (2.2). We also note that the GEVP renders excited state matrix elements accessible. In our numerical application, the interpolating fields are constructed from Gaussian smearing operators where ∆ is the gauge-covariant spatial Laplace operator with APE-smeared links. The approximate width r i ≈ 2a √ κ G R i is chosen to keep the smearing radii at r i ≈ 0.2, 0.3, 0.7 fm for each lattice spacing. More details about the construction of these wave-functions can be found in Refs. [8,10].
From the correlation functions we construct M eff n (t) ≡ M eff n (t, t − a) and examine the large t behaviour. Beginning with t ≈ r 0 ≈ 0.5fm, we increase t until the asymptotic corrections due to excited states are small enough so that where δt = 1 ∆ N,n and σ(t) is the statistical (1σ) error on M eff n (t). We call the first t at which this condition is satisfied t min . Under the assumption that the asymptotic decay eq. (2.7) has roughly set in at t min , our requirement of eq. (2.9) means that statistical errors exceed systematic ones by a factor e − 1 ≈ 2 at t = t min . We then define our estimate ofĝ nn as the weighted average of M eff (t) over the range [t min , t max ] with t max chosen to avoid points with excessive statistical errors. The estimate for the statistical error on M eff (t) will be discussed in future sections, in particular for the N f = 2 results where autocorrelations must be treated with care.

Use of random sources
In order to reduce statistical fluctuations, we use full translation invariance everywhere. We achieve this by a stochastic estimation of one of the spatial sums, using [38] a random U (1) source on each time-slice of the lattice ('time-dilution' [39]) and a 'sequential inversion' (eq. (2.13) below) for the insertion of the axial current.
The explicit expression for the two-point function is where η r (y; t 0 ) ∝ δ(y 0 − t 0 ) ∈ U (1) is a random U(1) field on timeslice t 0 and vanishes otherwise, S b (x, y) is the easily computed static quark propagator of the b-quark and D is the Dirac operator of the light quarks. The subscript r enumerates the random source fields. Similarly, the three-point functions with summed current insertions read In the numerical application, we average over r = 1 . . . N r , all values 0 ≤ t 0 ≤ T − a and over k = 1, 2, 3. Hence, eq. (2.11) and eq. (2.13) add up to N r × T /a × 4 Dirac equations which need to be solved on each gauge configuration. The ensemble average . indicates an average over the random fields η r as well as the gauge fields.

Results
We now discuss numerical results for N f = 0 and N f = 2. In both cases we employ the usual periodic and anti-periodic temporal boundary conditions for the gauge and fermion fields, respectively.

N f = 0 continuum limit
We first apply the methods discussed in Sec. 2 to a set of three ensembles of quenched gauge configurations used previously in the ALPHA collaboration HQET program [9], with the goal of taking the continuum limit ofĝ 11 andĝ 22 . Details of the ensembles and measurements are given in Tab. 1. The valence quark mass on each of these ensembles was tuned to reproduce the physical strange quark mass [40]. The effective matrix elements M eff n (t), n = 1, 2 for the L/a = 16 ensemble are shown together with their plateau averages in the left plot of Fig. 2. As autocorrelations play no role here, the errors are estimated using 100 single-elimination jackknife bins after first averaging over the N r sources on each configuration. The range of the plateau averages is chosen according to the criteria in Sec. 2.
Finally, we perform continuum extrapolations of the renormalized axial current matrix elements. These extrapolations forĝ 11 andĝ 22 are shown in the right plot of Fig. 2 and suggest that cutoff effects are small for these quantities. Simple constant extrapolation of   as our final quenched results, where a linear term in a 2 is allowed in the fit formula. It should be noted that our result forĝ 11 is compatible with the previous result of Ref. [23], but utilizes a more robust treatment of the systematic errors due to excited states, namely that of Sec. 2.  Details of the N f = 2 ensembles. The pseudoscalar meson masses and lattice spacings are taken from Ref. [31] and N r is as defined in Sec. 2. We also list our estimate of the exponential autocorrelation time τ exp in units of the separation between configurations.

N f = 2 results
We next use ensembles of the Coordinated Lattice Simulations (CLS) community effort. Details are tabulated in Tab. 2. While we follow the same procedure to calculate the bare matrix elements, the large autocorrelations present in HMC simulations with periodic boundary conditions must be taken into account in order to safely estimate the statistical errors. To this end, we follow the procedure of Ref. [42] and attach an exponential 'tail' to the autocorrelation functions of the matrix elements, with a fall-off ∼ exp(−t MC /τ exp ), where τ exp has been estimated roughly in Ref. [42].
A selection of the effective matrix elements is shown in Fig. 3. This figure also compares the use of the summed insertion with and without the GEVP for M eff 1 (t). A clear picture for the difference is not easily seen for the few points available. However, the figure indicates that while corrections are not necessarily smaller with the GEVP at small times t ≈ r 0 /2, the approach to the plateau is then accelerated soon after. Indeed, this is needed for our criteria eq. (2.9) to apply.
For the excited state matrix element our statistical errors are not small enough to apply eq. (2.9) to fix the start of the plateaux. Here we simply inspect the figures and choose a fixed t min ≈ 0.5 fm in physical units. The large statistical errors seem to dominate over the systematic errors due to excited state contributions. The renormalized matrix elements together with the statistical errors estimated using the additional exponential tail are collected in Tab. 3. . The left column also shows the effective matrix elements obtained without the GEVP represented by the gray points as in Fig. 2   Next, the combined chiral and continuum extrapolation is performed. Let us first concentrate on the ground state matrix element, which yieldsĝ χ . We parameterize the quark mass dependence by the pion mass (as in Ref. [8]) through the variable 2 y = m 2 π /(8π 2 f 2 π ). As already evident from Fig. 1, the data is rather linear in m 2 π (or y), while chiral perturbation theory predicts a significant logarithmic modification [18,43]. We therefore perform two extrapolations to the chiral limit. Namely we fit to the two formŝ g lin 11 (y, a) =ĝ χ + By + Ca 2 , (3.2)

11
(y, a) =ĝ χ 1 − (1 + 2(ĝ χ ) 2 )y log y + By + Ca 2 , (3.3) where B, C, and the desiredĝ χ are fit parameters. In both forms, the Ca 2 term can probe cutoff effects and we consider terms of order ya 2 as too small to be relevant, as we do with a 3 and y 2 . In both fits, the results for C are compatible with zero. As our central values we take the fit results forĝ χ with C set to zero. The linear fit, eq. (3.2), then yieldŝ g χ = 0.513(8) while the correct asymptotic form, eq. (3.3), extrapolates further down tô g χ = 0.469 (7). We combine these numbers to our central result g χ = 0.492 (29) .
The error is by far dominated by the difference of the two chiral extrapolations. We have chosen a range which encompasses the linear and the NLO extrapolation and their errors. Allowing for non-vanishing C changes rather little concerning this result. Of course the situation is far from perfect: the theoretically well motivated functional form is not verified by the data; a linear dependence fits somewhat better. However, in the end we are interested in the extrapolated value and it seems very safe to assume that it lies in the range eq. (3.4). Indeed, if NLO chiral behaviour sets in at masses which are below the ones in Fig. 4, the downward bend will happen later and the result will be in between the two values shown in the figure and used to form eq. (3.4).
Forĝ 22 the functional form including chiral log's is not known. Also the data are much less precise. We thus perform a simple linear extrapolation both with and without an a 2 term. They are shown on the right side of Fig. 4. A range covering both results iŝ g 22 = 0.425(70). (3.5)

Conclusions
In this paper we have presented a precise N f = 2 determination ofĝ χ , the leading low energy constant appearing in HMχPT parametrizing the coupling of heavy-light mesons to pions. We have calculated the bare matrix elements using solutions of the GEVP together with the summed insertion technique, resulting in a precision which exceeds previous ones by an order of magnitude. We renormalized these matrix elements non-perturbatively.
We have taken the continuum and chiral limits assuming both a phenomenological linear behaviour in the square of the pion mass and next-to-leading-order continuum HMχPT. Two discretizations of the static quark action serve as a further check on lattice spacing effects. These two discretizations give statistically compatible results for all quantities. Our central result is eq.  [20]. Within the cited overall uncertainties all previous numbers are in agreement with our more precise value. In fact the agreement is better than one might have expected given that some numbers come from extrapolations from rather large pion masses and lattice spacings.
As discussed in Sec. 2, we have treated the systematic errors due to excited states in a conservative manner. Similarly, we have also safely estimated the statistical errors by including tails in the autocorrelation functions. At our finest lattice spacing, these are significant, see Sec. 3.2. However, in the end, the dominating uncertainty comes from the fact that the data does not appear to be at such small pion masses where NLO chiral behaviour can be seen. Instead an approximately linear behaviour in m 2 π prevails down to m π = 270 MeV. We thus take a final range which also covers the result of a simple linear extrapolation.
In comparison to the quenched result, we have to take into account that eq. (3.1) is for a light quark mass set to the strange mass. This corresponds to y ≈ 0.2, outside the range of Fig. 4. The figure then suggests that N f = 2 and the quenched number agree within at least 5% precision. We do not see any sea quark effects at the strange mass. It thus appears safe to use eq. (3.4) with its more than 5% error also for the three (or more) flavor theory.
Despite our limited control over the chiral limit, our determination ofĝ χ is precise enough to help the chiral extrapolation of many quantities of phenomenological interest in heavy meson physics. For example, this result is used broadly in the ALPHA collaboration HQET program. the N6 ensemble from Ref. [22], and Patrick Fritzsch for useful comments on an earlier version of this manuscript. This work is supported by the Deutsche Forschungsgemeinschaft in the SFB/TR 09 and by the European community through EU Contract No. MRTN-CT-2006-035482, "FLAVIAnet". We are grateful to NIC and to the Norddeutsche Rechnerverbund for allocating computing resources to this project. Some of the correlation function measurements were performed on the PAX cluster at DESY, Zeuthen.