# Theoretical analysis of the leptonic decays
$\bm{B\to\ell\ell\ell^{\prime}\bar{\nu}_{\ell^{\prime}}}$

Mikhail A. Ivanova and Dmitri Melikhova,b,c aJoint Institute for Nuclear
Research, Dubna, 141980, Russia
bD. V. Skobeltsyn Institute of Nuclear Physics, M. V. Lomonosov Moscow State
University, 119991, Moscow, Russia
cFaculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090 Vienna,
Austria

###### Abstract

We discuss the general properties of the amplitude of the $B\to
l^{+}l^{-}l^{\prime}\nu^{\prime}$ decays and calculate the related kinematical
distributions $d^{2}\Gamma/dq^{2}dq^{\prime 2}$, $q$ the momentum of the
$l^{+}l^{-}$ pair emitted from the electromagnetic vertex and $q^{\prime}$ the
momentum of the $l^{\prime}\nu^{\prime}$ pair emitted from the weak vertex. We
emphasize that electromagnetic gauge invariance imposes essential constraints
on the $B\to\gamma^{*}l^{\prime}\nu^{\prime}$ amplitude at small $q^{2}$ which
in the end yield the behaviour of the differential branching fraction as
$d\Gamma(B\to l^{+}l^{-}l^{\prime}\nu^{\prime})/dq^{2}\propto 1/q^{2}$ and a
mild logarithmic dependence of $\Gamma(B\to l^{+}l^{-}l^{\prime}\nu^{\prime})$
on the lepton mass $m_{l}$. Consequently, (i) the main contribution to the
decay rate $\Gamma(B\to\mu^{+}\mu^{-}e\nu_{e})$ comes from the region of light
vector resonances $\rho^{0}$ and $\omega$, $q^{2}\simeq
M_{\rho}^{2},M_{\omega}^{2}$ and (ii) the decay rate $\Gamma(B\to
e^{+}e^{-}\mu\nu_{\mu})$ receives comparable contributions from the region of
small $q^{2}$ and from the resonance region. As the result, the decay rate
$\Gamma(B\to e^{+}e^{-}\mu\nu_{\mu})$ is only a factor $\sim 2$ larger than
$\Gamma(B\to\mu^{+}\mu^{-}e\nu_{e})$. We perform a detailed analysis of the
uncertainties in the theoretical predictions for the decays $B\to
l^{+}l^{-}l^{\prime}\nu^{\prime}$ in the Standard Model. We found that the
theoretical expectations for such decays in the Standard Model are only
marginally compatible with the recent upper limits of the LHCb collaboration
exp1 .

##  1 Introduction

This paper aims at studying the general properties of the amplitude
$B\to\gamma^{*}l^{\prime}\nu^{\prime}$ with the emphasis on the constraints
imposed by the electromagnetic gauge invariance, and, eventually, obtaining
predictions for the differential distributions in the $B$-meson decays into
four leptons in the final state, $B\to l^{+}l^{-}l^{\prime}\nu^{\prime}$. The
latter reactions are being studied experimentally exp1 ; exp2 ; exp3 ; exp4 ,
thus requiring a proper theoretical understanding of the $B$-meson form
factors into two currents. By now, there have been a few theoretical papers
sehgal ; nikitin ; beneke2 , where $B$-decays into two lepton pairs have been
studied. Here we emphasize that essential constraints on the
$B\to\gamma^{*}l^{\prime}\nu^{\prime}$ amplitude imposed by the
electromagnetic gauge invariance escaped the attention of the previous
analyses. We demonstrate that an essential consequence of this constraint is
the behaviour of the differential distribution at small $q^{2}$ as

$\displaystyle d\Gamma(B\to l^{+}l^{-}l^{\prime}\nu^{\prime})/dq^{2}\propto
1/q^{2},$ (1.1)

and not the behaviour $\propto 1/q^{4}$ as may be obtained naively nikitin ;
beneke2  if the gauge-invariance constraints are not taken into account. This
essential qualitative feature of the $q^{2}$-distribution at small $q^{2}$
leads to important consequences for the theoretical estimates of the $B\to
l^{+}l^{-}l^{\prime}\nu^{\prime}$ branching fractions: (i) The main
contribution to ${\rm Br}(B\to\mu^{+}\mu^{-}e\nu_{e})$ comes from the region
of $q^{2}$ around the light vector resonances whereas the region of small
$q^{2}$ gives a comparatively smaller contribution; (ii) The branching
fraction of ${\rm Br}(B\to l^{+}l^{-}l^{\prime}\nu^{\prime})$ has a mild
logarithmic dependence $\propto\log(m_{l}^{2})$ on the lepton mass $m_{l}$
(one would obtain a much stronger $1/m_{l}^{2}$ dependence if the gauge-
invariance constraint is not taken into account), and as the result ${\rm
Br}(B\to e^{+}e^{-}\mu\nu_{\mu})$ turns out to be only a factor $\sim 2$
larger than ${\rm Br}(B\to\mu^{+}\mu^{-}e\nu_{e})$ in contrast to the
estimates presented in nikitin .

The $B\to\gamma^{*}l^{\prime}\nu^{\prime}$ amplitude (see Fig. 1) may be
parameterized via Lorenz-invariant form factors as follows:

$\displaystyle T_{\alpha\nu}(q,q^{\prime}|p)=i\int dx\,e^{iqx}\langle
0|T\\{j^{\rm e.m.}_{\alpha}(x),\bar{u}(0){\cal
O_{\nu}}b(0))\\}|B(p)\rangle=\sum_{i}L^{(i)}_{\alpha\nu}(q,q^{\prime})F_{i}(q^{\prime
2},q^{2})+\ldots,\quad p=q+q^{\prime},$ (1.2)

with $q^{\prime}$ the momentum of the weak $b\to u$ current, and $q$ the
momentum of the electromagnetic current. In Eq. (1.2), ${\cal
O_{\nu}}=\gamma_{\nu},\gamma_{\nu}\gamma_{5}$ and $j_{\alpha}^{\rm e.m.}$ is
the conserved electromagnetic current

$\displaystyle j_{\alpha}^{\rm
e.m.}(0)=eQ_{b}\bar{b}(0)\gamma_{\alpha}b(0)+eQ_{u}\bar{u}(0)\gamma_{\alpha}u(0).$
(1.3)

The quantities $L^{(i)}_{\alpha\nu}(q,q^{\prime})$ represent the transverse
Lorents structures, $q^{\alpha}L^{(i)}_{\alpha\nu}(q,q^{\prime})=0$, and the
dots stand for the longitudinal part which is fully determined by the
conservation of the electromagnetic current, $\partial_{\alpha}j_{\alpha}^{\rm
e.m.}=0$, and the equal-time commutation relations.

The form factors $F_{i}(q^{\prime 2},q^{2})$ are complicated functions of two
variables, $q^{\prime 2}$ and $q^{2}$; the general properties of these objects
in QCD have been studied recently in ims . Noteworthy, gauge invariance
provides essential constraints on some of the form factors describing the
transition of the $B$-meson into the real photon, i.e., at $q^{2}=0$ bmns2001
; m2002 ; kruger ; kmn2016 .

In the past, theoretical analyses focused on a family of similar reactions,
namely, the $B\to\gamma l^{+}l^{-}$ and $B\to\gamma l\nu$ decays (see, e.g.,
aliev ; korchemsky ; kou ; mn2004 ; beneke1 ; kmn2018 ; ivanov ); these
processes are described by the same form factors as four-lepton $B$-decays,
but evaluated at a zero value of one of the momenta squared. The corresponding
form factors depend on one variable $q^{\prime 2}$, $q^{\prime}$ the momentum
of the weak current; for instance, for radiative leptonic decays $B\to\gamma
l^{\prime}\nu^{\prime}$, one needs the form factors $F_{i}(q^{\prime
2},q^{2}=0)$.

The four-lepton decay of interest, $B\to l^{+}l^{-}l^{\prime}\nu^{\prime}$,
requires the form factors $F_{i}(q^{\prime 2},q^{2})$ for $0<q^{2},q^{\prime
2}<M_{B}^{2}$. The dependence of the form factors on the variable $q^{\prime
2}$ can be predicted reasonably well: there are no hadron resonances in the
full decay region $0<q^{\prime 2}<M_{B}^{2}$, and the $q^{\prime
2}$-dependence of the form factors is determined to a large extent by the
influence of the beauty mesons with the appropriate quantum numbers; all these
mesons are heavier than the $B$-meson and therefore lie beyond the physical
decay region of the variable $q^{\prime 2}$. The calculation of the
$q^{2}$-dependence of the form factors is a much more difficult task: light
vector mesons lie in the physical decay region and should be properly taken
into account. At $q^{2}$ in the region of light vector meson resonances, which
turns out to give the dominant contribution to the four-meson decay rate, the
form factors cannot be obtained directly in pQCD ims ; here model
considerations based on the explicit account of these light vector
resonances—including their finite width effects—are mandatory. We discuss in
great detail form factors in the full kinematical region of four-lepton
$B$-decays and obtain the corresponding differential distributions. Still, any
consistent model for the form factors should obey the gauge-invariance
constraints. We will demonstrate that such constraints determine some
essential features of the differential distributions in the region of small
values of $q^{2}$.

Fig. 1:  Feynman diagrams describing the amplitude (1.2).

##  2 Constraints on the transition form factors

We now discuss the requirements imposed by the electromagnetic gauge
invariance on the $B\to\gamma^{*}$ transition amplitudes
$\langle\gamma^{*}(q)|\bar{u}\left\\{\gamma_{\nu},\gamma_{\nu}\gamma_{5}\right\\}b|\bar{B}_{u}(p)\rangle$
induced by the vector and the axial-vector charged currents. The corresponding
form factors are functions of two variables, $q^{2}$ and $q^{\prime 2}$, where
$q^{\prime}$ is the momentum of the weak $b\to u$ current, and $q$ is the
momentum of the electromagnetic current, and $p=q+q^{\prime}$. Gauge
invariance provides constraints on some of the form factors describing the
transition of $B_{u}$ to the real photon, $q^{2}=0$.

###  2.1 Form factors of the vector weak current

In case of the vector charged quark current $\bar{u}\gamma_{\nu}b$, the gauge-
invariant amplitude contains one Lorentz structure and one dimensionless form
factor $F_{V}(q^{\prime 2},q^{2})$:

$\displaystyle T_{\alpha,\nu}=i\int dxe^{iqx}\langle 0|T\left\\{j^{\rm
e.m.}_{\alpha}(x),\bar{u}\gamma_{\nu}b(0)\right\\}|\bar{B}_{u}(p)\rangle=e\,\epsilon_{\mu\alpha
q^{\prime}q}\frac{F_{V}(q^{\prime 2},q^{2})}{M_{B}}.$ (2.4)

The amplitude is transverse, $q^{\alpha}T_{\alpha,\nu}=0$, and contains no
contact term. It is free of the kinematical singularities so gauge invariance
provides no constraints on $F_{V}(q^{\prime 2},q^{2}=0)$. The contribution of
the vector charged quark current to the amplitude of the
$B\to\gamma^{*}l^{\prime}\nu^{\prime}$ decay reads

$\displaystyle A_{\rm
vector}(B\to\gamma^{*}l^{\prime}\nu^{\prime})=e\frac{G_{F}}{\sqrt{2}}\,V_{ub}\,\bar{l}^{\prime}\gamma_{\nu}(1-\gamma_{5})\nu^{\prime}\,\varepsilon_{\alpha}^{*}(q)\epsilon_{\nu\alpha
q^{\prime}q}\frac{F_{V}(q^{\prime 2},q^{2})}{M_{B}}.$ (2.5)

###  2.2 Form factors of the axial-vector weak current

For the axial-vector current, $\bar{u}\gamma_{\nu}\gamma_{5}b$, the
corresponding amplitude is more complicated: it contains three independent
gauge-invariant structures and three form factors, $f_{1A}(q^{\prime
2},q^{2}),f_{2A}(q^{\prime 2},q^{2}),f_{3A}(q^{\prime 2},q^{2})$, and in
addition it has the contact term which is fully determined by the conservation
of the electromagnetic current, $\partial^{\alpha}j_{\alpha}^{e.m.}=0$:

$\displaystyle T^{5}_{\alpha,\nu}$ $\displaystyle=$ $\displaystyle i\int
dxe^{iqx}\langle 0|T\left\\{j^{\rm
e.m.}_{\alpha}(x),\bar{u}\gamma_{\nu}\gamma_{5}b(0)\right\\}|\bar{B}_{u}(p)\rangle$
(2.6) $\displaystyle=$ $\displaystyle
ie\,\left(g_{\mu\alpha}-\frac{q_{\alpha}q_{\mu}}{q^{2}}\right)\,q^{\prime}q\,f_{1A}+ie\,\left(q^{\prime}_{\alpha}-\frac{q^{\prime}q}{q^{2}}q_{\alpha}\right)\bigg{\\{}p_{\nu}f_{2A}+q_{\mu}f_{3A}\bigg{\\}}+ie\,Q_{B}f_{B}\frac{q_{\alpha}p_{\mu}}{q^{2}}.$

Here $Q_{B}\equiv Q_{\bar{B}_{u}}=Q_{b}-Q_{u}$ is the electric charge of the
$\bar{B}_{u}$ meson and $f_{B}>0$ is defined according to

$\displaystyle\langle
0|\bar{u}\gamma_{\mu}\gamma_{5}b|\bar{B}_{u}(p)\rangle=if_{B}p_{\mu}.$ (2.7)

The last term in (2.6) is just the longitudinal contact term mentioned above.
Let us briefly recall the standard way this term is obtained: We calculate
$q_{\alpha}T^{5}_{\alpha,\nu}$, represent
$q_{\alpha}e^{iqx}=-i\frac{\partial}{\partial x^{\alpha}}e^{iqx}$, and perform
parts integration moving the derivative to the $T$-product. Making use of the
conservation of the electromagnetic current
$\partial^{\alpha}j^{e.m.}_{\alpha}=0$, the only nonzero contribution comes
from the differentiation of the $\theta$-functions defining the $T$-product,
leading to the equal-time commutator. In the end, we obtain ($\hat{Q}$ is the
time-independent electric charge operator, $\hat{Q}=\int
d^{3}x\,j_{0}(x^{0},\vec{x})$)

$\displaystyle q_{\alpha}T^{5}_{\alpha,\nu}=-\langle
0|[\hat{Q},\bar{u}\gamma_{\mu}\gamma_{5}b(0)]|B(p)\rangle=iQ_{B}f_{B}p_{\nu}.$
(2.8)

This relation does not determine the longitudinal Lorentz structure in the
unique way: one can, e.g., choose this structure in the form
$p_{\alpha}p_{\nu}/pq$ or in the form $q_{\alpha}p_{\nu}/q^{2}$. However, only
the latter form, which is implemented in (2.6), corresponds to the
longitudinal part in the form of a contact term.111By definition, a contact
term is a quantity represented by $\delta$-function and its derivatives in
configuration space; therefore $q_{\alpha}p_{\nu}/q^{2}$ is a contact term
whereas $p_{\alpha}p_{\nu}/pq$ is not a contact term according to the standard
definition. For further details we refer to currentalgebra . Obviously,
different choices of the longitudinal part lead to redefinitions of the form
factors $f_{iA}$ in the transverse part of the amplitude. The choice of the
longitudinal structure in the form of a contact term $q_{\alpha}p_{\nu}/q^{2}$
is suggested by the structure of the quark electromagnetic vertex and is
preferable with respect to the analytic properties of the form factors
$f_{iA}$ m2002 .

The projectors in (2.6) contain kinematical singularities at $q^{2}=0$; these
kinematical singularities however should not be the singularities of the
physical amplitude, as the spectrum of physical states does not contain a
massless particle in the $q^{2}$-channel. Therefore, as the consequence of
gauge invariance, we obtain the following relations between the form factors
at $q^{2}=0$:

$\displaystyle\left[f_{1A}(q^{\prime 2},q^{2})+f_{3A}(q^{\prime
2},q^{2})\right]_{q^{2}=0}=0,$ (2.9)
$\displaystyle[q^{\prime}q\,f_{2A}(q^{\prime 2},q^{2})]_{q^{2}=0}=Q_{B}f_{B}.$
(2.10)

To implement these constraints at $q^{2}=0$, we write down dispersion
representations for the form factors $f_{1A},f_{2A},f_{3A}$ in the variable
$q^{2}$ with one subtraction and determine the subtraction terms to satisfy
(2.9) and (2.10). Such representations have the following form

$\displaystyle f_{1A}(q^{\prime 2},q^{2})$ $\displaystyle=$
$\displaystyle\xi(q^{\prime 2})+q^{2}\int\frac{ds}{\pi
s(s-q^{2})}\rho_{1A}(q^{\prime 2},s),$ (2.11) $\displaystyle f_{2A}(q^{\prime
2},q^{2})$ $\displaystyle=$
$\displaystyle\frac{2Q_{B}f_{B}}{M_{B}^{2}-q^{\prime
2}}+q^{2}\int\frac{ds}{\pi s(s-q^{2})}\rho_{2A}(q^{\prime 2},s),$ (2.12)
$\displaystyle f_{3A}(q^{\prime 2},q^{2})$ $\displaystyle=$
$\displaystyle-\xi(q^{\prime 2})+q^{2}\int\frac{ds}{\pi
s(s-q^{2})}\rho_{3A}(q^{\prime 2},s).$ (2.13)

The form factor $\xi(q^{\prime 2})$ is related to the form factor of the
$B\to\gamma l^{\prime}\nu^{\prime}$ transition, and for the spectral densities
$\rho_{Ai}(q^{\prime 2},s)$ we will construct phenomenological expressions
based on the contributions of the light vector resonances $\rho_{0}$ and
$\omega$.

Next, we should add the Bremsstrahlung contribution (i.e., the photon emitted
from the lepton $l^{\prime}$ in the final state) that in the limit of a
massless lepton $m_{l^{\prime}}=0$ reads

$\displaystyle A_{\rm
Brems}(B\to\gamma^{*}l^{\prime}\nu^{\prime})=ie\,Q_{l}\frac{G_{F}}{\sqrt{2}}\,V_{ub}\,\bar{l}^{\prime}\gamma_{\nu}(1-\gamma_{5})\nu^{\prime}\,\varepsilon^{*}_{\alpha}(q)f_{B}\left(-g_{\alpha\nu}\right),\qquad
Q_{l}=Q_{B},$ (2.14)

The axial part of weak-transition amplitude
$B\to\gamma^{*}l^{\prime}\nu^{\prime}$ then takes the form

$\displaystyle A_{\rm axial}(B\to\gamma^{*}l^{\prime}\nu^{\prime})$ (2.15)
$\displaystyle=$ $\displaystyle
ie\frac{G_{F}}{\sqrt{2}}\,V_{ub}\,\bar{l}^{\prime}\gamma_{\nu}(1-\gamma_{5})\nu^{\prime}\,\varepsilon_{\alpha}^{*}(q)\left\\{\left(g_{\alpha\nu}-\frac{q_{\alpha}q_{\nu}}{q^{2}}\right)\,q^{\prime}q\,f_{1A}+\left(q^{\prime}_{\alpha}-\frac{q^{\prime}q}{q^{2}}q_{\alpha}\right)\bigg{[}p_{\nu}f_{2A}+q_{\nu}f_{3A}\bigg{]}+Q_{B}\,f_{B}\frac{q_{\alpha}p_{\nu}}{q^{2}}\right\\}$
$\displaystyle+ie\,Q_{B}\frac{G_{F}}{\sqrt{2}}\,V_{ub}\,\bar{l}^{\prime}\gamma_{\nu}(1-\gamma_{5})\nu^{\prime}\,\varepsilon^{*}_{\alpha}(q)f_{B}\left(-g_{\alpha\nu}\right),$

the last term being the Bremsstrahlung contribution (2.14).

The amplitude may be simplified by taking into account that
$q^{\alpha}\varepsilon^{*}_{\alpha}(q)=0$, yielding

$\displaystyle A_{\rm axial}(B\to\gamma^{*}l^{\prime}\nu^{\prime})$
$\displaystyle=$ $\displaystyle
ie\frac{G_{F}}{\sqrt{2}}V_{ub}\bar{l}^{\prime}\gamma_{\nu}(1-\gamma_{5})\nu^{\prime}\varepsilon_{\alpha}^{*}(q)$
(2.16)
$\displaystyle\times\bigg{\\{}(g_{\alpha\nu}\,q^{\prime}q-q^{\prime}_{\alpha}q_{\nu})\left[f_{1A}-\frac{Q_{B}f_{B}}{q^{\prime}q}\right]+q^{\prime}_{\alpha}q_{\nu}\left[f_{2A}+f_{3A}+f_{1A}-\frac{Q_{B}f_{B}}{q^{\prime}q}\right]+q^{\prime}_{\alpha}q^{\prime}_{\nu}f_{2A}\bigg{\\}}.$

For the massless lepton, $m_{l}^{\prime}=0$, the Lorentz structure
$q^{\prime}_{\alpha}q^{\prime}_{\nu}$ does not contribute and we find the
final expression for the contribution of the axial-vector part of the quark
current, $\bar{u}\gamma_{\nu}\gamma_{5}b$, to the amplitude

$\displaystyle A_{\rm
axial}(B\to\gamma^{*}l^{\prime}\nu^{\prime})=ie\frac{G_{F}}{\sqrt{2}}\,V_{ub}\bar{l}^{\prime}\gamma_{\nu}(1-\gamma_{5})\nu^{\prime}\varepsilon_{\alpha}^{*}(q)\left\\{(g_{\alpha\nu}q^{\prime}q-q^{\prime}_{\alpha}q_{\nu})\frac{F_{1A}}{M_{B}}+q^{\prime}_{\alpha}q_{\nu}\frac{F_{2A}}{M_{B}}\right\\},$
(2.17)

where dimensionless form factors $F_{1A,2A}$ are introduced

$\displaystyle\frac{F_{1A}(q^{\prime 2},q^{2})}{M_{B}}$ $\displaystyle=$
$\displaystyle f_{1A}(q^{\prime 2},q^{2})-\frac{Q_{B}f_{B}}{q^{\prime}q},$
$\displaystyle\frac{F_{2A}(q^{\prime 2},q^{2})}{M_{B}}$ $\displaystyle=$
$\displaystyle f_{2A}(q^{\prime 2},q^{2})+f_{3A}(q^{\prime
2},q^{2})+f_{1A}(q^{\prime 2},q^{2})-\frac{Q_{B}f_{B}}{q^{\prime}q}.$ (2.18)

Notice that as follows from Eqs. (2.9) and (2.10), the form factor $F_{2A}$
vanishes for the real photon in the final state:

$\displaystyle F_{2A}(q^{\prime 2},q^{2}=0)=0.$ (2.19)

The latter relation is important as it has strong impact on the behavior of
the differential distributions in $B\to lll^{\prime}\nu^{\prime}$ at small
$q^{2}$.

##  3 The $B\to\gamma l^{\prime}\nu^{\prime}$ transition

We now illustrate the way the well-known formulas for the amplitude and the
differential distribution in the leptonic radiative $B$ decay, $B\to\gamma
l^{\prime}\nu^{\prime}$, emerge.

Since, as the consequence of gauge invariance, $F_{2A}(q^{\prime
2},q^{2}=0)=0$, only the form factors $F_{1A}$ and $F_{V}$ contribute to the
amplitude for the real photon and the massless lepton in the final state, and
one finds the $B\to\gamma l^{\prime}\nu^{\prime}$ amplitude:

$\displaystyle A(B\to\gamma
l^{\prime}\nu^{\prime})=ie\frac{G_{F}}{\sqrt{2}}\,V_{ub}\bar{l}^{\prime}\gamma_{\nu}(1-\gamma_{5})\nu^{\prime}\varepsilon_{\alpha}^{*}(q)\bigg{\\{}(g_{\alpha\nu}q^{\prime}q-q^{\prime}_{\alpha}q_{\nu})\frac{F_{A}(q^{\prime
2})}{M_{B}}+i\epsilon_{\nu\alpha q^{\prime}q}\frac{F_{V}(q^{\prime
2})}{M_{B}}\bigg{\\}},$ (3.20)

where $F_{A}(q^{\prime 2})\equiv F_{1A}(q^{\prime 2},q^{2}=0)$ and
$F_{V}(q^{\prime 2})\equiv F_{V}(q^{\prime 2},q^{2}=0)$.

The differential decay rate (for the massless lepton $m_{l^{\prime}}=0$) takes
a simple form:

$\displaystyle\frac{d\Gamma(B\to\gamma
l^{\prime}\nu^{\prime})}{dE_{\gamma}}=\frac{G_{F}^{2}V_{ub}^{2}}{48\pi^{2}}M_{B}^{4}\alpha_{\rm
e.m.}\,x_{\gamma}^{3}(1-x_{\gamma})\big{(}|F_{A}|^{2}+|F_{V}|^{2}\big{)},\quad
x_{\gamma}=2E_{\gamma}/M_{B},$ (3.21)

where $M_{B}E_{\gamma}=pq=q^{\prime}q$, $E_{\gamma}$ the photon energy in the
$B$-meson rest frame.

##  4 The $B^{-}\to l^{+}l^{-}l^{\prime-}\bar{\nu}^{\prime}$ transition

The amplitude of the $B\to lll^{\prime}\nu^{\prime}$ transition is readily
obtained from the amplitudes of the $B\to\gamma^{*}l^{\prime}\nu^{\prime}$
transitions induced by the vector and the axial quark current by performing
the replacement

$\displaystyle\varepsilon^{*}_{\alpha}(q)\to-e\frac{\bar{l}\gamma_{\alpha}l}{q^{2}},$
(4.22)

so we obtain

$\displaystyle A(B\to
lll^{\prime}\nu^{\prime})=ie^{2}\frac{G_{F}}{\sqrt{2}}\,V_{ub}\cdot\bar{l}\gamma_{\alpha}l\cdot\bar{l}^{\prime}\gamma_{\nu}(1-\gamma_{5})\nu^{\prime}\frac{1}{q^{2}}\bigg{\\{}(g_{\alpha\nu}q^{\prime}q-q^{\prime}_{\alpha}q_{\nu})\frac{F_{1A}}{M_{B}}+q^{\prime}_{\alpha}q_{\nu}\frac{F_{2A}}{M_{B}}+i\epsilon_{\nu\alpha
q^{\prime}q}\frac{F_{V}}{M_{B}}\bigg{\\}},$ (4.23)

where the dimensionless form factors $F_{1A,2A,V}$ are functions of $q^{2}$
and $q^{\prime 2}$, and $F_{2A}(q^{\prime 2},q^{2}=0)=0$. As the next step, by
summing over the lepton polarizations and integrating over the phase space of
the $l^{+}l^{-}$ pair and the $l^{\prime}\nu^{\prime}$ pair, one easily
obtains the explicit analytic expression for the double differential
distribution.

$\displaystyle|A|^{2}=\frac{G_{F}^{2}}{2}\,V_{ub}^{2}\,\frac{4}{3}\frac{8}{3}\frac{e^{4}}{(q^{2})^{2}}\frac{1}{M_{B}^{2}}$
$\displaystyle\bigg{[}$
$\displaystyle(F_{1A}F^{*}_{1A}+F_{V}F_{V}^{*})\left\\{2(q^{\prime}q)^{2}+q^{2}q^{\prime
2}\right\\}q^{2}q^{\prime 2}$ (4.24)
$\displaystyle+(F_{1A}F_{2A}^{*}+F_{1A}^{*}F_{2A})\left\\{(q^{\prime}q)^{2}-q^{\prime
2}q^{2}\right\\}q^{2}q^{\prime 2}$
$\displaystyle+F_{2A}F_{2A}^{*}\left\\{(q^{\prime 2}q^{2})^{2}-2q^{\prime
2}q^{2}(q^{\prime}q)^{2}+(q^{\prime}q)^{4}\right\\}\bigg{]}.$

The factors $4/3$ and $8/3$ result from the summation over polarizations of
massless leptons coming from the electromagnetic vertex and from the weak
vertex, respectively. Noteworthy, the expression in large square brackets
behaves at small $q^{2}$ as $\propto q^{2}$, because of the constraint
$F_{2}(q^{\prime 2},q^{2}=0)=0$ 222One may write for $F_{2}(q^{\prime
2},q^{2})$ dispersion representation in $q^{2}$ with one subtraction leading
to $F_{2}(q^{\prime 2},q^{2})=O(q^{2})$ at $q^{2}\to 0$. Neglecting all lepton
masses in the amplitude, but keeping them in the phase-space volume, the
double differential distribution then takes the form (we display separately
all numerical factors according to the definition of the differential
distribution):

$\displaystyle\frac{d^{2}\Gamma(B\to
lll^{\prime}\nu^{\prime})}{dq^{2}dq^{\prime
2}}=\frac{(2\pi)^{4}}{2M_{B}}\frac{1}{(2\pi)^{12}}\frac{\pi\lambda^{1/2}(q^{2},m_{l}^{2},m_{l}^{2})}{2q^{2}}\frac{\pi\lambda^{1/2}(q^{\prime
2},m_{l^{\prime}}^{2},0)}{2q^{\prime
2}}\frac{\pi\lambda^{1/2}(M_{B}^{2},q^{2},q^{\prime 2})}{2M_{B}^{2}}|A|^{2},$
(4.25)

where $\lambda(a,b,c)=(a-b-c)^{2}-4bc$ is the triangle function. The
kinematical constrains on the variables $q^{2}$ and $q^{\prime 2}$ come from
the $\lambda$-functions in Eq. (4.24) and read

$\displaystyle 4m_{l}^{2}\leq q^{2},\quad m_{l^{\prime}}^{2}\leq q^{\prime
2},\quad\sqrt{q^{2}}+\sqrt{q^{\prime 2}}\leq M_{B}.$ (4.26)

The first remark we would like to make is the following: because of the gauge-
invariance constrains on the form factors, $F_{2}(q^{\prime 2},q^{2}=0)=0$,
and therefore $|A|^{2}\propto 1/q^{2}$, and not $\propto 1/(q^{2})^{2}$ as may
be naivley obtained when the gauge-invariance constrains are not taken into
account. As the result, the decay rate $\Gamma(B\to lll^{\prime}\nu^{\prime})$
has a mild logarithmic dependence $\log(m_{l}^{2})$ on the leptom mass
$m_{l}$.333 In nikitin , the relation $\Gamma(B\to
e^{+}e^{-}\mu\nu_{\mu})/\Gamma(B\to\mu^{+}\mu^{-}e\nu_{e})\propto(m_{\mu}/m_{e})^{2}$
was reported. This relation is however misleading as it arises when the
constraints coming from gauge invariance of the $B\to
lll^{\prime}\nu^{\prime}$ amplitude are not taken into account. As our
calculation in the next Section demonstrates, a substantial contribution to
the decay rate $\Gamma(B\to lll^{\prime}\nu^{\prime})$ comes from the
$q^{2}$-region of light vector resonances which is not sensitive to the lepton
masses; the region of soft $q^{2}\simeq 4m_{l}^{2}$ gives not more than a half
of the full rate even for $m_{l}=m_{e}$. In the end, $\Gamma(B\to
e^{+}e^{-}\mu\nu_{\mu})$ turns out to be only a factor $\sim 2$ larger than
$\Gamma(B\to\mu^{+}\mu^{-}e\nu_{e})$. The relation $\Gamma(B\to
e^{+}e^{-}\mu\nu_{\mu})/\Gamma(B\to\mu^{+}\mu^{-}e\nu_{e})\propto(m_{\mu}/m_{e})^{2}$
proposed in nikitin  is not justified and contradicts to the constraints
coming from gauge invariance of the $B\to lll^{\prime}\nu^{\prime}$ amplitude.

It should be emphasized that the double differential distribution
$d^{2}\Gamma(B\to l^{+}l^{-}l^{\prime}\nu^{\prime})/dq^{2}dq^{\prime 2}$ is
easily calculable due to the fact that the leptons emitted from the
electromagnetic vertex and the lepton emitted from the weak vertex have
different flavours; no exchange diagrams emerge in this case and one obtains
the explicit analytic expression for the double differential distribution.

##  5 The form factors and the differential distributions

###  5.1 Modelling the form factors

The form factors $F_{A1,A2}(q^{\prime 2},q^{2})$ are obtained from the form
factors $f_{1A,2A,3A}$, using for the latter the $q^{2}$-dispersion
representation with one subtraction at $q^{2}=0$, Eqs. (2.11), (2.12), and
(2.13). Similarly, for the for factor $F_{V}(q^{\prime 2},q^{2})$ a single-
subtracted dispersion representation in $q^{2}$ is used. The subtraction
procedure allows us to incorporate the constraints imposed by gauge invariance
and the known properties of the form factors at $q^{2}=0$ (i.e., corresponding
to the real photon in the final state).

Furthermore, we assume that the spectral densities are saturated by the light
vector-meson resonances $\rho^{0}$ and $\omega$ in the $q^{2}$-channel and,
since these resonances emerge in the physical region of the $B$-decay of
interest, we take into account the finite-width effects of these resonances.
In the end, we come to the following expressions

$\displaystyle F_{A1}(q^{\prime 2},q^{2})$ $\displaystyle=$ $\displaystyle
F_{A}(q^{\prime
2})-\frac{Q_{B}f_{B}}{q^{\prime}q}-q^{2}\sum\limits_{V=\rho^{0},\omega}\bigg{(}\frac{1}{M_{V}^{2}}\frac{2M_{B}}{M_{B}^{2}-M_{V}^{2}-q^{\prime
2}}\frac{M_{V}f_{V}}{M_{V}^{2}-q^{2}-i\Gamma_{V}(q^{2})M_{V}}(M_{B}+M_{V})A_{1}^{B\to
V}(q^{\prime 2})\bigg{)},$ $\displaystyle F_{A2}(q^{\prime 2},q^{2})$
$\displaystyle=$
$\displaystyle-q^{2}\sum\limits_{V=\rho^{0},\omega}\frac{1}{M_{V}^{2}}\frac{M_{V}f_{V}}{M_{V}^{2}-q^{2}-i\Gamma_{V}(q^{2})M_{V}}\left[\frac{2M_{B}}{M_{B}^{2}-M_{V}^{2}-q^{\prime
2}}(M_{B}+M_{V})A_{1}^{B\to V}(q^{\prime 2})-\frac{A_{1}^{B\to V}(q^{\prime
2})}{(M_{B}+M_{V})}\right]$ (5.28)
$\displaystyle+Q_{B}f_{B}\left(\frac{2M_{B}}{M_{B}^{2}-q^{\prime
2}}-\frac{2M_{B}}{M_{B}^{2}-q^{\prime 2}-q^{2}}\right),$ $\displaystyle
F_{V}(q^{\prime 2},q^{2})$ $\displaystyle=$ $\displaystyle F_{V}(q^{\prime
2})-q^{2}\sum\limits_{V=\rho^{0},\omega}\bigg{(}\frac{1}{M_{V}^{2}}\frac{M_{V}f_{V}}{M_{V}^{2}-q^{2}-i\Gamma_{V}(q^{2})M_{V}}\frac{2V^{B\to
V}(q^{\prime 2})}{M_{B}+M_{V}}\bigg{)},$ (5.29)

The expressions above contain several ingredients. Let us briefly discuss
them:

$\bullet$ The form factors $F_{A}(q^{\prime 2})$ and $F_{V}(q^{\prime 2})$
describe the $B\to\gamma l^{\prime}\nu^{\prime}$ transition; they emerge as
subtraction terms at $q^{2}=0$ in the $q^{2}$-disperison representations for
the form factors $F_{1A,V}(q^{\prime 2},q^{2})$ at $q^{2}=0$. The form factors
$F_{A}(q^{\prime 2})$ and $F_{V}(q^{\prime 2})$ are equal to each other at the
leading order of the double $1/E_{\gamma}$
($2M_{B}E_{\gamma}=M_{B}^{2}-q^{\prime 2}$) and $1/M_{B}$ expansions in QCD
korchemsky  but differ in the subleading orders mn2004 ; beneke1 :

$\displaystyle F_{A,V}(q^{\prime 2})$ $\displaystyle=$
$\displaystyle-\frac{Q_{u}M_{B}f_{B}}{2E_{\gamma}\lambda_{B}}+O(Q_{b}f_{B}/E\gamma)+O(Q_{u}f_{B}M_{B}/E_{\gamma}^{2}).$
(5.30)

The magnitude of the form factors $F_{A}(q^{\prime 2})$ and $F_{V}(q^{\prime
2})$ is determined to a large extent by the parameter $\lambda_{B}$, the
inverse moment of the $B$-meson light-cone distribution amplitude $\phi_{B}$
korchemsky . The value of $\lambda_{B}$ presently has a large uncertainty: for
instance, beneke1  makes use of $\lambda_{B}(1\,{\rm GeV})=0.35$ GeV; the sum-
rule estimate of kou  led to $\lambda_{B}(1\,{\rm GeV})=0.57$ GeV; the results
of calculating $F_{A}(q^{2})$ and $F_{V}(q^{\prime 2})$ mn2004  using the
dispersion approach m  correspond to a relatively large value
$\lambda_{B}(1\,{\rm GeV})=0.657$ GeV. Obviously, the uncertainty in the
parameter $\lambda_{B}$ dominates the uncertainty in the differential
distributions $d\Gamma(B\to lll^{\prime}\nu^{\prime})/dq^{2}$ at small
$q^{2}$. Still, the most important qualitative feature of the distributions at
low $q^{2}$, namely its singular behaviour as $d\Gamma(B\to
lll^{\prime}\nu^{\prime})/dq^{2}\propto 1/q^{2}$ is fully determined by the
constraints coming from the gauge-invariance of the amplitude (2.6).

$\bullet$ In the region $0.4\leq q^{2}({\rm GeV}^{2})\leq 0.9$, where light
vector meson resonances show up in the differential distributions, the form
factors of interest cannot be calculated using perturbative QCD ims . To
calculate the form factors $F_{1A,2A,V}(q^{\prime 2},q^{2})$ in this region of
$q^{2}$ and for any $q^{\prime 2}$ appropriate for the four-lepton decay, we
make use of the dispersion representations and assume kmn2018 ; nikitin  that
they may be saturated by the intermediate $\rho^{0}$ and $\omega$-states in
the $q^{2}$-channel.444We would like to notice that no relative phase between
the $\rho^{0}$ and $\omega$ contributions to the form factors $F_{i}(q^{\prime
2},q^{2})$ as proposed in nikitin  may emerge: these form factors contain a
sum over the intermediate states $|V\rangle\langle V|$, so even if one
introduces arbitrary complex phases in the states $|V\rangle$, these phases
appear both in the decay constants $f_{V}$ and the $B\to V$ weak form factors
such that they finally drop out from $F_{i}(q^{\prime 2},q^{2})$. Since the
light neutral vector mesons lie in the physical decay region of $q^{2}$, it is
necessary to take into account their finite $q^{2}$-dependent width
$\Gamma_{V}(q^{2})$. For a relatively broad $\rho$-meson the function
$\Gamma_{V}(q^{2})$ takes into account the effects of the $2\pi$ intermediate
states; the appropriate formulas are given in nachtmann . In practical
calculations, we use a simplified expression which takes into account the
correct threshold behaviour of the $\rho\to\pi\pi$ phase space:
$\Gamma_{\rho^{0}}(q^{2})=\theta(q^{2}-4m_{\pi}^{2})(1-4m_{\pi}^{2}/q^{2})^{3/2}/(1-4m_{\pi}^{2}/M_{\rho}^{2})^{3/2}\Gamma_{\rho^{0}}$.
For a narrow $\omega$-meson, we take an approximation of constant width (This
approximation is not fully theoretically clean: the imaginary part of the
propagator of the vector meson should vanish below the threshold in the
corresponding decay channel. This means that $\Gamma_{V}(q^{2})$ should vanish
below the corresponding light-meson threshold. But for a narrow $\omega$-meson
the effect is tiny). Table 1 gives the meson parameters entering the form
factors Eqs. (5.1), (5.28), and (5.29).

We notice that making use of the contribution of the continuum of light
pseudoscalar mesons to $\Gamma_{V}(q^{2})$, we effectively take into account
the contribution of this hadron continuum in the spectral densities of the
form factors $F_{1A,2A,V}(q^{\prime 2},q^{2})$ nachtmann . In the end, one
finds that the nonresonance $q^{2}$-region, $q^{2}\geq$1.0 GeV2 gives a small
contribution to the decay width of the $B\to lll^{\prime}\nu^{\prime}$ decay.
This agrees with the expectations of the analysis of beneke2 .

Table 1:  Meson parameters entering the expressions for the form factors, Eqs. (5.1), (5.28), and (5.29). Data from ball2007 ; pdg . $f_{B}$ | $\sqrt{2}f_{\rho}^{0}$ | $\sqrt{2}f_{\omega}$ | $\Gamma_{\rho^{0}}$ | $\Gamma_{\omega}$
---|---|---|---|---
190 MeV | 216 MeV | 190 MeV | 150 MeV | 8.49 MeV

$\bullet$ The contributions of the light vector mesons are expressed via the
form factors $A_{1}^{B\to V}(q^{\prime 2})$, $A_{2}^{B\to V}(q^{\prime 2})$,
and $V^{B\to V}(q^{\prime 2})$ describing the weak decay $B\to V$. In spite of
many efforts to calculate these form factors in a broad kinematical decay
region $0<q^{\prime 2}<M_{B}^{2}$, our knowledge of these quantities is not
very accurate. Table 2 presents some selected results for the relevant form
factors: although the central values of the form factors at $q^{\prime 2}=0$
from different approaches are in reasonable agreement with each other, the
uncertainties vary from an “educated guess” of 10% for ms2000  to almost 50%
in gubernari2019 . The uncertainties in these form factors, along with the
uncertainty in the parameter $\lambda_{B}$, is the second main source of the
uncertainty in the theoretical predictions for $B\to
l^{+}l^{-}l^{\prime}\nu^{\prime}$ decays.

Table 2:  Selected theoretical predictions for the weak form factors describing $B$ decays into light vector mesons. The form factors from ms2000 ; ballzwicky2005 ; ivanov2  are expected to have a 10-15% uncertainty. To obtain the form factors for $B^{-}\to\rho^{0}$ and $B^{-}\to\omega$ decays, the numbers given in this Table should be multiplied by the isotopic factor $1/\sqrt{2}$. Ref. | $A_{1}^{B\to\rho}(0)$ | $A_{1}^{B\to\omega}(0)$ | $A_{2}^{B\to\rho}(0)$ | $A_{2}^{B\to\omega}(0)$ | $V^{B\to\rho}(0)$ | $V^{B\to\omega}(0)$
---|---|---|---|---|---|---
ms2000  | 0.26 | $-$ | 0.24 | $-$ | 0.31 | $-$
ballzwicky2005  | 0.24 | 0.22 | 0.22 | 0.20 | 0.32 | 0.29
ivanov2  | 0.26 | $-$ | 0.24 | $-$ | 0.28 | $-$
gubernari2019  | $0.22\pm 0.1$ | $-$ | $0.19\pm 0.11$ | $-$ | $0.27\pm 0.14$ | $-$

###  5.2 The differentail distributions

With the analytic expressions for the form factors, Eqs. (5.1), (5.28), and
(5.29) at hand, Eqs. (4.25) give the differential distributions in
$B\to\mu^{+}\mu^{-}e\nu_{e}$ decays, Fig. 2. Here we use $V_{ub}=0.004$ and
$\tau_{B^{-}}=1.638\,10^{-12}$ s. The Plots show the impact of the parameter
$\lambda_{B}$ on the differential distributions $d\Gamma/dq^{2}$ and
$d\Gamma/dq^{\prime 2}$.

|
---|---
(a) | (b)
Fig. 2:  The differential distrbutions
$d\Gamma(B\to\mu^{+}\mu^{-}e\nu_{e})/dq^{2}$ (a) and
$d\Gamma(B\to\mu^{+}\mu^{-}e\nu_{e})/dq^{\prime 2}$ (b), for the weak
transition $B\to V$ form factors from ms2000 . Solid lines corresponds to
$\lambda_{B}=0.35$, dashed lines correspond to $\lambda_{B}=0.65$.

Figure 3 shows the double differential distributions calculated for
$\lambda_{B}=0.65$ GeV and the form factors from gubernari2019 .

Fig. 3:  The double differential distribution
$d^{2}\Gamma(B\to\mu^{+}\mu^{-}e\nu_{e})/dq^{2}dq^{\prime 2}$.

##  6 Conclusions

We presented a theoretical analysis of $B\to l^{+}l^{-}l^{\prime}\nu^{\prime}$
decays with special emphasis on the constraints imposed by gauge invariance.
Our results are summarized below:

1\. Gauge invariance provides essential constraints on the amplitude of Eq.
(2.6):

$\displaystyle T^{5}_{\alpha,\nu}$ $\displaystyle=$ $\displaystyle i\int
dxe^{iqx}\langle 0|T\left\\{j^{\rm
e.m.}_{\alpha}(x),\bar{u}\gamma_{\nu}\gamma_{5}b(0)\right\\}|\bar{B}_{u}(p)\rangle.$
(6.31)

We emphasize that for a consistent analysis of the $B\to
l^{+}l^{-}l^{\prime}\nu^{\prime}$ amplitude it is necessary to start with the
amplitude (6.31) and properly parameterize this amplitude taking into account
all constraints imposed by the electromagnetic gauge invariance and
analyticity. Taking into account these constraints leads to the relation

$\displaystyle F_{2A}(q^{\prime 2},q^{2}=0)=0.$ (6.32)

The latter relation is crucial as the form factor $F_{2A}(q^{\prime 2},q^{2})$
determines the behaviour of the decay rate $d^{2}\Gamma(B\to
lll^{\prime}\nu^{\prime})/dq^{2}dq^{\prime 2}$ at small $q^{2}$ (see Eq.
(4.24)). Making use of the constraint (6.32), one finds the following
small-$q^{2}$ behaviour of the differential distribution:

$\displaystyle d\Gamma(B\to lll^{\prime}\nu^{\prime})/dq^{2}\propto 1/q^{2}$
(6.33)

and, consequently, a mild logarithmic dependence of the decay rate on the
leptom mass $m_{l}$

$\displaystyle\Gamma(B\to lll^{\prime}\nu^{\prime})\propto\log(m_{l}^{2}).$
(6.34)

If the constraint (6.32) is not taken into account, one comes to a wrong
behaviour $d\Gamma(B\to lll^{\prime}\nu^{\prime})/dq^{2}\propto 1/q^{4}$ and,
consequently, $\Gamma(B\to lll^{\prime}\nu^{\prime})\propto 1/m_{l}^{2}$.

2\. We derived an explicit analytic expression for the differential
distributions $d\Gamma/dq^{2}dq^{\prime 2}$ in $B\to lll^{\prime}\nu$ decays.
For the form factors $F_{V,1A,2A}(q^{\prime 2},q^{2})$, describing the
amplitude of $B\to lll^{\prime}\nu^{\prime}$ decay, we obtained dispersion
representations in $q^{2}$ with one subtraction; this allows us to take into
account properly all constraints imposed by gauge invariance. Assuming that
the spectral densities are saturated by the contributions of light vector
mesons $\rho^{0}$ and $\omega$, we obtained analytic representations for the
form factors in a broad range of $q^{2}$ and $q^{\prime 2}$; these
representations involve the weak form factors describing $B\to\rho$ and
$B\to\omega$ decays.

3\. We obtained numerical predictions for the differential and integrated
branching ratios of the $B\to l^{+}l^{-}l^{\prime}\nu^{\prime}$ decays. The
uncertainties in these predictions come from the two main sources: (i) The
uncertainty in the parameter $\lambda_{B}$, which mainly governs the behaviour
of the differential distributions at small $q^{2}\leq 0.4$ GeV2. We allow the
parameter $\lambda_{B}$ to vary in the range $\lambda_{B}=0.35-0.65$ GeV (the
lower values of this range has been advocated in several analyses korchemsky ;
beneke1  whereas the upper values of $\lambda_{B}$ is obtained in explicit
model calculations kou ; mn2004 ). (ii) The uncertainties in the
$B\to\omega,\rho$ weak form factors $V,A_{1},A_{2}$, which mainly govern the
differential distributions in the region of $q^{2}=0.4-0.9$ GeV2. In the
numerical estimates we make use the recent results from gubernari2019  and
took into account the suppression of the $B\to\omega$ form factors compared to
the corresponding $B\to\rho$ form factors according to ballzwicky2005 . Taking
into account the uncertainties mentioned above, we obtain the following
estimates

$\displaystyle{\rm Br}(B\to\mu^{+}\mu^{-}e\nu_{e})=(4.0\pm
0.5_{\lambda_{b}}\pm 2.3_{\rm weak\,ffs})\,10^{-8},$ $\displaystyle{\rm
Br}(B\to e^{+}e^{-}\mu\nu_{\mu})=(7.0\pm 2.0_{\lambda_{b}}\pm 2.7_{\rm
weak\,ffs})\,10^{-8}.$ (6.35)

Both decay rates have the same order of magnitude and ${\rm Br}(B\to
e^{+}e^{-}\mu\nu_{\mu})$ is expected to be only a factor 2 larger than ${\rm
Br}(B\to\mu^{+}\mu^{-}e\nu_{e})$.

4\. For the decay $B\to l^{+}l^{-}l^{+}\nu_{l}$ ($l=\mu,e$) with identical
positive-charged leptons in the funal state, the amplitude is given by the sum
of direct and exchange diagrams, $A=A_{\rm dir}+A_{\rm exchange}$, and the
phase space includes a factor 1/2 because of the presence of the identical
particles in the final state. The phase-space integration of both $|A_{\rm
dir}|^{2}$ and $|A_{\rm exchange}|^{2}$ leads to the same result, $\Gamma(B\to
lll^{\prime}\nu^{\prime})|_{m_{l^{\prime}}=m_{l}}$, and one can write (see,
e.g., beneke2 ):

$\displaystyle\Gamma(\to lll\nu)=\Gamma(B\to
lll^{\prime}\nu^{\prime})|_{m_{l^{\prime}}=m_{l}}+\Gamma_{\rm
interference}(B\to lll\nu),$ (6.36)

where $\Gamma_{\rm interference}(B\to lll\nu)$ is the phase-space integral of
$(A_{\rm dir}A^{*}_{\rm exchange}+A^{*}_{\rm dir}A_{\rm exchange})$. The
inteference term should be calculated numerically as the integral over the
phase space. A simple analytic result similar to Eqs. (4.24) and (4.25) cannot
be obtained. We did not perform explicit calculations of the interference
branching fraction. However, making use of the estimates for the exchange-
diagram contributions from nikitin ; beneke2 , the interference branching
fraction leads to a reduction of ${\rm Br}(B\to lll^{\prime}\nu^{\prime})$ by
10-15%. So we expect that

$\displaystyle{\rm Br}(B^{+}\to\mu^{+}\mu^{-}\mu^{+}\bar{\nu}_{\mu})=(3.5\pm
0.4_{\lambda_{b}}\pm 2.0_{\rm weak\,ffs})10^{-8}.$ (6.37)

This value agrees with a semi-qualitative estimates from beneke2  and is only
marginally compatible with the upper limits obtained by the LHCb Collaboration
exp1  ${\rm Br}(B^{+}\to\mu^{+}\mu^{-}\mu^{+}\bar{\nu}_{\mu})\leq
1.6\dot{\,}10^{-8}$.

###### Acknowledgements.

D. M. gratefully acknowledges support from the Russian Foundation of Basic
Research under joint RFBR-CNRS grant 19-52-15022.

## References

  * (1) LHCb Collaboration (R. Aaij et al.), Search for the rare decay $B_{s}^{0}\to\mu^{+}\mu^{-}\mu^{+}\mu^{-}$, Phys. Rev. Lett. 110, 211801 (2013).
  * (2) ATLAS Collaboration (M. Aaboud et al.), Study of the rare decays of $B^{0}_{s}$ and $B_{0}$ into muon pairs from data collected during the LHC Run 1 with the ATLAS detector, Eur. Phys. J. C76, 513 (2016).
  * (3) LHCb Collaboration (R. Aaij et al.), Search for decays of neutral beauty mesons into four muons, JHEP 1703, 001 (2017).
  * (4) LHCb Collaboration (R. Aaij et al.), Search for the rare decay $B^{+}\to\mu^{+}\mu^{-}\mu^{+}\nu_{\mu}$, Eur. Phys. J. C79, 675 (2019).
  * (5) Y. Dincer and L. M. Sehgal, Electroweak effects in the double Dalitz decay $B(s)\to l^{+}l^{-}l^{\prime+}l^{\prime-}$, Phys. Lett. B 556, 169 (2003).
  * (6) A. V. Danilina and N. V. Nikitin, Four-Leptonic Decays of Charged and Neutral $B$ Mesons within the Standard Model, Phys. Atom. Nucl. 81, 347 (2018), Yad. Fiz. 81, 331 (2018); A. Danilina N. Nikitin, and K. Toms, Decays of charged $B$-mesons into three charged leptons and a neutrino, Phys. Rev. D101, 096007 (2020).
  * (7) M. Beneke, P. Böer, P. Rigatos, K. K. Vos, QCD factorization of the four-lepton decay $B\to lll\nu$, arXiv:2102.10060.
  * (8) M. A. Ivanov, D. Melikhov, S. Simula, Form factors for $B\to j_{1}j_{2}$ decays into two currents in QCD, Phys. Rev. D 101, 094022 (2020).
  * (9) M. Beyer, D. Melikhov, N. Nikitin, B. Stech, Weak annihilation in the rare radiative $B\to\rho\gamma$ decay, Phys. Rev. D64, 094006 (2001).
  * (10) D. Melikhov, Dispersion approach to quark binding effects in weak decays of heavy mesons, Eur. Phys. Journal direct 4, 2 (2002) [hep-ph/0110087].
  * (11) F. Kruger and D. Melikhov, Gauge invariance and form-factors for the decay $B\to\gamma l^{+}l^{-}$, Phys. Rev. D67, 034002 (2003).
  * (12) A. Kozachuk, D. Melikhov, N. Nikitin, Annihilation type rare radiative $B_{(s)}\to V\gamma$ decays, Phys. Rev. D93, 014015 (2016).
  * (13) T. M. Aliev, A. Ozpineci, and M. Savci, $B_{q}\to l^{+}l^{-}\gamma$ decays in light cone QCD, Phys. Rev. D55, 7059 (1997).
  * (14) G. Korchemsky, D. Pirjol, T.-M. Yan, Radiative leptonic decays of $B$ mesons in QCD, Phys. Rev. D61, 114510 (2000).
  * (15) P. Ball and E. Kou, $B\to\gamma e\nu$ transitions from QCD sum rules on the light cone, JHEP 0304, 029 (2003).
  * (16) D. Melikhov and N. Nikitin, Rare radiative leptonic decays $B_{(d,s)}\to\gamma l^{+}l^{-}$, Phys. Rev. D70, 114028 (2004).
  * (17) M. Beneke, J. Rohrwild, B meson distribution amplitude from $B\to\gamma l\nu$, Eur. Phys. J.bf C71, 1818 (2011).
  * (18) A. Kozachuk, D. Melikhov, and N. Nikitin, Rare FCNC radiative leptonic $B_{s,d}\to\gamma l^{+}l^{-}$ decays in the Standard Model, Phys. Rev. D97, 053007 (2018).
  * (19) S. Dubnicka, A. Z. Dubnickova, M. A. Ivanov, A. Liptaj, P. Santorelli, and C. T. Tran, Study of $B_{s}\to l^{+}l^{-}\gamma$ decays in covariant quark model, Phys. Rev. D99, 014042 (2019).
  * (20) S. Treiman, R. Jackiw, D. Gross, Lectures on current algebra and its applications, Princeton University Press, Princeton, New Jersey, 1972.
  * (21) D. Melikhov, Form factors of meson decays in the relativistic constituent quark model, Phys. Rev. D53, 2460 (1996); Heavy quark expansion and universal form-factors in the quark model Phys. Rev. D56, 7089 (1997).
  * (22) D. Melikhov, O. Nachtmann, V. Nikonov, T. Paulus, Masses and couplings of vector mesons from the pion electromagnetic, weak, and $\pi\gamma$ transition form-factors, Eur. Phys. J. C 34, 345 (2004).
  * (23) P. Ball, G. W. Jones and R. Zwicky, $B\to V\gamma$ beyond QCD factorisation, Phys. Rev. D 75, 054004 (2007).
  * (24) Particle Data Group, C. Patrignani et al., Chin. Phys. C, 40, 100001 (2016).
  * (25) D. Melikhov, B. Stech, Weak form-factors for heavy meson decays: an update, Phys. Rev. D62, 014006 (2000).
  * (26) P. Ball and R. Zwicky, $B_{d,s}\to\rho,\omega,K^{*},\phi$ decay form factors from light-cone sum rules reexamined, Phys. Rev. D71, 014029 (2005).
  * (27) M. A. Ivanov, J. G. Korner, S. G. Kovalenko, P. Santorelli, G. G. Saidullaeva, Form factors for semileptonic, nonleptonic and rare $B(B_{s})$ meson decays, Phys. Rev. D85, 034004 (2012).
  * (28) N. Gubernari, A. Kokulu and D. van Dyk, $B\to P$ and $B\to V$ form factors from $B$-meson light-cone sum rules beyond leading twist, JHEP 01, 150 (2019).