Branching Fraction of the Decay $B^+ \to \pi^+ \tau^+ \tau^-$ and Lepton Flavor Universality Test via the Ratio $R_\pi (\tau/\mu)$

Among (semi)leptonic rare $B$-decays induced by the $b \to d$ flavor changing neutral current, the decay $B^+ \to \pi^+ \mu^+ \mu^-$ is the only one observed so far experimentally. Related decays involving the $e^+e^-$ and $\tau^+ \tau^-$ pairs are the targets for the ongoing experiments at the LHC, in particular LHCb, and Belle II. The muonic and electronic semileptonic decays have almost identical branching fractions in the Standard Model (SM). However, the tauonic decay $B^+ \to \pi^+ \tau^+ \tau^-$ differs from the other two due to the higher reaction threshold which lies slightly below the $\psi (2S)$-resonance. We present calculations of the ditauon ($\tau^+ \tau^-$) invariant-mass distribution and the branching fraction ${\rm Br} (B^+ \to \pi^+ \tau^+ \tau^-)$ in the SM based on the Effective Electroweak Hamiltonian approach, taking into account also the so-called long-distance contributions. The largest theoretical uncertainty in the short-distance part of the decay rates is due to the $B \to \pi$ form factors, which we quantify using three popular parametrizations. The long-distance contribution can be minimized by a cut on the ditauon mass $m_{\tau^+ \tau^-}>M_{\psi (2S)}$. Once available, the branching fractions in the tauonic and muonic (and electronic) modes provide stringent test of the lepton flavor universality in the $b \to d$ transitions. We illustrate this by calculating the ratio $R_\pi (\tau/\mu) \equiv {\rm Br} (B^+ \to \pi^+ \tau^+ \tau^-)/{\rm Br} (B^+ \to \pi^+ \mu^+ \mu^-)$ in the SM for the total and binned ratios of the branching fractions.


Introduction
Rare bottom-hadron decays induced by the quarklevel Flavor Changing Neutral Current (FCNC) transitions b → s and b → d are of special interest as they allow us to test the SM precisely and search for possible deviations from the Standard Model (SM). FCNC processes in the SM are governed by the GIM mechanism [1], which allows such transitions only through higher-order electroweak (loop) diagrams. In particular, semileptonic rare decays are a very useful tool for testing the Lepton Flavor Universality (LFU), a linchpin of the electroweak sector of the SM. Semileptonic decays due to the b → s currents such as B ± → K ( * )± µ + µ − , B 0 → K ( * )0 µ + µ − , and B 0 s → ϕµ + µ − and their electronic counterparts, while suppressed by the loops, are favored by the quark mixing Cabibbo-Kobayashi-Maskawa (CKM) matrix [2,3]. Hence, there is plenty of data available on their branching fractions and decay characteristics, such as the leptonpair invariant-mass and angular distributions [4][5][6][7][8][9][10]. Some of these measurements were found to be not in accord with the SM-based predictions, triggering searches for better models incorporating physics beyond the Standard Model (BSM) [11][12][13][14][15][16].
An important issue in these decays is the interference between the short (perturbative)-and long (non-perturbative)-distance contributions. The standard experimental procedure is to exclude the dilepton invariant-mass squared (q 2 ) spectrum close to the J/ψand ψ(2S )-resonances, and extract the short-distance part of the spectrum from the rest. Measurements of the phase difference between the short-and longdistance amplitudes in the B + → K + µ + µ − decay has been undertaken by the LHCb collaboration based on data collected in 2011 and 2012 [17]. Their analysis shows that the phases of the J/ψ-and ψ(2S )-mesons are important near the resonance masses, due to their small decay widths, but their influence on the dilepton invariant mass spectrum in other regions is small. In addition, the branching fractions of the higher charmonium states: ψ(3770), ψ(4040), ψ(4160), and ψ(4415), were measured. This analysis is potentially helpful for studies of other semileptonic decays, B + → π + ℓ + ℓ − , in particular. For the B → πe + e − and B → πµ + µ − , also light vector mesons, ρ 0 , ω and ϕ, give sizable contributions to the lepton invariant-mass distribution around q 2 ∼ 1 GeV 2 .
Dedicated searches of possible LFU-violations in rare decays due to the b → s currents have been undertaken by the LHCb collaboration [18][19][20] in terms of the ratios R K ( * ) (µ/e) ≡ Br(B → K ( * ) µ + µ − )/Br(B → K ( * ) e + e − ), measured in selected bins in the dilepton invariant mass squared. This data hinted at LFUviolation, typically reaching three standard deviation from the SM. Similar analysis by the Belle collaboration [21,22], on the other hand, yielding R K ( * ) (µ/e) = 1.03 +0. 28 −0.24 ± 0.01 for q 2 ∈ (1.0, 6.0) GeV 2 , while consistent with the SM is less conclusive due to larger experimental errors. However, recent measurements of the ratios R K ( * ) (µ/e) in the low-and central-q 2 parts of the spectrum by the LHCb collaboration [23,24] are found in almost perfect agreement with the SM predictions [25,26]. This data, based on 9 fb −1 integrated luminosity, with improved understanding of the background and tighter electron particle identification, supersedes the earlier LHCb data.
Data on the FCNC semileptonic b → d transitions is rather sparse. For decays induced by the b → dℓ + ℓ − transition, where ℓ = e, µ, τ, the B + → π + µ + µ − decay is so far the only mode observed in the B-meson sector, first reported by the LHCb col-laboration in 2012 [41] and analyzed in detail in 2015 [42]. The measured dimuon invariant mass distribution in the B + → π + µ + µ − decay is in good agreement with theoretical predictions in the SM [43][44][45] in almost all regions of the spectrum, except the lowest q 2 -part. In this region, experimental data significantly exceeds theoretical predictions based on the short-distance contribution [42]. Taking into account the sub-leading weak annihilation (WA) and long-distance (LD) contributions, however, gives better agreement between theoretical predictions and experimental data [46][47][48]. We also note the evidences for the B 0 → π + π − µ + µ − decay with a significance of 4.8σ [49], the B 0 s → K * 0 µ + µ − decay at 3.4σ [50], reported by the LHCb collaboration, and the observation of the Λ 0 b → pπ − µ + µ − decay in the bottombaryon sector by the same collaboration [51], all of them are mediated by the b → dℓ + ℓ − transition. A model-independent analysis of the |∆b| = |∆d| = 1 processes to test the SM and probe flavor patterns of new physics was undertaken in [52]. Constraints on Wilson coefficients are obtained from global fits from data on exclusive B + → π + µ + µ − , B s →K * 0 µ + µ − , B 0 → µ + µ − , and inclusive radiative B → X d γ decays. Being consistent with the SM, these fits leave a sizable room for new physics. The ratio R π (µ/e), involving the branching ratios of B ± → π ± ℓ + ℓ − for ℓ ± = e ± , µ ± , has been studied theoretically at great length to search for the LFU-violation [53] in the semileptonic b → d sector. At present, however, there is no data available on the ratio R π (µ/e).
Precision tests of LFU involving the decays b → (s, d) τ + τ − remain to be undertaken. Compared to the b → (s, d) e + e − and b → (s, d) µ + µ − modes, they have a reduced phase space, in addition to the experimental difficulty of reconstructing the τ ± -leptons. These modes will be targeted at the LHCb and Belle II experiments. Theoretically, they have the advantage of being relatively free of the LD-contributions and the form factors involved in the SD-piece can eventually be calculated precisely on the lattice. In the b → s sector, the decays B → K ( * ) τ + τ − have recently been studied theoretically from the BSM physics point of view [54][55][56][57]. Currently only weak experimental upper limits on these decays are available, with Br(B + → K + τ + τ − ) < 2.25 × 10 −3 by the BaBar collaboration [58] and Br(B 0 → K * 0 τ + τ − ) < 2.0 × 10 −3 by Belle [59], both obtained at 90% CL.

Effective Hamiltonian for the
Our analysis is carried out in the Effective Electroweak Hamiltonians approach [63,64], where the SM heavy degrees of freedom (W ± , Z 0 , t) are absent. This effective theory also does not contain photons and gluons with energies exceeding the mass of the b-quark, m b , which represents the largest energy scale of the theory. Photons and gluons with lower energies are included using the QED and QCD Lagrangians. Rare semileptonic decays of the B-mesons involving the b → s and b → d FCNC transitions are calculated in this framework, of which the b → d part has the form: where G F is the Fermi constant, V q 1 q 2 are the CKM matrix elements satisfying the unitary condition V ud V * ub + V cd V * cb + V td V * tb = 0 which can be used to eliminate one of their products, and C j (µ) are Wilson coefficients determined at the scale µ. For the operators P j (µ), the following basis is chosen [64,65]: where p = u, c is the quark flavor, are the generators of the color S U(3) C -group, L, R = (1 ∓ γ 5 ) /2 are the left-and right-handed fermionic projectors, F µν and G A µν are the electromagnetic and gluon field strength tensors, respectively, m b and m d are the b-and d-quark masses of which the d-quark mass is neglected, σ µν = i γ µ γ ν − γ ν γ µ /2, and α em = e 2 /(4π) is the fine structure constant. The summation over q and ℓ denotes sums over all quarks (except the t-quark) and charged leptons, respectively. The Wilson coefficients C j (µ), which depend on the renormalization scale µ, are calculated at the matching scale µ W ∼ m W , where m W is the W-boson mass, as a perturbative expansion in the strong coupling constant α s (µ W ) [65]: which are evolved to a lower scale µ b ∼ m b using the anomalous dimensions of the above operators. They have been calculated to the next-next-leadinglog (NNLL) accuracy [65]: Numerical values of the Wilson coefficient, calculated to NLL accuracy, are presented in Table 1, where one can see that the Wilson coefficients of the QCD penguin operators, C j (m b ) with j = 3, 4, 5, 6, have much smaller values than the others.
Feynman diagrams for the B + → π + ℓ + ℓ − decay are shown in Fig. 1, where the left one denotes the P 7γ contribution, and the right one denotes the P 9ℓ and P 10ℓ contributions. The matrix elements for the B → P transition, where P is a pseudo-scalar meson, are expressed in terms of three transition form factors [66]: vector f + (q 2 ), scalar f 0 (q 2 ), and tensor transferred to the lepton pair: where m B and m P are the B-and pseudo-scalar meson masses, respectively.
Taking into account the sub-leading contributions, the differential branching fraction is as follows [47]: where S P is the isospin factor of the final meson (S π ± = 1 and S π 0 = 1/2 for the π-mesons, the case of our interest in this paper), C eff 7,9,10 are the effective Wilson coefficients including the NLO QCD corrections [67], L BP A (q 2 ) is the Weak-Annihilation (WA) contribution, ∆C BP V (q 2 ) is the Long-Distance (LD) contribution, and is the kinematical function encountered in three-body decays (the triangle function). Note that the differential branching fraction for the decay with the τ + τ −pair production differs from its counterparts with e + e − and µ + µ − due to the important role of the scalar form factor, f 0 (q 2 ). In the electronic and muonic modes, its contribution is chirally suppressed by m 2 e and m 2 µ , respectively, while this no longer holds for the τ + τ − case.
The WA contribution is calculated in the so-called Large Energy Effective Theory (LEET) [68] and has a significant effect for q 2 ≲ 1 GeV 2 only, so its inclusion makes sense for the B ± → π ± e + e − and B ± → π ± µ + µ − decays, but it is irrelevant for the B ± → π ± τ + τ − case having the q 2 -threshold above 12 GeV 2 .
Two-particle decays, B → Vπ, where V is a neutral vector meson, followed by the leptonic decay V → ℓ + ℓ − determine the LD-contributions. They can be represented as follows [46]: where m V , f V and Γ tot V are the mass, decay constant and total decay width of the vector meson, respectively, k V is a valence quark content factor, A p BVπ (p = u, c) are the transition amplitudes, and the free parameter q 2 0 = −1.0 GeV 2 is chosen to achieve a better convergence in the denominator of (19). The differential branching fraction (16) involves three B → P form factors. They are scalar functions of q 2 , discussed for the B → π case in the next section.

Form Factor Parametrizations
Among the available parametrizations of the B → π transition form factors (FF) known in the literature, we chose those which are based on analyticity, crossing symmetry and the QCD dispersion relations. They are represented as a series in powers of the function z(q 2 , q 2 0 ) projecting q 2 into the unit ellipse in the complex plane 1 .
The second one is the Bourrely-Caprini-Lellouch (BCL) parametrization [61] (i = +, T ): Here, z(q 2 , q 2 0 ) is the same as in Eq. (21) and the form factors are calculated by truncating the series at N = 4.
The third type is the modified Bourrely-Caprini-Lellouch (mBCL) parametrization [62] (i = +, T ): wherez k (q 2 , q 2 0 ) = z k (q 2 , q 2 0 ) − z k (0, q 2 0 ). The function z(q 2 , q 2 0 ) is defined in Eq. (21) and q 2 0 takes the optimal value (24). Here, unlike other types of the f 0 (q 2 ) parametrizations, this form factor has a pole but at higher q 2 -at the scalar B 0 -meson mass squared, m 2 B 0 . This state is not yet observed experimentally and its mass is taken from theory. We set m B 0 = 5.54 GeV, as was used in the determination of the expansion coefficients b (0) k [62]. Note that the Dispersion Matrix (DM) method was suggested in [69] to describe the FFs by using also analyticity, crossing symmetry and the QCD dispersion relations. This method is based on the nonperturbative determination of the dispersive bounds and describes in a model-independent way the FFs in the full kinematical range, starting from existing Lattice QCD data at large momentum transfer, without a series expansion in powers of z(q 2 , q 2 0 ). It was already applied to the semileptonic B → πℓν ℓ decays [70] and can be also used for the analysis of semileptonic FCNC B-meson decays. 4. Numerical Analysis of the B + → π + τ + τ − Decay

Perturbative Contribution
The distribution in the tau-pair invariant mass calculated in perturbation theory for three types of form factor parametrizations is presented in Fig. 2. The spread shown in these distributions reflect the convoluted uncertainties in the scale parameter µ, entering via the Wilson coefficients by varying it in the range m b /2 ≤ µ ≤ 2m b , and the input value of the CKM matrix element V td = (8.54 ± 0.30) × 10 −3 [36]. Numerical results for the total branching fraction for the three FF parametrizations used in this work are consistent with each other within uncertainties as shown in Table 2. In working out the numerical values, the expansion coefficients of the BGL parametrization are taken from [43], where the data on the CCprocess B → πℓν ℓ decay are fitted, and the relations between the B → K and B → π form factors are used. The values of the BCL parametrization coefficients were obtained within the framework of Lattice QCD (LQCD) [71,72]. The values of the mBCL parametrization coefficients are obtained by the combined use of the Light-Cone Sum Rules (LCSR) and LQCD [62]. All the expansion coefficients are collected in Appendix A. The entries in Table 2 are also consistent with the existing theoretical predictions in the literature, for example, Br FG th (B + → π + τ + τ − ) = (7.0 ± 0.7) × 10 −9 [44] and Br WX th (B + → π + τ + τ − ) = (6.0 +2.6 −2.1 ) × 10 −9 [73]. It is customary to compare data and theoretical distributions in bins of q 2 : To that end, we plot the theoretical results for the partial branching ratio (∆Br) τ π (q 2 1 , q 2 2 ) in bins of the ditauon invariant-mass squared using the three FF parametrization in Fig. 3, and collect the corresponding values of the partial branching fractions, integrated over the indicated ranges, in Table 3. For comparison, the Lattice results [72] are also shown in the last column. The errors shown are from the CKM matrix element, form factors, variation of the high and low matching scales, and the quadrature sum of all other contributions, respectively. We note that the BGL parametrization predictions are in good agreement with the Lattice-based estimates, both of which are, however, systematically higher than the predictions based on the BCL and mBCL ones in each bin. The same also holds for the total branching fraction (see Table 2).
There are yet more vector charmonium resonances with masses above the ψ(3S )-meson mass, ψ(4040), ψ(4160), ψ(4230), ψ(4360), and ψ(4415), which also have purely leptonic decay modes. However, as they decay strongly into the DD-pair etc., their electronic or muonic decay rates are also of order of 10 −5 [36], similar to the case of the ψ(3S )-meson, as shown in Table 4. It follows that their contributions to the B + → π + τ + τ − branching fraction are of the same order of magnitude as from the ψ(3S )-meson. Consequently, we drop the contribution from all the strongly decaying charmonium resonances (ψ(3S ) and higher), and consider the LD-contribution from the narrow ψ(2S )meson only.
Since the LD contributions (19) depend on the choice of the amplitude phases δ (u) ψ(2S ) and δ (c) ψ(2S ) , we present the total branching fraction of the B + → π + τ + τ − decay, including the ψ(2S ) LD-contribution, and its dependence on the assumed values of the strong phases in Tables 5, 6 and 7 for the BGL, BCL and mBCL parametrizations of the form factors, respectively. As can be seen, the variation of the branching fraction on the strong phases is not very marked, and is similar to the errors shown from the SD contribution. The central value including the LD contribution is given for δ (u) ψ(2S ) = 0 and δ (c) ψ(2S ) = 3π/4. The ditauon invariant mass distribution including the LD contribution from the ψ(2S )-meson is presented in Fig. 4 for the BGL, BCL and mBCL form factors. The vertical solid line at q 2 = 15 GeV 2 in the plots indicates the kinematical cut to exclude the dominant ψ(2S ) contribution.
To study the dependence of theoretical results on the choice of the FF parametrization, we plot the partial ratio R (τ/µ) π (q 2 min , q 2 max ) in bins of the dilepton invariant mass squared, shown in Fig. 5. Numerical values of Table 7: The total branching fraction for the B + → π + τ + τ − decay in the mBCL parametrization including the LD contribution from the ψ(2S )-meson for the various assumed values of the amplitude phases. SDC means the short-distance (perturbative) contribution. this ratio are shown in Table 8, obtained by integrating the partial branching ratio over the indicated q 2ranges. The errors shown take into account the uncertainties due to the factorization-scale, the CKM matrix element V td , and the form factor errors. Theoretical predictions for the total ratio for the BGL, BCL and mBCL parametrizations are as follows: They agree with each other within the indicated uncertainties. The central values for R π (τ/µ) lie in the range 0.30 − 0.45. The main uncertainty on R π (τ/µ), as opposed to the very precise ratio R K ( * ) (µ/e), is due to the form factors. These results are potentially useful in testing the lepton flavor universality in the b → dℓ + ℓ − sector.

Summary and outlook
We have presented theoretical predictions for the branching ratio Br(B + → π + τ + τ − ) and the ditauon invariant-mass distribution at NLO accuracy in the SM, using three popular parametrizations of the B → π form factors, known in the literature as the BGL [60], BCL [61], and mBCL [62]. In the SM, LFU holds, which relates the decay B + → π + τ + τ − to the observed decay B + → π + µ + µ − . In extensions of the SM, the LFU hypothesis can be easily violated, of which the leptoquark models are the primary candidates [75,76]. In view of this, we have worked out the ratio of the branching ratios R π (τ/µ). Together with the corresponding ratios R π (µ/e) and R π (τ/e), their measurement will provide a precision test of the lepton flavor universality in the FCNC b → d transitions. Suffice to say that none of these ratios have been subjected to experimental scrutiny so far. Of these, the ratio R π (µ/e) has been theoretically investigated in [53]. We have concentrated here on the ratios R π (τ/µ) and R π (τ/e).
The decay B + → π + τ + τ − involves all three form factors, f + (q 2 ), f 0 (q 2 ), and f T (q 2 ). The uncertainties in Br(B + → π + τ + τ − ) arise from the FF parametrizations, scale-dependence of the Wilson coefficients, and the CKM matrix element. Numerical values of the SD (perturbative) contribution to Br(B + → π + τ + τ − ) are listed in Table 2. Partial branching fractions of the B + → π + τ + τ − decay, (∆Br) τ π (q 2 min , q 2 max ), in bins of ditauon invariant mass squared, are shown in Fig. 3 and displayed in Table 3. In addition, the LD-contribution from the process B → πV → πℓ + ℓ − is calculated. Experimental data on the masses, partial and total decay widths of the charmonium states are given in Table 4. Due to the small decay width of the ψ(2S )-meson, its LD-contribution is concentrated close to the B + → π + τ + τ − reaction threshold and can be largely eliminated by a cut on the ditauon invariant mass squared (q 2 ≥ 15 GeV 2 ) (see Fig. 4). We have given arguments why the contribution from the higher charmonium resonances are numerically small, and hence they do not change the SDcontribution in this region perceptibly. Taking into account the LD contribution from the ψ(2S )-resonance, numerical results for the branching ratio Br(B + → π + τ + τ − ) are given in Table 5 for the BGL form factors. Various entries in this table correspond to using the indicated strong phases. The corresponding results for the BCL form factors are given in Table 6 and for the mBCL form factors in Table 7. Since the BGL parametrization and the Lattice-QCD based estimates are rather close to each other, our estimate of the total branching fraction is Br th (B + → π + τ + τ − ) = 7.5×10 −9 , with an uncertainty of about 10%.
Numerical values for the ratio R π (τ/µ) are given in Eqs. (37)-(39) for three parametrizations chosen. The central values lie in the range 0.30 − 0.45. The main uncertainty on R π (τ/µ), as opposed to the very precise ratio R K ( * ) (e/µ), is due to the form factors. Partial ratios R (τ/µ) π (q 2 min , q 2 max ) in bins of ditauon invariant mass squared are shown in Fig. 5. These ratios can be more precisely calculated in the future by progress in Lattice QCD.
The branching ratio for B + → π + τ + τ − , integrated over the region q 2 ≥ 15 GeV 2 , has a factor of about 3 suppression compared with the B + → π + µ + µ − total branching fraction, which has already been measured. Of course, one has to factor in the experimental efficiency of reconstructing the τ + τ − pair, but a precise measurement is still feasible for the anticipated integrated luminosities at Belle II and LHCb experiments. Once sufficient data is collected, also the measurements of the various asymmetries, such as the isospin-asymmetry involving B 0 → π 0 ℓ + ℓ − and B ± → π ± τ + τ − , and CP-violating asymmetries involv- ing B + → π + τ + τ − and B − → π − τ + τ − become interesting, which should be looked at theoretically.