Fermilab Lattice and MILC Collaborations

Motivation — Hadron decays that proceed through flavor-changing neutral currents may be sensitive to new physics, because their leading Standard-Model contributions are loop suppressed. Here we study the semileptonic decay $B\to\pi\ell^{+}\ell^{-}$, which proceeds through a $b\to d$ transition. Hadronic effects in this decay are parametrized by three form factors. In this Letter, we present the first ab-initio QCD calculation of the tensor form factor $f_{T}$, based on lattice-QCD work that also yielded the vector and scalar form factors, $f_{+}$ and $f_{0}$ Bailey et al. (2015a). Lattice QCD has several advantages over other approaches to the form factors Ball and Zwicky (2005); Wang et al. (2008); Duplančić et al. (2008); Wu and Huang (2009); Faustov and Galkin (2014); Ali et al. (2014); Li et al. (2014); Hambrock et al. (2015), particularly in providing a path to controlled uncertainties that can be systematically reduced Kronfeld (2012).

The LHCb experiment recently made the first observation of $B^{+}\to\pi^{+}\mu^{+}\mu^{-}$ Aaij et al. (2012), while the $B$-factories have set limits on the $e^{+}e^{-}$ and $\tau^{+}\tau^{-}$ channels Wei et al. (2008); Lees et al. (2013); Lutz et al. (2013). Below we present the first calculations of $B\to\pi\ell^{+}\ell^{-}$ ($\ell=e,\mu,\tau$) observables in the Standard Model using form factors with fully controlled uncertainties.

The form factors $f_{+}$, $f_{0}$, and $f_{T}$ suffice to parameterize $B\to\pi$ decays in all extensions of the Standard Model. New physics from heavy particles—such as those appearing in models with supersymmetry Bobeth et al. (2001); Demir et al. (2002); Choudhury and Gaur (2002); Wang et al. (2008), a fourth generation Hou et al. (2013), or extended Aliev and Savci (1999); Iltan (1999); Bobeth et al. (2001); Erkol and Turan (2002); Erkol et al. (2005); Song et al. (2008) or composite Gripaios et al. (2015) Higgs sectors—alter Wilson coefficients in the effective Hamiltonian pertaining to particle physics below the electroweak scale Grinstein et al. (1989); Buras et al. (1994); Huber et al. (2006); Altmannshofer et al. (2009). Whatever these unknown particles may be, the hadronic physics remains the same.

Lattice-QCD calculation — Our work on $f_{T}(q^{2})$ was carried out in parallel with $f_{+}(q^{2})$ and $f_{0}(q^{2})$. Our aim in Ref. Bailey et al. (2015a) was a precise determination of the Cabibbo-Kobayashi-Maskawa (CKM) element $|V_{ub}|$, and every step of the analysis was subjected to many tests. Further, two of the authors applied a multiplicative offset to the numerical data at an early stage. This “blinding” factor was disclosed to the others only after finalizing the error analysis. Full details of the simulation parameters, analysis, and cross-checks are given in Ref. Bailey et al. (2015a).

Our calculation uses ensembles of lattice gauge-field configurations Aubin et al. (2015a); *asqtad:en06b; *asqtad:en05a; *asqtad:en05b; *asqtad:en04a; *asqtad:en15a; *asqtad:en15b; *asqtad:en14a; *asqtad:en13a; *asqtad:en13b; *asqtad:en12a; *asqtad:en23a; *asqtad:en23b; *asqtad:en20a; *asqtad:en20b; *asqtad:en19a; *asqtad:en18a; *asqtad:en18b; *asqtad:en24a from the MILC Collaboration Bernard et al. (2001); Aubin et al. (2004); Bazavov et al. (2010), which are generated with a realistic sea of up, down, and strange quarks. In practice, the up and down sea quarks have the same mass, and the strange-quark mass is tuned close to its physical value. The statistics are high, with 600–2200 gauge-field configurations per ensemble. The physical volume is large enough that we can repeat the calculation in different parts of the lattice, thereby quadrupling the statistics. We use four lattice spacings ranging from 0.12 fm to 0.045 fm to control the extrapolation to zero lattice spacing.

The tensor form factor is defined via the matrix element of the $b\to d$ tensor current $i\bar{d}\sigma^{\mu\nu}b$:

where $p_{B}$ and $p_{\pi}$ are the particles’ momenta and $q=p_{B}-p_{\pi}$ is the momentum carried off by the leptons. The Lorentz invariant $q^{2}$ is related to the pion energy in the $B$-meson rest frame via $E_{\pi}=(M_{B}^{2}+M_{\pi}^{2}-q^{2})/2M_{B}$. In the finite volume that can be simulated on a computer, $E_{\pi}$ takes discrete values, dictated by the spatial momenta $\bm{p}_{\pi}$ compatible with periodic boundary conditions. Because statistical and discretization errors increase with pion momentum, we restrict $|\bm{p}_{\pi}|\leq|2\pi(1,1,1)/L|$. The resulting simulation range of $E_{\pi}\lesssim 1$ GeV is significantly smaller than the kinematically allowed range of $E_{\pi}\leq 2.5$ GeV. Extending this discrete set of calculations into the full $q^{2}$ dependence is the central challenge of this work, and is met in two steps.

The two light quarks (up and down) have a mass larger than it should be, but the range simulated is wide and the smallest pion mass is 175 MeV, close to Nature’s 140 MeV. Therefore, we can apply an effective field theory of pions—chiral perturbation theory—to extrapolate the simulation data to the physical point. We use a form of chiral perturbation theory adapted to lattice QCD, with additional terms describing the lattice-spacing dependence Aubin and Bernard (2007); Bailey et al. (2015b) and with modifications needed for energetic final-state pions Bijnens and Jemos (2010). As discussed in Ref. Bailey et al. (2015a), we try several fit variations. For example, we replace the loop integrals with momentum sums appropriate to the finite volume, finding negligible changes in the results. Our final fit includes next-to-next-to-next-to-leading order analytic terms and terms to model the discretization errors of the heavy quark. The latter come from an effective field theory for heavy $b$ quarks Kronfeld (2000); Harada et al. (2002); Oktay and Kronfeld (2008).

Figure 1 shows the $q^{2}$ dependence of the errors after the chiral-continuum extrapolation just described.

Table 1 gives a numerical error budget for $f_{T}(q^{2}=20~\text{GeV}^{2})$.

The largest uncertainty comes from the statistical errors, as increased during the chiral-continuum extrapolation. This error is under good control for $q^{2}$ corresponding to the spatial momenta that we simulate, but grows large elsewhere.

The subdominant errors are as follows. To convert from lattice units to physical units, we introduce a physical distance $r_{1}$, which is defined via the force between static quarks Bernard et al. (2000); Sommer (1994). We use it to form physical, dimensionless quantities, which are the input data for the chiral-continuum fit. At the end, we set $r_{1}=0.3117\pm 0.0022$ fm Bazavov et al. (2012) based on a related lattice-QCD calculation of $r_{1}f_{\pi}$ Bazavov et al. (2009) and the pion decay constant $f_{\pi}=130.41$ MeV Olive et al. (2014). To propagate the parametric uncertainty in $r_{1}$ to $f_{T}$, we repeat the fit shifting $r_{1}$ by $\pm 1\sigma_{r_{1}}$, leading to the second line in Table 1.

In lattice gauge theory, the tensor current does not have the normalization used in QCD phenomenology. We obtain most of the normalization nonperturbatively El-Khadra et al. (2001) from $b\to b$ and $d\to d$ transitions with the vector current, with statistical errors below 1%. Another matching factor $\rho_{T}$ remains, but, by design and in practice, it is close to unity. We calculate $\rho_{T}$ at the renormalization scale $\mu=m_{b,\text{pole}}$ through first order in the QCD coupling $\alpha_{s}$. We estimate the resulting error of order $\alpha_{s}^{2}$ after removing a logarithmic dependence on the matching scale $\mu$, which is present in continuum QCD too. We then examine how the one-loop coefficient depends on heavy-quark mass, identifying the largest value, $\rho^{[1]}_{T,\,\text{max}}$. Finally, we estimate the error in $\rho_{T}$ to be $2\alpha_{s}^{2}|\rho^{[1]}_{T,\,\text{max}}|$, evaluating $\alpha_{s}$ on the second-finest lattice with $a\approx 0.06$ fm. This yields the 2% perturbative-matching uncertainty in Table 1.

The last two uncertainties arise as follows. When generating data, we choose the simulation quark masses based on short runs and previous experience. The full analysis yields better estimates. To correct the simulation $b$-quark mass a posteriori, we recompute $f_{T}$ on one ensemble with two additional values of the bare $b$-quark mass. Using the slope from all three mass values, we interpolate the data for $f_{T}$ slightly from the production $b$-quark mass to the physical value. This leaves an error due to the uncertainty in the size of the $b$-quark mass correction. The details for $f_{T}$ are nearly identical to those for $f_{+}$ Bailey et al. (2015a), leading to the same estimate, 0.4%, for this error. The light-quark mass dependence is embedded in the chiral-continuum extrapolation, described above. The parametric uncertainty from the input light-quark mass Bazavov et al. (2010) is propagated to $f_{T}$ by repeating the fit with $\pm 1\sigma_{m_{q}}$ shifts to these parameters, and is given in the penultimate line of Table 1.

The final line in Table 1 and the upper edge of the stack in Fig. 1 represent the quadrature sum of the systematic uncertainties with the chiral-continuum fit error.

Extension to all $q^{2}$ — To extend $f_{T}$ in the chiral-continuum limit from the range of simulated lattice momenta to the full kinematic range, $0<q^{2}\leq(M_{B}-M_{\pi})^{2}$, with controlled errors, we use a method based on the analytic structure of the form factor.

In the complex $q^{2}$ plane, $f_{T}(q^{2})$ has a cut for timelike $q^{2}\geq t_{+}\equiv(M_{B}+M_{\pi})^{2}$ and a pole at $q^{2}=M_{B^{*}}^{2}$ but is analytic elsewhere. The variable

maps the whole $q^{2}$ plane into the unit disk, with the cut mapped to the boundary and the semileptonic region mapped to an interval on the real axis. Unitarity bounds then guarantee that an expansion of $f_{T}$ in $z$ (with the $B^{*}$ pole removed) converges for $|z|<1$ Bourrely et al. (1981); Boyd et al. (1995); Lellouch (1996); Boyd and Savage (1997). Following Bourrely, Caprini, and Lellouch (BCL) Bourrely et al. (2009), we factor out the $B^{*}$ pole and expand in $z$:

choosing $t_{0}=(M_{B}+M_{\pi})(\sqrt{M_{B}}-\sqrt{M_{\pi}})^{2}$ to minimize $|z|$ in the semileptonic region. Although Eq. (3) was derived for the vector form factor $f_{+}$, we use it for the tensor form factor $f_{T}$ because the two form factors are proportional to each other at leading order in the $1/m_{b}$ expansion.

We determine the $b_{n}^{T}$ with a functional method connecting the independent functions of the chiral-continuum fit with the first several powers of $z$ Bailey et al. (2015a). Our preferred fit uses $N_{z}=4$; adding higher-order terms in $z$ does not significantly change the central value. Table 2 presents our final result for $f_{T}$ as coefficients of the $N_{z}=4$ BCL $z$ fit and the correlation matrix between them, where the errors include statistical and all systematic uncertainties.

This information can be used to reconstruct $f_{T}(q^{2})$ over the full kinematic range. Table 2 also provides the (mostly statistical) correlations between $f_{T}$, $f_{+}$, and $f_{0}$. Figure 2 shows the extrapolation of $f_{T}$ to $q^{2}=0$.

Table 2 and Fig. 2 represent the first main result of this Letter.

Implications — The largest contribution in the Standard Model to the amplitude for $B\to\pi\ell^{+}\ell^{-}$ is proportional to the vector form factor. Assuming that new physics does not contribute significantly to the tree-level decay $B\to\pi\ell\nu$, one can use experimental measurements of this process to constrain the shape of $f_{+}(q^{2})$, especially at low $q^{2}$. In Ref. Bailey et al. (2015a), we obtain the CKM element $|V_{ub}|$ from a combined $z$ fit to our lattice-QCD results for $f_{+}$ and $f_{0}$ and measurements of $\tau_{B}d\Gamma(B\to\pi\ell\nu)/dq^{2}$ from BaBar del Amo Sanchez et al. (2011); Lees et al. (2012) and Belle Ha et al. (2011); Sibidanov et al. (2013). This joint fit also yields the most precise current determinations of $f_{+}$ and $f_{0}$. To enable them to be combined with the results for $f_{T}$ from Table 2, Table 3 provides the correlations between the $z$-expansion coefficients for all three form factors.

The correlations are small, because $f_{+}$ contains independent experimental information.

Using $f_{T}$ from this work and $f_{+}$ and $f_{0}$ just described, we show the Standard-Model partial branching fractions for $B\to\pi\ell^{+}\ell^{-}$ in Fig. 3.

Other ingredients are needed besides the form factors. We take the Wilson coefficients from Ref. Huber et al. (2006), the CKM elements from Ref. Charles et al. (2005), the meson masses and lifetimes from Ref. Olive et al. (2014), and the $b$- and $c$-quark masses from Ref. Ali et al. (2014). To calculate contributions that cannot be parameterized by the form factors, we employ QCD factorization at low $q^{2}$ Beneke et al. (1999, 2000); Beneke and Feldmann (2001); Asatrian et al. (2001); Beneke et al. (2001); Asatryan et al. (2002); Asatrian et al. (2004); Beneke et al. (2005); Bobeth et al. (2007) and an operator product expansion (OPE) in powers of $E_{\pi}/\sqrt{q^{2}}$ at large $q^{2}$ Grinstein and Pirjol (2002); Seidel (2004); Grinstein and Pirjol (2004); Greub et al. (2008); Bobeth et al. (2010); Beylich et al. (2011); Bobeth et al. (2011, 2012). Full details will be provided in Ref. Du et al. (2015).

Table 4 presents numerical predictions for selected $q^{2}$ bins. The last error in parenthesis contains effects of parametric uncertainties in $\alpha_{s}$, $m_{t}$, $m_{b}$, $m_{c}$; of missing power corrections, taking 10% of contributions not directly proportional to the form factors; and of violations of quark-hadron duality, estimated to be 2% at high-$q^{2}$ Beylich et al. (2011). At low $q^{2}$, the uncertainty predominantly stems from the form factors; at high $q^{2}$, the CKM elements $|V_{td}^{*}V_{tb}|$ and form factors each contribute similar errors. Figure 3 and Table 4 represent the second main result of this Letter.

In the regions $q^{2}\lesssim 1~\text{GeV}^{2}$ and $6~\text{GeV}^{2}\lesssim q^{2}\lesssim 14~\text{GeV}^{2}$, $u\bar{u}$ and $c\bar{c}$ resonances dominate the rate. To estimate the total BR, we simply disregard them and interpolate linearly in $q^{2}$ between the QCD-factorization result at $q^{2}\approx 8.5~\text{GeV}^{2}$ and the OPE result at $q^{2}\approx 13~\text{GeV}^{2}$. While this treatment does not yield the full branching ratio, it does enable a comparison with LHCb’s published result, $\text{BR}(B^{+}\to\pi^{+}\mu^{+}\mu^{-})=23(6)\times 10^{-9}$ Aaij et al. (2012), which was obtained from a similar interpolation over these regions. Our result $\text{BR}(B^{+}\to\pi^{+}\mu^{+}\mu^{-})=20.4(2.1)\times 10^{-9}$ agrees with LHCb, and is more precise than the best previous theoretical estimate Ali et al. (2014) because we use $f_{T}$ directly, which avoids a large uncertainty from varying the matching scale $\mu$.

Outlook – The largest uncertainty in our determination of the $B\to\pi$ form factors is the combined error from statistics with chiral-extrapolation and discretization effects included. We will be able to reduce these with calculations on the MILC Collaboration’s recently generated four-flavor ensembles with physical light-quark masses Bazavov et al. (2013). LHCb’s measurement of $\text{BR}(B^{+}\to\pi^{+}\mu^{+}\mu^{-})$ will improve, and Belle II expects to observe the neutral decay mode $B^{0}\to\pi^{0}\ell^{+}\ell^{-}$. If a deviation from the Standard Model is observed, our form factors can be used to compute other observables such as asymmetries, thereby providing information about new heavy particles, such as their masses, spin, and couplings.

Note added: after this paper was submitted for publication, the LHCb experiment announced a new measurement for the $B^{+}\to\pi^{+}\mu^{+}\mu^{-}$ differential decay rate Tekampe (2015). The new results are shown in Fig. 3. The large difference in the lowest $q^{2}$ bin is due to the presence of light ($\rho,\omega,\phi$) resonances, whose effects are important but cannot be estimated in a model-independent manner. Given the present experimental and theoretical uncertainties, it is too early to discern possible new-physics contributions to this process.

Acknowledgements – We thank Ulrik Egede and Tobias Tekampe from LHCb for useful correspondence. Computations for this work were carried out with resources provided by the USQCD Collaboration, the Argonne Leadership Computing Facility, the National Energy Research Scientific Computing Center, and the Los Alamos National Laboratory, which are funded by the Office of Science of the United States Department of Energy; and with resources provided by the National Institute for Computational Science, the Pittsburgh Supercomputer Center, the San Diego Supercomputer Center, and the Texas Advanced Computing Center, which are funded through the National Science Foundation’s Teragrid/XSEDE Program. This work was supported in part by the U.S. Department of Energy under Grants No. DE-FG02-91ER40628 (C.B.), No. DE-FC02-06ER41446 (C.D., L.L., S.-W.Q.), No. DE-SC0010120 (S.G.), No. DE-FG02-91ER40661 (S.G.), No. DE-FG02-13ER42001 (D.D., A.X.K.), No. DE-FG02-91ER40664 (Y.M.), and No. DE-FG02-13ER41976 (D.T.); by the National Science Foundation under Grants No. PHY-1067881, No. PHY-10034278 (C.D., L.L., S.-W.Q..), No. PHY-1417805 (J.L., D.D.), and No. PHY-1316748 (R.S.); by the URA Visiting Scholars’ program (A.X.K., Y.M.); by the MINECO (Spain) under Grant FPA2013-47836-C3-1-P, and the Ramón y Cajal program (E.G.); by the Junta de Andalucía (Spain) under Grants FQM-101 and FQM-6552 (E.G.); by the European Commission under Grant No. PCIG10-GA-2011-303781 (E.G.); by the German Excellence Initiative, the European Union Seventh Framework Programme under grant agreement No. 291763, and the European Union’s Marie Curie COFUND program (A.S.K); and by the Basic Science Research Program of the National Research Foundation of Korea (NRF) funded by the Ministry of Education (No. 2014027937) and the Creative Research Initiatives Program (No. 2014001852) of the NRF grant funded by the Korean government (MEST) (J.A.B). Brookhaven National Laboratory is supported by the U. S. Department of Energy under Contract No. DE-SC0012704. Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy.