1 Introduction

Recently significant experimental progress has been achieved in studying weak heavy-to-light decays of \(B\) mesons. For the semileptonic \(B\rightarrow \pi l\nu _l\) and \(B\rightarrow \rho l\nu _l\) decays not only total decay branching fractions were measured by Belle and BaBar Collaborations [13] rather precisely but also differential distributions in rather narrow \(q^2\) bins. Such measurements are very important since they provide the test of the momentum dependence of the weak differential decay branching fractions and thus significantly constrain the theoretical models. It also allows us to extract the Cabibbo–Kobayashi–Maskawa (CKM) matrix element \(V_{ub}\) from exclusive decay channels with better precision and confront it with the value obtained from inclusive semileptonic decays. Moreover, recently the LHCb Collaboration [4] reported first observation of the rare \(B^+\rightarrow \pi ^+\mu ^+\mu ^-\) decay with the branching fraction \(Br(B^+\rightarrow \pi ^+\mu ^+\mu ^-)=(2.3\pm 0.6\pm 0.1)\times 10^{-8}\). Such decays are governed by the flavor changing neutral current and thus are very sensitive to the contributions of new intermediate particles and interactions. Therefore the study of rare \(B\) decays is important for constraining the theories which go beyond the standard model. Since these decays are induced by loop diagrams they are suppressed by at least three orders of magnitude compared to corresponding heavy-to-light semileptonic \(B\) decays. The observation of the rare \(B\rightarrow \pi \mu ^+\mu ^-\) decay signifies an important progress since this decay is stronger CKM suppressed compared with better studied \(B\rightarrow K^{(*)}\mu ^+\mu ^-\) decays.

Theoretical investigation of weak decays requires the determination of the decay matrix elements of the weak current between meson states. It is convenient to parametrize these decay matrix elements in terms of the invariant form factors. The calculation of these form factors demands application of the nonperturbative methods. Since these decays are governed by the heavy-to-light quark transitions they have a very broad kinematical range. Various theoretical approaches were used to calculate these form factors. However most of the employed methods allow determination of the momentum transfer dependence of the form factors only in a rather limited \(q^2\) range. For example light-cone QCD sum rules are applicable in the large recoil region (\(q^2\approx 0\)) while lattice QCD provides results at small recoil (\(q^2\approx q^2_\mathrm{max}\)). Therefore some model assumptions and/or parametrizations should be used to extrapolate the results in the whole range of \(q^2\), thus introducing additional theoretical uncertainties.

In our previous papers [5, 6] we investigated the rare \(B\rightarrow K^{(*)}l^+l^-\) and rare \(B_s\) decays in the framework of the relativistic quark model with the QCD-motivated quasipotential of the quark–antiquark interaction. Calculating the decay form factors we systematically took into account relativistic effects including transformations of the meson wave functions from the rest to the moving reference frame and contributions of the intermediate negative-energy states. All form factors are expressed through the usual overlap integrals of the meson wave functions, which are known from the previous mass spectrum considerations [711]. The important advantage of the developed method consists in the fact that it provides explicit calculation of the momentum transfer dependence of the form factors in the whole kinematical range, thus improving the reliability of the obtained results. Here we apply this approach for studying the rare \(B\rightarrow \pi (\rho )l\bar{l}\) decays. Evaluation of such decay branching fractions requires calculation of the matrix elements of weak vector, axial vector and tensor currents. In Ref. [12] we calculated the form factors parametrizing matrix elements of the weak vector and axial vector currents for the \(B\rightarrow \pi \) and \(B\rightarrow \rho \) transitions and on this basis studied the corresponding semileptonic decays \(B\rightarrow \pi (\rho )l\nu _l\). We first confront the predictions of our model with new detailed Belle and BaBar data [13] on heavy-to-light semileptonic decays. Such comparison provides additional test of the \(q^2\) dependence of the form factors. Then we apply the model for the calculation of the tensor form factors. On this basis the differential distributions and total decay branching fractions as well as the forward–backward asymmetry and the \(\rho \) polarization fractions are calculated and confronted with available experimental data and other theoretical predictions.

2 Form factors of the weak \(B\) meson transitions to \(\pi \) and \(\rho \) mesons in the relativistic quark model

The matrix elements of the weak current for the heavy-to-light \(b\rightarrow q\) (\(q=u,d\)) weak transitions between the initial \(B\) meson and final pseudoscalar \(\pi \) or vector \(\rho \) mesons are usually parametrized by the following set of invariant form factors.

  1. (a)

    \(B \rightarrow \pi \) weak decays

    $$\begin{aligned}&\langle \pi (p_{\pi })|\bar{q} \gamma ^\mu b|B(p_{B})\rangle \nonumber \\&\quad =f_+(q^2)\left[ p_{B}^\mu + p_{\pi }^\mu - \frac{M_{B}^2-M_{\pi }^2}{q^2}\ q^\mu \right] \nonumber \\&\qquad +\, f_0(q^2)\frac{M_{B}^2-M_{\pi }^2}{q^2}\ q^\mu , \end{aligned}$$
    (1)
    $$\begin{aligned} \langle \pi (p_{\pi })|\bar{q} \gamma ^\mu \gamma _5 b|B(p_{B})\rangle \!&= \!\langle \pi (p_\pi )|\bar{q} \sigma ^{\mu \nu }\gamma _5 q_\nu b|B(p_{B})\rangle \nonumber \\ \!&= \! 0, \end{aligned}$$
    (2)
    $$\begin{aligned}&\langle \pi (p_\pi )|\bar{q} \sigma ^{\mu \nu }q_\nu b|B(p_{B})\rangle = \frac{if_T(q^2)}{M_{B}+M_\pi }\nonumber \\&\quad \times [q^2(p_{B}^\mu +p_\pi ^\mu )-(M_{B}^2-M_\pi ^2)q^\mu ], \end{aligned}$$
    (3)
  2. (b)

    \(B \rightarrow \rho \) weak decays

    $$\begin{aligned}&\langle {\rho }(p_{\rho })|\bar{q} \gamma ^\mu b|B(p_{B})\rangle =\frac{2iV(q^2)}{M_{B}+M_{\rho }} \epsilon ^{\mu \nu \tau \sigma }\epsilon ^*_\nu p_{B\tau } p_{{\rho }\sigma },\end{aligned}$$
    (4)
    $$\begin{aligned}&\langle {\rho }(p_{\rho })|\bar{q} \gamma ^\mu \gamma _5 b|B(p_{B})\rangle =2M_{\rho } A_0(q^2)\frac{\epsilon ^*\cdot q}{q^2}\ q^\mu \nonumber \\&\quad +\,(M_{B}+M_{\rho })A_1(q^2)\left( \epsilon ^{*\mu }-\frac{\epsilon ^*\cdot q}{q^2}\ q^\mu \right) \nonumber \\&\quad -\,A_2(q^2)\frac{\epsilon ^*\cdot q}{M_{B}+M_{\rho }}\left[ p_{B}^\mu + p_{\rho }^\mu -\frac{M_{B}^2-M_{\rho }^2}{q^2}\ q^\mu \right] ,\nonumber \\ \end{aligned}$$
    (5)
    $$\begin{aligned} \langle \rho (p_{\rho })|\bar{q} i\sigma ^{\mu \nu }q_\nu b|B(p_{B})\rangle =2T_1(q^2) \epsilon ^{\mu \nu \tau \sigma } \epsilon ^*_\nu p_{{\rho }\tau } p_{{B}\sigma },\nonumber \\ \end{aligned}$$
    (6)
    $$\begin{aligned}&\langle \rho (p_{\rho })|\bar{q} i\sigma ^{\mu \nu }\gamma _5q_\nu b|B(p_{B})\rangle \!=\! T_2(q^2)[(M_{B}^2-M_{\rho }^2)\epsilon ^{*\mu } \nonumber \\&\quad -(\epsilon ^*\cdot q)(p_{B}^\mu + p_{\rho }^\mu )]\nonumber \\&\quad +T_3(q^2)(\epsilon ^*\cdot q)\left[ q^\mu -\frac{q^2}{M_{B}^2-M_{\rho }^2} (p_{B}^\mu +p_{\rho }^\mu )\right] , \end{aligned}$$
    (7)

where \(q=p_{B}- p_{\pi (\rho )}\) is the four-momentum transfer, \(M_{B,\pi (\rho )}\) are the initial and final meson masses, and \(\epsilon _\mu \) is the polarization vector of the final \(\rho \) meson.

At the maximum recoil point (\(q^2=0\)) these form factors satisfy the following conditions:

$$\begin{aligned}&f_+(0)=f_0(0),\\&A_0(0)=\frac{M_{B}+M_\rho }{2M_{\rho }}A_1(0) -\frac{M_{B}-M_{\rho }}{2M_{\rho }}A_2(0),\\&T_1(0)=T_2(0). \end{aligned}$$

In this paper we use the relativistic quark model based on the quasipotential approach and QCD for the calculation of the form factors of weak \(B\) decays to final \(\pi \) or \(\rho \) mesons. The meson is described by the covariant single-time wave function which satisfies the three-dimensional relativistically invariant Schrödinger-like equation with the QCD-motivated interquark potential [7, 8]

$$\begin{aligned} \left( \frac{b^2(M)}{2\mu _{R}}-\frac{{\mathbf p}^2}{2\mu _{R}}\right) \Psi _{M}({\mathbf p}) =\int \frac{d^3 q}{(2\pi )^3} V({\mathbf p,q};M)\Psi _{M}({\mathbf q}), \end{aligned}$$
(8)

where the relativistic reduced mass is

$$\begin{aligned} \mu _{R}=\frac{E_1E_2}{E_1+E_2}=\frac{M^4-(m^2_1-m^2_2)^2}{4M^3}, \end{aligned}$$
(9)

with the on-mass-shell energies

$$\begin{aligned} E_1=\frac{M^2-m_2^2+m_1^2}{2M}, \quad E_2=\frac{M^2-m_1^2+m_2^2}{2M}, \end{aligned}$$

and \(M=E_1+E_2\) is the meson mass, \(m_{1,2}\) are the quark masses, and \({\mathbf p}\) is their relative three-momentum. In the center-of-mass system the on-mass-shell relative momentum squared \(b^2(M)\) is expressed through the meson and quark masses:

$$\begin{aligned} b^2(M) =\frac{[M^2-(m_1+m_2)^2][M^2-(m_1-m_2)^2]}{4M^2}. \end{aligned}$$
(10)

The kernel \(V({\mathbf p,q};M)\) of Eq. (8) is the quasipotential operator of the quark–antiquark interaction. It is constructed with the help of the off-mass-shell scattering amplitude, projected onto the positive-energy states. The explicit expression for the corresponding quasipotential \(V({\mathbf p,q};M)\) can be found in Refs. [7, 8].

The constituent quark masses \(m_c=1.55\) GeV, \(m_b=4.88\) GeV, \(m_u=m_d=0.33\) GeV, \(m_s=0.5\) GeV and the parameters of the linear confining potential \(A=0.18\) GeV\(^2\) and \(B=-0.3\) GeV have been fixed previously and have values typical for quark models. The value of the mixing coefficient of vector and scalar confining potentials \(\varepsilon =-1\) has been determined from the consideration of charmonium radiative decays [7, 8] and heavy quark effective theory. The universal Pauli interaction constant \(\kappa =-1\) has been fixed from the analysis of the fine splitting of heavy quarkonia \({ }^3P_J\) - states [7, 8]. In this case, the long-range chromomagnetic quark moment \((1+\kappa )\) vanishes in accordance with the flux-tube model.

The matrix element of the weak current between the \(B\) meson with mass \(M_{B}\) and momentum \(p_{B}\) and a final (\(F=\pi \) or \(\rho \)) meson with mass \(M_{F}\) and momentum \(p_{F}\) is given [7, 8] by the expression

$$\begin{aligned}&\langle F(p_{F}) \vert J^W_\mu \vert B(p_{B})\rangle \nonumber \\&\quad = \int \frac{d^3p\, d^3q}{(2\pi )^6} \bar{\Psi }_{{F}\,{\mathbf p}_F}({\mathbf p})\Gamma _\mu ({\mathbf p},{\mathbf q})\Psi _{B\,{\mathbf p}_{B}}({\mathbf q}), \end{aligned}$$
(11)

where \({\mathbf p},{\mathbf q}\) are relative quark momenta, \(\Gamma _\mu ({\mathbf p},{\mathbf q})\) is the two-particle vertex function. Here \(\Psi _{M\,{\mathbf p}_M}({\mathbf p})\) are the meson (\(M=B,{F})\) wave functions projected onto the positive-energy states of quarks (\(q_1=b,u,d\) and \(q_2=u,d\)) and boosted to the moving reference frame with momentum \({\mathbf p}_M\)

$$\begin{aligned} \Psi _{M\,{\mathbf p}_M}({\mathbf p})=D_{q_1}^{1/2}(R_{L_{{\mathbf p}_M}}^W)D_{q_2}^{1/2}(R_{L_{{ \mathbf p}_M}}^W)\Psi _{M\,{\mathbf 0}}({\mathbf p}), \end{aligned}$$
(12)

where \(\Psi _{M\,{\mathbf 0}}({\mathbf p})\equiv \Psi _M({\mathbf p})\) is the wave function at rest, \(R^W\) is the Wigner rotation, \(L_{\varvec{\Delta }}\) is the Lorentz boost from the meson rest frame to a moving one, and \(D^{1/2}(R)\) is the spin rotation matrix.

Calculating the weak decay matrix elements we take into account both the leading (spectator) term \(\Gamma ^{(1)}({\mathbf p},{\mathbf q})\) of the vertex function and subleading term \(\Gamma ^{(2)}({\mathbf p},{\mathbf q})\) which takes into account contributions of the intermediate negative-energy states. The diagrams and explicit expressions for these terms can be found in Refs. [5, 6]. The previously developed methods [12] allow us to express the relativistic decay matrix elements through the usual overlap integrals of the initial and final meson wave functions in their rest frames. These wave functions are known from the meson-mass spectrum calculations [911]. Note that all considerations were done completely relativistically without employing \(v/c\) expansion. It is important to point out that the obtained expressions for the decay matrix elements are valid in the whole kinematical \(q^2\) range accessible in weak decays. This fact allows one to explicitly determine the \(q^2\) dependence of meson form factors without additional model assumptions and extrapolations, thus increasing reliability of results. The analytic expressions for the form factors are given in Refs. [5, 13]. It is important to emphasize that these form factors in the heavy quark and large recoil limits satisfy all model-independent relations imposed by heavy quark and large energy effective theories [14, 15].

The numerical values of the form factors \(f_+(q^2)\), \(f_0(q^2)\), \(V(q^2)\), \(A_i(q^2)\) (\(i=0,1,2\)) parametrizing matrix elements of vector and axial vector weak currents for the \(B\rightarrow \pi \) and \(B\rightarrow \rho \) transitions were previously calculated in Ref. [12]. Here we further extend our calculation to the form factors \(f_T(q^2)\) and \(T_i(q^2)\) (\(i=1,2,3\)) parametrizing matrix elements of the tensor and pseudotensor currents responsible for the rare \(B\rightarrow \pi (\rho )l\bar{l}\) decays. These form factors are plotted in Figs. 1 and 2.

Fig. 1
figure 1

Form factors of the weak \(B\rightarrow \pi \) transition

Fig. 2
figure 2

Form factors of the weak \(B\rightarrow \rho \) transition

To simplify the comparison of the obtained form factors with experiment and other theoretical calculations it is useful to have approximate analytic expressions for them. Our analysis shows that the weak \(B\rightarrow \pi (\rho )\) transition form factors can be well fitted by the following formulas [12, 16]:

  1. (a)

    \(F(q^2)= \{f_+(q^2),f_T(q^2),V(q^2),A_0(q^2),T_1(q^2)\}\)

    $$\begin{aligned} F(q^2)=\frac{F(0)}{\displaystyle \left( 1-\frac{q^2}{M^2}\right) \left( 1-\sigma _1 \frac{q^2}{M_{B^*}^2}+ \sigma _2\frac{q^4}{M_{B^*}^4}\right) }, \end{aligned}$$
    (13)
  2. (b)

    \(F(q^2)=\{f_0(q^2), A_1(q^2),A_2(q^2),T_2(q^2),T_3(q^2)\}\)

    $$\begin{aligned} F(q^2)=\frac{F(0)}{\displaystyle \left( 1-\sigma _1 \frac{q^2}{M_{B^*}^2}+ \sigma _2\frac{q^4}{M_{B^*}^4}\right) }, \end{aligned}$$
    (14)

where \(M=M_{B^*}\) for the form factors \(f_+(q^2),f_T(q^2), V(q^2),T_1(q^2)\) and \(M=M_{B}\) for the form factor \(A_0(q^2)\). The obtained values of \(F(0)\) and \(\sigma _{1,2}\) are given in Table 1. The accuracy of such approximation is rather high, the deviation from the calculated form factors does not exceed 1 %. The rough estimate of the total uncertainty of the form factors within our model is of order of 5 %.

Table 1 The form factors of weak \(B\rightarrow \pi (\rho )\) weak transitions

In Table 2 we compare the predictions of our model for the weak \(B\rightarrow \pi (\rho )\) decay form factors at the maximum recoil point \(q^2=0\) with other theoretical calculations [1622]. The different versions of light-cone sum rules are employed in Refs. [17, 18]. Calculations of Ref. [16] are based on the constituent quark model within relativistic dispersion approach while in Ref. [19] the covariant constituent quark model with the infrared confinement is applied. The perturbative QCD factorization approach with the inclusion of the leading and next-to-leading-order corrections is used in Refs. [20, 21]. The authors of Ref. [22] extract the form factors combining available experimental data on semileptonic \(B\rightarrow \pi l\nu _l\) decays and lattice QCD calculations of the corresponding form factors for the rare \(B\rightarrow K\) transitions within the \(SU(3)\)-breaking Ansatz. Comparison of the results presented in this table shows that, although there are some differences between the central values of predictions, in general there is a reasonable agreement between the values of these form factors at zero recoil calculated using significantly different theoretical methods.

Table 2 Comparison of theoretical predictions for the form factors of weak \(B\rightarrow \pi (\rho )\) transitions at the maximum recoil point \(q^2=0\)

Note that most of the discussed theoretical approaches allow one to calculate the form factors at a single point only or in some limited range of the recoil momentum, then some model extrapolation to the whole kinematical range should be used. The important advantage of our approach consists in the fact that it determines various decay form factors through the overlap integrals of hadron wave functions in the whole kinematically accessible range without additional assumptions and extrapolations. These wave functions are obtained by numerical solving Eq. (1) with the nonperturbative treatment of relativistic effects.

We can further test our model confronting the form factor \(f _0\) at zero recoil point (\(q^2=q^2_\mathrm{max}\)) with the Callan–Treiman-type normalization condition

$$\begin{aligned} f_0(q^2_\mathrm{max})\approx \frac{f_B}{f_\pi }, \end{aligned}$$

derived using the soft pion limit \(p\rightarrow 0\) and \(M_\pi ^2\rightarrow 0\) [23, 24]. Taking our prediction for the decay constant \(f_B=189\) MeV [25], which is well consistent with the averaged theoretical value given in Ref. [26], and the experimental value of \(f_\pi \) [26] we get \(f_0(q^2_\mathrm{max})\approx 1.45\) in good agreement with the value 1.32 given in Table 1 (see also Ref. [12]).

3 Semileptonic \(B\rightarrow \pi l\nu _l\) and \(B\rightarrow \rho l\nu _l\) decays

We start from the consideration of the semileptonic \(B\rightarrow \pi l\nu _l\) and \(B\rightarrow \rho l\nu _l\) decays. They were investigated in detail in Ref. [12], where all necessary formulas and values of branching fractions can be found. Using the new data from Belle [1, 2] and BaBar [3] on the exclusive charmless semileptonic \(B\) decays we can update our analysis in Ref. [12] as follows.

The branching fractions of such decays predicted by our model are given [12] by

$$\begin{aligned} Br(\bar{B}^0\rightarrow \pi ^+ l^-\nu )&= 8.36|V_{ub}|^2,\nonumber \\ Br(B^-\rightarrow \pi ^0 l^-\nu )&= 4.51|V_{ub}|^2,\nonumber \\ Br(\bar{B}^0\rightarrow \rho ^+ l^-\nu )&= 20.12|V_{ub}|^2,\nonumber \\ Br(B^-\rightarrow \rho ^0 l^-\nu )&= 10.87|V_{ub}|^2. \end{aligned}$$
(15)

Comparing these predictions with recent experimental data [13]

$$\begin{aligned} Br(\bar{B}^0\rightarrow \pi ^+ l^-\nu )&= (1.49\pm 0.04\pm 0.07)\times 10^{-4}\ [1],\nonumber \\ Br(\bar{B}^0\rightarrow \pi ^+ l^-\nu )&= (1.49\pm 0.09\pm 0.07)\times 10^{-4}\ [2],\nonumber \\ Br(\bar{B}^0\rightarrow \pi ^+ l^-\nu )&= (1.47\pm 0.05\pm 0.06)\times 10^{-4}\ [3],\nonumber \\ Br(B^-\rightarrow \pi ^0 l^-\nu )&= (0.80\pm 0.08\pm 0.04)\times 10^{-4}\ [2],\nonumber \\ Br(B^-\rightarrow \pi ^0 l^-\nu )&= (0.77\pm 0.04\pm 0.03)\times 10^{-4}\ [3],\nonumber \\ Br(B^0\rightarrow \rho ^+ l^-\nu )&= (3.22\pm 0.27\pm 0.24)\times 10^{-4}\ [2],\nonumber \\ Br(B^-\rightarrow \rho ^0 l^-\nu )&= (1.83\pm 0.10\pm 0.10)\times 10^{-4}\ [2],\nonumber \\ \end{aligned}$$
(16)

we find the values of the CKM matrix element \(|V_{ub}|\) presented in Table 3.

Table 3 Values of the CKM matrix element \(|V_{ub}|\times 10^3\) extracted in our model from the recent experimental data. Only experimental errors are given

Averaging these values we get the following exclusive value for the CKM matrix element \(|V_{ub}|\):

$$\begin{aligned} |V_{ub}|=(4.15\,{\pm }\,0.09_\mathrm{exp}\,{\pm }\,0.21_\mathrm{theor})\times 10^{-3} \qquad (\mathrm{exclusive}) , \end{aligned}$$
(17)

where the last error is the rough (conservative upper) estimate of the theoretical uncertainties within our model. Note that this value is consistent with the one extracted from the inclusive charmless semileptonic \(B\) decays [26]

$$\begin{aligned} |V_{ub}|=(4.41\pm 0.15^{+0.15}_{-0.17})\times 10^{-3} \qquad (\mathrm{inclusive}) . \end{aligned}$$
(18)

It is important to point out that recent data provide us not only the total branching fractions but also the partial branching fractions \(\Delta Br/\Delta q^2\) averaged over rather small \(q^2\) bins. This allows one to test rather precisely the \(q^2\) dependence of decay form factors. The comparison of our results for the partial decay rates with recent data is given in Figs. 3, 4, and 5. In Fig. 3 we confront our predictions for the semileptonic decay of the neutral \(\bar{B}^0\) meson to the charged \(\pi ^+\) meson with recent untagged and tagged data from Belle [1, 2] and data from BaBar [3], while in Fig. 4 the corresponding predictions and the Belle data for the decay of the charged \(B^-\) meson to the neutral \(\pi ^0\) meson are presented. Differential branching fractions for decays of the charged and neutral \(B\) mesons to \(\rho \) mesons are plotted in Fig. 5. Here only Belle data are available [2]. From these figures we see that the reasonable agreement of our theoretical results and data is observed both for the semileptonic \(B\) decays to the pseudoscalar \(\pi \) and vector \(\rho \) mesons. In most cases our predictions agree with data within error bars or lie just in-between individual measurements. This comparison ensures the reliability of our approach utilizing the model form factors.

Fig. 3
figure 3

Comparison of predictions of our model with recent experimental data for the \(B^0\rightarrow \pi ^+ l^-\nu \) decay [13]

Fig. 4
figure 4

Comparison of predictions of our model with recent experimental data for the \(B^-\rightarrow \pi ^0 l^-\nu \) decay (left figure—Belle data [1], right figure—BaBar data [3])

Fig. 5
figure 5

Comparison of predictions of our model with recent Belle data [1] for the \(B\rightarrow \rho l\nu \) decay

4 Rare semileptonic \(B\rightarrow \pi (\rho ) l^+l^-\) and \(B\rightarrow \pi (\rho ) \nu \bar{\nu }\) decays

Now we apply the calculated weak decay form factors to the consideration of the rare \(B\) decays to light \(\pi \) or \(\rho \) mesons. Such decays are significantly less studied experimentally. Their theoretical description is usually based on the effective Hamiltonian \({\mathcal {H}}_\mathrm{eff}\) in which heavy degrees of freedom (gauge bosons and top quark) are integrated out. The operator product expansion allows the separation of short- and long-distance effects which are assumed to factorize. The short-distance contributions are described by the Wilson coefficients \(c_i\) which are calculated within perturbation theory, while the long-distance part is attributed to the set of the standard model operators \({\mathcal {O}}_i\).

The effective Hamiltonian for the \(b\rightarrow d l^+ l^-\) transitions can be presented [27] in the following form taking into account the unitarity of the CKM matrix:

$$\begin{aligned} {\mathcal {H}}_\mathrm{eff}&= -\frac{4G_F}{\sqrt{2}}\left[ V_{td}^*V_{tb}\sum _{i=1}^{10}c_i{\mathcal {O}}_i \right. \nonumber \\&\left. + V_{ud}^*V_{ub}\sum _{i=1}^{2}c_i\left( {\mathcal {O}}_i-{\mathcal {O}}_i^{(u)}\right) \right] , \end{aligned}$$
(19)

where \(G_F\) is the Fermi constant, \(V_{tj}\) and \(V_{uj}\) are the CKM matrix elements, \(c_i\) are the Wilson coefficients and \({\mathcal {O}}_i({\mathcal {O}}_i^{(u)})\) comprise the four-quark operator basis. Then the resulting transition amplitude is given by

$$\begin{aligned} {\mathcal {M}}&= \frac{G_F\alpha }{\sqrt{2}\pi } |V_{td}^*V_{tb}|\nonumber \\&\times \Biggl \{(\bar{d}\left[ c_9^\mathrm{eff}\gamma _\mu (1\!-\!\gamma _5)\!-\! \frac{2m_b}{q^2}c_7^\mathrm{eff}i\sigma _{\mu \nu }q^\nu (1\!+\!\gamma _5)\right] b)\nonumber \\&\times (\bar{l}\gamma ^\mu l)+c_{10}(\bar{d}\gamma _\mu (1-\gamma _5)b) (\bar{l}\gamma ^\mu \gamma _5l)\Biggr \}. \end{aligned}$$
(20)

The values of the Wilson coefficients \(c_i\) and of the effective Wilson coefficient \(c_7^\mathrm{eff}\) are taken from Ref. [28]. The effective Wilson coefficient \( c_9^\mathrm{eff}\) contains additional perturbative and long-distance contributions. It can be written as

$$\begin{aligned} c_9^\mathrm{eff}&= c_9+h^\mathrm{eff}\left( \frac{m_c}{m_b},\frac{q^2}{m_b^2}\right) \nonumber \\&\times (3c_1+c_2+3c_3+c_4+3c_5+c_6)\nonumber \\&+\lambda _u\left[ h^\mathrm{eff}\left( \frac{m_c}{m_b},\frac{q^2}{m_b^2}\right) -h^\mathrm{eff}\left( \frac{m_u}{m_b},\frac{q^2}{m_b^2}\right) \right] \nonumber \\&\times \,(3c_1+c_2)- \frac{1}{2} h\left( 1,\frac{q^2}{m_b^2}\right) (4c_3+4c_4+3c_5+c_6)\nonumber \\&-\frac{1}{2} h\left( 0,\frac{q^2}{m_b^2}\right) (c_3\!+\!3c_4)\!+\!\frac{2}{9}(3c_3+c_4+3c_5+c_6),\nonumber \\ \end{aligned}$$
(21)

where \(\lambda _u=\frac{V_{ud}^*V_{ub}}{V_{td}^*V_{tb}}\) and

$$\begin{aligned}&\begin{aligned} h\left( \frac{m_c}{m_b}, \frac{q^2}{m_b}\right)&= \!-\! \frac{8}{9}\ln \frac{m_c}{m_b} \!+\! \frac{8}{27} \!+\! \frac{4}{9} x \!-\! \frac{2}{9}(2\!+\!x) |1\!-\!x|^{1/2} \nonumber \\&\quad \times \left\{ \begin{array}{ll} \ln \left| \frac{\sqrt{1-x} + 1}{\sqrt{1-x} - 1}\right| - i\pi , &{} x \equiv \frac{4 m_c^2}{ q^2} \!<\! 1, \\ 2 \arctan \frac{1}{\sqrt{x-1}}, &{} x \equiv \frac{4 m_c^2}{ q^2} \!>\! 1, \end{array} \right. \end{aligned} \\&h\left( 0, \frac{q^2}{m_b} \right) = \frac{8}{27} - \frac{4}{9} \ln \frac{q^2}{m_b} + \frac{4}{9} i\pi , \end{aligned}$$

while the function

$$\begin{aligned}&h^\mathrm{eff}\left( \frac{m_c}{m_b},\frac{q^2}{m_b^2}\right) =h\left( \frac{m_c}{m_b},\frac{q^2}{m_b^2}\right) \nonumber \\&\quad + \frac{3\pi }{\alpha ^2 c_0} \sum _{V_i=J/\psi ,\psi (2S)\dots } \frac{\Gamma (V_i\rightarrow l^+l^{-})M_{V_i}}{M_{V_i}^{2}-q^{2}-iM_{V_i}\Gamma _{V_i}} \end{aligned}$$
(22)

contains additional long-distance (nonperturbative) contributions which originate from the \(c\bar{c}\) mesons [\(J/\psi , \psi (2S)\dots \)]. We include contributions of the vector \(V_i(1^{--})\) charmonium states: \(J/\psi \), \(\psi (2S)\), \(\psi (3770)\), \(\psi (4040)\), \(\psi (4160)\) and \(\psi (4415)\), with their masses (\(M_{V_i}\)), leptonic [\(\Gamma (V_i\rightarrow l^+l^-)\)] and total (\(\Gamma _{V_i}\)) decay widths taken from PDG [26]. The coefficient \(c_0=3c_1+c_2+3c_3+c_4+3c_5+c_6\). Similar expression holds for the function \(h^\mathrm{eff}\small \left( \frac{m_u}{m_b},\frac{q^2}{m_b^2}\right) \), where the long-distance contributions now come from \(u\bar{u}\) states (\(\rho \) and \(\omega \)).

The matrix element of the \(b\rightarrow d l^+l^-\) transition amplitude between meson states can be expressed through the helicity amplitudes \(H^{(i)}_m\) (where the superscript \(i=1,2\) corresponds to the first and second terms in the amplitude (20), while the subscript \(m=\pm ,0,t\) denotes transverse, longitudinal, and time helicity components, respectively). The explicit formulas for the helicity amplitudes in terms of the decay form factors defined in Eqs. (1)–(7) are given in our papers [5, 6].

Then the differential decay rate can be written in terms of the helicity amplitudes [5, 29] as follows:

  1. (a)

    for the \(B\rightarrow \pi (\rho )l^+l^-\) decays

    $$\begin{aligned}&\frac{d\Gamma (B\rightarrow \pi (\rho )l^+l^-)}{dq^2}=\frac{G_F^2}{(2\pi )^3}\left( \frac{\alpha |V_{td}^*V_{tb}|}{2\pi }\right) ^2\nonumber \\&\qquad \times \, \frac{\lambda ^{1/2}q^2}{48M_{B}^3} \sqrt{1-\frac{4m_l^2}{q^2}}\Biggl [H^{(1)}H^{\dag (1)}\left( 1+\frac{2m_l^2}{q^2}\right) \nonumber \\&\qquad +\, H^{(2)}H^{\dag (2)}\left( 1-\frac{4m_l^2}{q^2}\right) +\frac{2m_l^2}{q^2}3 H^{(2)}_tH^{\dag (2)}_t\Biggr ],\nonumber \\ \end{aligned}$$
    (23)
  2. (b)

    for the \(B\rightarrow \pi (\rho )\nu \bar{\nu }\) decays

    $$\begin{aligned} \frac{d\Gamma (B\rightarrow \pi (\rho )\nu \bar{\nu })}{dq^2}&= 3\frac{G_F^2}{(2\pi )^3} \left( \frac{\alpha |V_{td}^*V_{tb}|}{2\pi }\right) ^2 \frac{\lambda ^{1/2}q^2}{24M_{B}^3} \nonumber \\&\times H^{(\nu )}H^{\dag (\nu )}, \end{aligned}$$
    (24)

    where \(m_l\) is the lepton mass and

    $$\begin{aligned} H^{(i)}H^{\dag (i)}\equiv H^{(i)}_+H^{\dag (i)}_++H^{(i)}_-H^{\dag (i)}_-+H^{(i)}_0H^{\dag (i)}_0.\quad \end{aligned}$$
    (25)

The forward–backward asymmetry for the \(B\rightarrow \rho \mu ^+\mu ^-\) decay can be expressed in terms of the helicity amplitudes in the following way:

$$\begin{aligned} A_{ FB}\!&= \!\frac{3}{4} \sqrt{1-\frac{4m_l^2}{q^2}} \nonumber \\&\times \frac{ \mathrm{Re}(H^{(1)}_+H^{\dag (2)}_+)-\mathrm{Re}(H^{(1)}_-H^{\dag (2)}_-)}{H^{(1)}H^{\dag (1)}\left( 1\!+\!\frac{2m_l^2}{q^2}\right) \!+\! H^{(2)}H^{\dag (2)}\left( 1\!-\!\frac{4m_l^2}{q^2}\right) \!+\!\frac{2m_l^2}{q^2}3 H^{(2)}_tH^{\dag (2)}_t},\nonumber \\ \end{aligned}$$
(26)

while the longitudinal polarization fraction of the vector \(\rho \) meson is given by

$$\begin{aligned} F_L\!&= \!\frac{H^{(1)}_0H^{\dag (1)}_0\left( 1\!+\!\frac{2m_l^2}{q^2}\right) \!+\! H^{(2)}_0H^{\dag (2)}_0\left( 1\!-\!\frac{4m_l^2}{q^2}\right) \!+\!\frac{2m_l^2}{q^2}3 H^{(2)}_tH^{\dag (2)}_t }{H^{(1)}H^{\dag (1)}\left( 1\!+\!\frac{2m_l^2}{q^2}\right) \!+\! H^{(2)}H^{\dag (2)}\left( 1\!-\!\frac{4m_l^2}{q^2}\right) \!+\!\frac{2m_l^2}{q^2}3 H^{(2)}_tH^{\dag (2)}_t}.\nonumber \\ \end{aligned}$$
(27)

Substituting the form factors of the \(B\rightarrow \pi \) and \(B\rightarrow \rho \) weak transitions calculated in Sect. 2 in the above expressions we get predictions for the differential decay rates, the forward–backward asymmetry and longitudinal polarization fraction of the vector \(\rho \) meson. The obtained differential distributions for the \(B^+\rightarrow \pi ^+\mu ^+\mu ^-(\tau ^+\tau ^-)\) and \(B^+\rightarrow \rho ^+\mu ^+\mu ^-(\tau ^+\tau ^-)\) decays are plotted in Figs. 6, 7, 8, and 9. The dashed and solid lines in these figures correspond to the so-called resonant and nonresonant results which were obtained with and without inclusion of the long-distance contributions originating from the \(c\bar{c}\) and \(u\bar{u}\) resonances [see Eq. (22)] in the effective coefficient \(c_9^\mathrm{eff}\) (21). The regions of the highest \(J/\psi \) and \(\psi (2S)\) peaks are usually vetoed in experiment in order to resolve the signal against their huge background. Other asymmetries both time-independent and time-dependent in these decays are discussed in detail in Ref. [30].

Fig. 6
figure 6

Theoretical predictions for the differential branching fractions \(d Br(B^+ \rightarrow \pi ^+ l^+l^-)/d q^2\). Nonresonant and resonant results are plotted by solid and dashed lines, respectively

Fig. 7
figure 7

Same as in Fig. 6 but for \(d Br(B \rightarrow \rho l^+l^-)/d q^2\) decay

Fig. 8
figure 8

Same as in Fig. 6 but for the longitudinal polarization \(F_L\) (left) and muon forward–backward asymmetry \(A_{FB}\) (right) for the rare \(B^+ \rightarrow \rho ^+ \mu ^+\mu ^-\) decay

Fig. 9
figure 9

Theoretical predictions for the differential branching fractions \(d Br(B^+ \rightarrow \pi ^+\nu \bar{\nu })/d q^2\) (left) and \(d Br(B^+ \rightarrow \rho ^+\nu \bar{\nu })/d q^2\) (right) (in \(10^{-6}\))

In Table 4 we present our predictions for the differential branching fractions of the rare semileptonic \(B^+\rightarrow \pi ^+ \mu ^+\mu ^-\) and \(B^+\rightarrow \rho ^+ \mu ^+\mu ^-\) decays integrated over several bins of \(q^2\) which in principle can be measured experimentally. In this table we also give the recent theoretical estimates [22] for the \(B^+\rightarrow \pi ^+ \mu ^+\mu ^-\) decay which are based on the vector weak current form factors [\(f_+(q^2)\) and \(f_0(q^2)\)] extracted from the combined analysis of the available experimental data. For the tensor current form factor [\(f_T(q^2)\)] lattice QCD results for the \(B\rightarrow K\) transition and \(SU(3)_F\)-breaking Ansatz were used. We find that in most \(q^2\) bins theoretical predictions agree within error bars.

Integrating the differential branching fractions over \(q^2\) we get the total rare decay branching fractions for \(B^+\rightarrow \pi ^+ l^+l^-(\nu \bar{\nu })\) and \(B^+\rightarrow \rho ^+ l^+l^-(\nu \bar{\nu })\). In Table 5 we compare our results for the nonresonant branching fractions with other theoretical calculations. At present experimental data are only available for the \(B^+\rightarrow \pi ^+ \mu ^+\mu ^-\) decay [4], which was observed recently. We see that theoretical predictions agree with each other and the experimental value within errors.

Table 4 Comparison of theoretical predictions for the branching fractions of the rare semileptonic \(B^+\rightarrow \pi ^+ \mu ^+\mu ^- \) and \(B^+\rightarrow \rho ^+ \mu ^+\mu ^- \) decays in several bins of \(q^2\) (in \(10^{-8}\))
Table 5 Theoretical predictions for the nonresonant branching fractions of the rare semileptonic \(B\) decays and available experimental data (in \(10^{-8}\))

5 Conclusions

The form factors parametrizing the heavy-to-light \(B\rightarrow \pi \) and \(B\rightarrow \rho \) weak transition matrix elements were obtained in the framework of the relativistic quark model. The quasipotential approach was used to express these form factors through the overlap integrals of the initial and final meson wave functions, which are taken from the previous calculations of meson masses. All relativistic effects, including the wave function transformations from the rest to the moving reference frame as well as contributions of the intermediate negative-energy states, were consistently taken into account. Our approach allowed us to explicitly determine the form factor dependence on the momentum transfer \(q^2\) in the whole kinematical range without additional model assumptions or extrapolations.

First we confronted the predictions of our model for the differential branching fractions of the semileptonic \(B\rightarrow \pi l\nu _l\) and \(B\rightarrow \rho l\nu _l\) decays with recent detailed experimental data [13]. Good agreement for all observables was found. From this comparison we determined the exclusive value of the CKM matrix element \(|V_{ub}|=(4.15\,{\pm }\,0.09_\mathrm{exp} \pm 0.21_\mathrm{theor})\times 10^{-3}\), which is consistent with the one extracted from the inclusive semileptonic \(B\rightarrow X_ul\nu _l\) decays [26].

Then we considered the rare weak \(B^+\rightarrow \pi ^+ l^+l^-(\nu \bar{\nu })\) and \(B^+\rightarrow \rho ^+ l^+l^-(\nu \bar{\nu })\) decays. Calculations were done both with and without account of the long-range contributions of the heavy charmonium states and light \(\rho ,\omega \) resonances. Detailed predictions for the differential branching fractions of these decays were presented. The calculated total branching fraction for the rare decay \(B^+\rightarrow \pi ^+ \mu ^+\mu ^-\) agrees well with the recent measurement [4]. The LHCb Collaboration also measured the ratio of the \(B^+\rightarrow \pi ^+ \mu ^+\mu ^-\) and \(B^+\rightarrow K^+ \mu ^+\mu ^-\) branching fractions to be \(0.053\pm 0.014\pm 0.001\). Using our prediction for the \(B^+\rightarrow K^+ \mu ^+\mu ^-\) decay [5] we get the value of this ratio equal to \(0.048\pm 0.005\) which agrees with the experimental one within error bars. The ratio of the corresponding branching fractions involving vector \(\rho \) and \(K^*\) mesons is predicted to be \(Br(B^+\rightarrow \rho ^+ \mu ^+\mu ^-)/Br(B^+\rightarrow K^{*+}\mu ^+\mu ^-)=0.048\pm 0.005\).