QCD factorization of the four-lepton decay $B^-\rightarrow \ell \bar{\nu}_\ell \ell^{(\prime)} \bar{\ell}^{(\prime)}$

Motivated by the first search for the rare charged-current $B$ decay to four leptons, $\ell \bar{\nu}_\ell \ell^{(\prime)} \bar{\ell}^{(\prime)}$, we calculate the decay amplitude with factorization methods. We obtain the $B\to \gamma^*$ form factors, which depend on the invariant masses of the two lepton pairs, at leading power in an expansion in $\Lambda_{\rm QCD}/m_b$ to next-to-leading order in $\alpha_s$, and at $\mathcal{O}(\alpha_s^0)$ at next-to-leading power. Our calculations predict branching fractions of a few times $10^{-8}$ in the $\ell^{(\prime)} \bar{\ell}^{(\prime)}$ mass-squared bin up to $q^2=1~\text{GeV}^2$ with $n_+q>3~$GeV. The branching fraction rapidly drops with increasing $q^2$. An important further motivation for this investigation has been to explore the sensitivity of the decay rate to the inverse moment $\lambda_B$ of the leading-twist $B$ meson light-cone distribution amplitude. We find that in the small-$q^2$ bin, the sensitivity to $\lambda_B$ is almost comparable to $B^- \rightarrow \ell^- \bar{\nu}_\ell\gamma$ when $\lambda_B$ is small, but with an added uncertainty from the light-meson intermediate resonance contribution. The sensitivity degrades with larger $q^2$.


Introduction
The radiative decay B − → −ν γ has been extensively studied in the context of QCD factorization (QCDF) [1][2][3][4][5] when the energy of the photon E γ is large compared to the scale of the strong interaction Λ QCD . At leading power in a simultaneous expansion in Λ QCD /E γ and Λ QCD /m b , and at leading order (LO) in the strong coupling α s , the relevant B → γ transition form factor can be expressed in terms of only two hadronic parameters: the accurately known B meson decay constant f B , and the poorly constrained first inverse moment 1/λ B = ∞ 0 dω φ B + (ω)/ω of φ B + (ω), the leading-twist B meson light-cone distribution amplitude (LCDA). This hadronic parameter was introduced in the theoretical description of charmless hadronic decays [6] and appears in the QCD calculation of almost any other exclusive B decay to light particles. The radiative decay B − → −ν γ has been advocated as a means to determine λ B from data [5]. First significant measurements can be expected from the BELLE II experiment (see [7] for the most recent BELLE result).
This strategy is difficult to implement in the hadronic B experiment LHCb, since the photon in the radiative decay cannot be easily reconstructed. In this paper we investigate whether the four-lepton decay B − → ν γ * → ν ( )¯ ( ) , in which the real photon is replaced by a virtual one, which decays into a lepton-antilepton pair ( , = e, µ), retains sensitivity to λ B , and hence could provide an alternative measurement. We focus on the kinematic region, where the γ * , respectively the lepton pair, has large energy but small invariant mass q 2 6 GeV 2 . 1 The four-lepton decays have not been observed up to now, but the LHCb experiment [8] established an upper bound of Br (B + → µ +ν µ µ − µ + ) < 1.6 · 10 −8 on the branching fraction of the muonic mode under the assumption that the smaller of the two possible µ + µ − invariant masses is below 980 MeV, which is close to, in fact somewhat below, theoretical expectations [9,10].
The factorization theorem for the B → γ form factors in the regime where the photon is energetic, n + q = 2E γ Λ QCD , has been established long ago [3,4]. Its generalization to B → γ * form factors is straightforward, when q 2 is away from light-meson resonances. The present treatment follows the strategy applied to B − → −ν γ [5] and B s → µ + µ − γ [11]we compute the form factor in QCD factorization at leading power (LP) including O(α s ) QCD corrections, and include next-to-leading power (NLP) corrections at O(α 0 s ). The light-meson resonance contribution is included in the same fashion as for the "type-B" contribution to B s → µ + µ − γ [11]. Since the four-lepton final state is produced from a virtual W boson and photon, an extension of previous calculations is required to B → γ * form factors that depend on two non-vanishing virtualities. We note that a previous computation [12] of these B → γ * form factors includes either only the resonance contribution at small q 2 , or employs QCD sum rules that apply only to large q 2 ∼ m 2 b . No attempts have so far been undertaken to estimate the form factors for intermediate and small q 2 with factorization methods, as done here, which apply when n + q Λ QCD . With these kinematic restrictions the differential branching fraction of the four-lepton decay is expressed, at LP, in terms of generalized inverse moments of the B meson LCDA, which can be related to λ B .
We consider the case of non-identical lepton flavours, = , and identical ones, which requires additional kinematic considerations.

Basic definitions
Following the conventions of [5], we write the B − → ν ¯ decay amplitude to lowest non-vanishing order in the electromagnetic coupling as where q ≡ q 1 +q 2 and p = m B v = q+k, such that k = p +p ν is the momentum of the virtual W boson. In addition, we use the convention iD µ = i∂ µ − Q eA µ em for the electromagnetic covariant derivative, with Q = −1 for the lepton fields. The hadronic tensor with the electromagnetic current j µ em = q Q qq γ µ q + Q ¯ γ µ accounts for the emission of the virtual photon from the B meson constituents. The second term in the square brackets in (2.1) corresponds to the emission from the final-state lepton, see Fig. 2 below. It can be expressed in terms of the B meson decay constant 0 ūγ ν (1 − γ 5 ) b B − (p) = −if B p ν and constitutes a power correction relative to the T µν term in the kinematic region of interest.
The hadronic tensor T µν can be decomposed into six form factors F i (q 2 , k 2 ) of two kinematic invariants. Applying the Ward identity q µ T µν = f B p ν leaves four form factors and a contact term (see App. A for details). We write We neglect the lepton masses, in which case the q µ , k ν terms drop out after contracting T µν (p, q) with the lepton tensor. The contact term is fixed by the Ward identity to This can be rewritten as f B /(v · q) v µ (k ν + q ν ) and has been absorbed intoF A and the k ν terms in (2.3). The convention for the totally anti-symmetric tensor is 0123 = 1. The virtual photon emission from the final-state lepton in (2.1) is exactly cancelled by the redefinition Therefore, the term in square brackets in (2.1) can be expressed in terms of three form factors. To separate amplitudes corresponding to the different polarization states of the virtual photon, we shall use the decomposition which implies The form factor F A arises from a longitudinally polarized virtual photon and vanishes in the realphoton limit q 2 → 0. Without loss of generality we choose the three-momentum q to point in the positive z direction, such that its decomposition into light-cone vectors n µ ± reads with n µ ± = (1, 0, 0, ∓1) and q 2 = n + q n − q. The transverse metric tensor is then g µν ⊥ = g µν − (n µ + n ν − + n ν + n µ − )/2. The large component n + q of q µ is related to the invariant masses q 2 and k 2 via Finally, we define the left-and right-helicity form factors Since helicity is conserved in high-energy QCD processes, F R is power-suppressed relative to F L in the heavy-quark / large n + q limit. For non-identical lepton flavours = , the differential decay width can be obtained analytically by a straightforward calculation. The full angular distribution, i.e. the fivefold differential rate is given in App. B. Here we quote the double differential rate in the invariant masses q 2 and k 2 , for which we obtain the simple expression keeping the lepton masses m ( ) in the phase space integration (implying q 2 > 4m 2 ), which is relevant for muons. The case of identical lepton flavours = is more complicated, as an additional contribution from the interchange of the two final-state leptons arises. To clarify this point, let . For the decay rate, this results in an additional interference term between M a and M b , while the rates Γ a,b ∝ |M a,b | 2 from the squares of the individual diagrams are equal (as depicted in Fig. 1). Since Γ a + Γ b is equal to the rate for non-identical lepton flavours, we find (2.11) For the interference term d 2 Br int B − → ν ¯ /(dq 2 dk 2 ) can only be obtained numerically.

"Squared Diagrams"
Rigatos Panagiotis Figure 1: Graphical representation of the squared amplitude for the non-identical-lepton final-state (left) and of the interference term for the case of identical lepton flavours (right).

Calculation of the form factors
The amplitude can be factorized through an expansion in Λ QCD /m b and Λ QCD /n + q, if the quark propagator that connects the electromagnetic and the weak current is far offshell. This happens for very large q 2 of order m 2 B in which case the amplitude can be reduced to a hard matching coefficient times the B meson decay constant defined as a local matrix element in heavy-quark effective theory (HQET). The decay rate for such large q 2 is highly suppressed. The situation is more interesting when q 2 m 2 B , but q µ has still a large component n + q ∼ O(m B ), while n − q ∼ O(Λ QCD ), or even smaller. In this case the intermediate quark propagator becomes hard-collinear and the γ * probes the light-cone structure of the B meson. A factorization formula, which expresses the form factors as a convolution of the B meson LCDA with a perturbative scattering kernel, can then be derived for the LP contribution using soft-collinear effective theory [13][14][15][16] by matching QCD → SCET I → HQET. Since the derivation is very similar to the one for B − → −ν γ and B s → µ + µ − γ decays [5,11], we only sketch the main steps in the following.
Upon integrating out the hard scales m b , n + q, the flavour-changing weak current is represented in SCET I bȳ with hard matching coefficients C i = C i (n + q; µ). Here q hc = W † ξ hc is the hard-collinear quark field in SCET, multiplied with a hard-collinear Wilson line to ensure SCET collinear gauge-invariance. Fields without arguments live at space-time point x = 0. At LP, the index µ is transverse, since the LP SCET I electromagnetic current j µ q,SCET I (x) [3] contains only the transverse polarization of the virtual photon. We therefore only need C with α s ≡ α s (µ) in the MS scheme, and z = n + q/m B = 1 − k 2 /m 2 B + O(Λ QCD /m B ). The hadronic tensor is then expressed as in terms of the matching coefficient and the SCET I correlation function A discussion of the precise power counting of the individual terms in j µ q,SCET I and the possibility of extrapolating the above expressions to the large q 2 region with tree-level accuracy can be found in [11].
The SCET I correlation function is then matched at LP to HQET. This results in where n − q = q 2 /n + q. The hard-collinear matching function [17] J(n + q, q 2 , ω; µ) is convoluted with the leading-twist B meson LCDA φ B + (ω) defined through which contains the scale-dependent HQET B meson decay constant F B = F B (µ). Similar to the denominator in (3.5), n − q in (3.6) must be understood as n − q + i0 + . As a consequence of helicity conservation of QCD in the high-energy limit, the Lorentz structure in (3.5) gives a non-vanishing contribution only to the left-helicity form factor F L , which can be expressed as The form factors F R and F A vanish at leading power, F LP R = F LP A = 0. The −i0 + prescription in (3.8) generates a rescattering phase of the form factor F L for q 2 > 0.
A factorization formula for NLP corrections is presently not known. Following [5,11] we infer the leading NLP O(α 0 s ) contributions by a diagrammatic analysis of the tree diagrams of Fig. 2. In the hard-collinear region, diagrams (b) and (c) vanish at LP. Their NLP contribution can be expressed in terms of f B . In diagram (a), which gives the O(α 0  in the above LP factorization result, we now expand the light-quark propagator to NLP, and obtain where l denotes the spectator-quark momentum of order Λ QCD and terms suppressed by two powers of Λ QCD /{n + q, m b } are neglected. This expression reduces to the one [5] for real photons, when n − q = 0 and n + q = 2E γ . Proceeding as in [5], we find for the O(α 0 s ) NLP terms of the form factors. These expressions are written in a form such that the complete expression for F L,R including F LP L is valid for both, hard-collinear and hard q 2 . 2 Setting 2v · q = n + q + n − q → n + q and neglecting q 2 in the denominator of the second term in (3.11), one recovers the strict NLP expressions in the hard-collinear region. The q 2 -dependent inverse moment of the B LCDA is defined as Power corrections to F L from the photon emission off the spectator quark cannot be factorized and are parametrized by the "symmetry-preserving", power-suppressed form factor ξ(q 2 , v · q). For q 2 → 0 we recover the result for the on-shell B − → ν γ form factors [5]. One important difference between the virtual and on-shell photon case concerns the second term in the square brackets in (3.9), which matches onto a hadronic matrix element with a transverse derivative acting on the spectator-quark field. This term contributes only to the longitudinal form factor F A and is hence irrelevant in B → γ ν. Since we do not consider explicitly the tree-level contributions proportional to the three-particle LCDAs of the B meson, we can compute this term in the so-called Wandzura-Wilczek (WW) approximation [18], in which case only the subleading two-particle LCDA φ B − (ω) appears. We then find, for hard-collinear q 2 , In the numerical analysis, we employ an expression for F A , which is accurate in both the hard and hard-collinear q 2 region. To this end, we use (2.6) together with and , which was already introduced for B → K * [19], is defined in analogy to (3.12). The finite invariant mass of the virtual photon regulates its endpoint divergence at ω → 0. Nevertheless, in the limit q 2 → 0 we find F A → 0 due to the additional n − q in the numerator, as it should be, since an on-shell photon has no longitudinal polarization. As in the case of F L we allow for a possible non-factorizable contribution by adding an unknown form factor ξ (q 2 , v · q), which must also vanish as q 2 → 0.
The power-suppressed form factor ξ that parameterizes the contribution from soft distances x ∼ 1/Λ QCD between the currents in T µν as well as the three-particle B LCDA contributions have been calculated with light-cone QCD sum rules [17,20], but this method can only be used for q 2 = 0 or space-like. We therefore follow the simple ansatz [11] which incorporates the observation that the power-suppressed form factors appear to reduce the LP ones by setting them to a fraction r LP = 0.2 of F L at tree level. Since there is no LP contribution to F A , ξ is set to 0 in this model. The branching fraction of the four-lepton decay is quite sensitive to the value of r LP . For the B s → µ + µ − γ decay the conservative estimate r LP = 0.2 ± 0.2 was adopted [11]. Below we also present results for r LP = 0.2 ± 0.1.
Since the B → γ * form factors are time-like, the heavy-quark / large-energy expansion is certainly upset locally by the lowest light-meson resonances, ρ and ω. However, as shown in [11], quark-hadron duality is also violated globally, such that for any q 2 bin that contains these resonances, the resonance contribution will be dominant. In order to describe the form factors in the entire region q 2 6 GeV 2 , we add the resonant process to the factorization expressions (3.10) and (3.11). As discussed in [11], this procedure can be justified parametrically, as the averaged resonance contribution is formally a power correction. Nevertheless, the existence of resonances implies that the QCD factorization calculation of the time-like form factors is not on as solid ground as at q 2 = 0 for B − → −ν γ. Writing the dispersion relation in q 2 at fixed n + q for the hadronic tensor T µν , and including only the ρ and ω resonances in the spectral function in the Breit-Wigner approximation, we find where the upper (lower) sign applies to F L (F R ). In addition, with c ρ = 1/2 and c ω = 1/6. 3 For the B → V transition form factors V, A 1 and A 2 we use the definition and numerical results of [21]. It follows from the heavy-quark symmetry relations for the heavy-to-light B → V form factors [18], (applicable since k 2 = k 2 (n + q, q 2 ) via (2.8) and n + q Λ QCD and q 2 = m 2 V in the argument of the form factors in (3.16)) that F res R is power-suppressed relative to F res L , hence it formally counts as a next-to-next-toleading power correction. We do not add a resonance contribution to the form factor F A for the longitudinal intermediate polarization state, since a simple Breit-Wigner ansatz as above would lead to a non-vanishing form factor at q 2 = 0 resulting in a 1/q 4 singularity in the rate, which is unphysical. 4

Numerical results
We combine the form factors at leading-power (LP) and next-to-leading power (NLP) calculated with QCD factorization with the resonance contribution into the final result We include renormalization group evolution to sum logarithms of the ratio of the hard, hard-collinear and soft scales in the LP term following [5], but not in F NLP X where we set F B = f B . Contrary to [11], we do not re-expand products of series expansions in α s . 3 Compared to [11], c ρ has opposite sign because it arises from the uū content of the ρ meson, while in [11] the dd component was the relevant one. 4 We note that such an ansatz has been used in [9,10]. We remark that F res L(R) retains an unphysical imaginary part at q 2 = 0 from (3.17), which, however, is even further suppressed by the small width  [25][26][27][28]. Here we quote the exclusive |V ub | value from HFLAV [22], which uses lattice inputs from [29,30]. Here α em = α We use the inputs specified in Table 1 and the exponential model for the B LCDA. We put ω 0 = λ B ≡ λ + B (n − q = 0) = 0.35 ± 0.15 GeV at the scale 1 GeV as our default value. In the LP terms, we evolve φ B + (ω) to the hard-collinear scale µ hc employing the analytic expression given in [20]. Previous analyses of B → γ ν showed that the shape of the B LCDA is also important when including power corrections [20]. For the time-like virtual photon form factors, there is less control over power corrections and we therefore content ourselves with the exponential model to present our main results and conclusions. We further study the dependence of λ ± B (n − q) and the branching fraction of the four-lepton decay on the shape of the B meson LCDA in Sec. 4.4 using three twoparameter models [20] for the B LCDA. In addition, as |V ub | is an overall factor we do not include its uncertainty in our error estimates, nor do we include the negligible uncertainties on the other input parameters in Table 1. We expect that eventual measurements of the four-lepton final states will be normalized to the decay rate of another, accurately known, exclusive b → u transition.

Form factors
In Fig. 3, we show |F LP L | at leading order (LP,LO) and next-to-leading order (LP,NLO) as a function of q 2 at fixed n + q = 4 GeV. The band describes the scale uncertainty of the hardcollinear scale µ hc = 1.5 ± 0.5 GeV (left) and that of the hard scale µ h = 5 +5 −2.5 GeV (right). Similar to the B → γ ν case, at small q 2 the form factor in the LO approximation has a large scale uncertainty, which is practically removed by the NLO correction. We conclude that the LP form factor is under very good control away from light-meson resonances, once the B LCDA input is specified. It is worth noting that the form factors do not fall off right away with increasing q 2 , but exhibit a maximum near q 2 ≈ 0.5 GeV 2 . The maximum is generated by the sizeable imaginary part πφ B + (n − q) of the q 2 -dependent B LCDA moment λ + B (n − q). These features of F L are largely independent of the chosen value of n + q. The breakdown of |F L | into its various contributions is shown in Fig. 4, starting with LP,NLO, then successively adding the NLP local (loc) contributions (defined as F NLP X without the ξ term), the ξ-contribution as defined in (3.15) and finally the resonance contribution (res). We observe that the NLP contribution is of similar size as the NLO correction at LP. In the small q 2 region, the form factor is locally dominated by the resonance contribution, as expected. However, also at larger q 2 the resonance contribution is comparable to the NLP local contribution. This is due to the fact that the fall-off of the form factors in QCD factorization with increasing q 2 is not faster than the 1/q 2 fall-off of the Breit-Wigner parametrization of the resonances. Note that we have fixed again n + q = 4 GeV, and only show the q 2 dependence of the form factors as the above observations are generic.
The q 2 dependence of the power-suppressed (NLP) form factors F R and F A is shown in the lower panel of Fig. 5. For F R , we show separately the local NLP contribution and the total by adding the resonance contribution. As we do not include a resonance contribution for F A , we only show the total form factor. In addition, we show the dependence on λ B by varying it from 200 MeV (dashed) to 500 MeV (dotted). Except for very small q 2 , the form factor relevant to the longitudinal polarization state of the virtual photon is significantly larger than the one for the right-helicity state.
For First, for all three form factors there is a crossing of the dashed and dotted lines, such that the lower value λ B = 200 MeV increases the form factors at small q 2 but decreases it for large q 2 , while for the upper value λ B = 500 MeV the situation is reversed. In the region where the crossing occurs (around 3.5 GeV 2 for F L ) all sensitivity to λ B is lost. Second, for F L at low q 2 the sensitivity to λ B is larger than the uncertainty coming from r LP . Finally, we comment on the contribution of the three form factors to the differential rate in (2.10). More precisely, we show in Fig. 6 the three terms in the round bracket in (2.10), that is, the form factors squared including their kinematic prefactors. It is remarkable that the longitudinal polarization term λ 2 q 2 |F A | 2 dominates the rate outside the resonance region, despite that fact that it is technically power-suppressed. Moreover, leaving out the resonance term, the longitudinal term would dominate even at small q 2 , although it vanishes for q 2 → 0 (since F A ∼ q 2 as q 2 → 0), while the other two terms approach constants in this limit.
This behaviour can be understood by comparing the analytic expressions for the three terms (without the resonance term) for small q 2 . For small q 2 , the first two terms in the round bracket of (2.10) combine to 16k 2 (m 2 B − k 2 ) 2 (|F L | 2 + |F R | 2 ), and we then estimate where the last line refers to the representative value n + q = 2m B /3. The parametric dependence identifies this ratio as power-suppressed in the hard-collinear region q 2 ∼ m B Λ QCD as it should be. However, the large numerical factor 27π 2 /4 implies that the longitudinal term dominates whenever q 2 is larger than the very small value 0.4 GeV 2 as seen in the Figure. The origin of the large factor is the π 2 that arises from the large imaginary part of the inverse B LCDA moment, in this case λ − B (n − q) in (3.13), since for values q 2 ∈ [0.1, 1] the logarithmic term ln 2 q 2 e γ E n + qλ + B is small.

Predictions for the branching ratios
In this section, we provide theoretical predictions for the branching ratio in various q 2 bins, integrated over n + q (alternatively, k 2 ). The factorization calculation of the form factors in (3.8), (3.10), (3.11) and (3.13) are valid only for n + q ∼ O(m B ). We therefore assume n + q > 3 GeV, which corresponds to E γ = 1.5 GeV at q 2 = 0, and integrate the doubledifferential branching fraction over n + q > 3 GeV before forming q 2 bins. A rough estimate, obtained by assuming that our results apply in the full phase space, shows that the n + q cut reduces the rate by O(20%) for the [1.5, 6] GeV 2 q 2 bin.
For non-identical lepton flavours, = , the required n + q cut can easily be applied as for each event n + q can be inferred from the reconstructed k 2 and q 2 using (2.8). For the q 2 bins, we consider the low bin [4m 2 µ , 0.96] GeV 2 , where the upper boundary of the bin is determined such that the large experimental background from φ mesons decaying into a lepton pair is avoided. This bin was also considered by the LHCb Collaboration [8]. Fig. 5 shows that in this bin, the ρ and ω resonances make a large contribution. As mentioned, we do not attribute an additional error due to our resonance model. This introduces an additional uncertainty in this region which is challenging to quantify. Above q 2 > 1 GeV 2 , the effect of the ρ and ω resonances (and thus a possible uncertainty associated with this) is significantly reduced. We consider three different q 2 bins: [1,6], [1.5, 6] and [2,6] GeV 2 . In these bins, the resonance contribution is approximately 10% only. In Table 2, we give the branching ratio in these q 2 bins, specifying the contributions which are successively added. In addition, we specify the uncertainties from variations of the scales µ h,hc , r LP = 0.2 ± 0.2 and λ B = 350 ± 150 MeV. We observe that in the three considered regions above q 2 > 1 GeV 2 , the effect of the resonances is smaller than the uncertainty from r LP .

Identical lepton flavours
A challenge arises when considering identical lepton flavours, = , because experimentally the two like-sign leptons cannot be distinguished. This results in the additional interference term (2.11). More challenging is the required cut on n + q, where q is the photon momentum, to ensure that the photon has hard-collinear momentum. Considering again B − → − (p 1 ) + (p 2 ) − (p 3 )ν (p ν ), with q 2 = (p 1 + p 2 ) 2 andq 2 = (p 2 + p 3 ) 2 , experimentally, q 2  Table 2: Branching ratio for the two non-identical lepton flavour cases (in 10 −8 ) integrated over different bins in q 2 and for n + q > 3 GeV. We show the individual contributions consecutively adding to the LP result the NLP local and ξ contributions and finally the resonances.
In addition, we quote the uncertainties from varying the scales µ h,hc , r LP = 0.2 ± 0.2 and λ B = 350 ± 150 MeV. The total uncertainty is obtained by adding them in quadrature. For electrons, we also consider a low bin with q 2 min = 0.0025 GeV 2 .
andq 2 cannot be distinguished. Instead, the invariant mass of two µ − µ + -pairs are defined as q 2 low < q 2 high . In this case, placing the required cut on n + q is not unambiguously possible as we cannot determine if the virtual photon has q 2 low or q 2 high associated with its momentum. To deal with this issue, several observations can be made: • for small q 2 low , the photon momentum can be associated with q 2 low most of the time. If this is the case, a cut on n + q low > 3 GeV suffices (similar to the non-identical lepton flavour case). In fact, a more detailed analysis shows that the cases falling outside this cut (i.e. the region which cannot be described in QCD factorization in which the photon has q 2 high but n + q small) is phase-space suppressed by two powers of 1/m b compared to the leading contribution.
• for q 2 bins above 1 GeV 2 , the situation is more complicated as the photon more often has q 2 high . Therefore, we have to ensure n + q > 3 GeV for both q 2 low and q 2 high .
We thus have to restrict both n + q low > 3 GeV and n + q high > 3 GeV. These quantities are now defined in the following way: 5 For each event, we specify q 2 low and q 2 high . We can then associate the remaining lepton plus neutrino as k 2 low and k 2 high , respectively. Here high and low are just labels and in this case k 2 low is not necessarily lower than k 2 high . Then using (2.8), both n + q low and n + q high can be calculated from their corresponding k 2 and q 2 . Alternatively, one could cut on k 2 low and k 2 high directly. 5 We remark that, unlike previously, here n + q high does not coincide with the component of q µ high in the n µ − direction as defined above (2.8), if the momentum of the γ * is q µ low . The reason is that we always align the z-axis with the three-momentum of the γ * , but we do not know which of q µ low and q µ high refers to the virtual photon momentum.  Table 3: Branching ratio for the two identical lepton cases (in 10 −8 ) integrated over different bins in q 2 low applying two cuts: n + q low > 3 GeV and n + q high > 3 GeV. We show the individual contributions consecutively adding to the LP result the NLP local and ξ contributions and finally the resonances. In addition, we quote the uncertainties from varying the scales µ h,hc , r LP = 0.2 ± 0.2 and λ B = 350 ± 150 MeV. The total uncertainty is obtained by adding these contribution in quadrature. For the total results, we also quote the result with only one cut: n + q low > 3 GeV in parenthesis. For electrons, we also consider a low bin with q 2 min = 0.0025 GeV 2 .
For our final results in Table 3, we thus include two cuts: n + q low and n + q high > 3 GeV for all bins. Our final results for the branching ratio for different q 2 low bins are given in Table 3. Again we present the different contributions added successively. We emphasize that placing these two cuts on n + q might be conservative, specifically for the low q 2 bin as discussed above, given the phase space suppression of the region in which q 2 high is associated with the photon. We confirm numerically that indeed this region is small, by calculating the rate with and without the cut on n + q high . For comparison, in Table 3 we also give the results for the total rate with only the n + q low cut in parenthesis. For identical lepton flavours, the branching ratio contains two contributions as defined in (2.11). With this convention, we find that Br int contributes positively to the rate but is suppressed by at least one order of magnitude compared to the non-identical lepton flavour rate.
A comment on the low-q 2 bin for B − → µ − µ + µ −ν µ is in order. Our prediction for Br(B + → µ − µ + µ −ν µ ) is 1.54 (1.77) · 10 −8 and includes cuts on n + q. Yet, it lies close to the upper limit < 1.6 · 10 −8 given by the LHCb collaboration for this decay mode in this bin [8]. Hence the LHCb result may already point towards a larger value of λ B .

Sensitivity to λ B
Our predictions for the branching ratio suffer from a large uncertainty due to λ B . Therefore, a measurement of the branching ratio may be used to obtain a bound on λ B . Figure 7 shows the rate as a function of λ B for the [4m 2 µ , 0.96] GeV 2 bin and for the [1,3] and [3, 6] GeV 2 bins. We have split the [1,6] GeV 2 to avoid integrating over the region where the sensitivity to λ B variation switches sign (see Fig. 5). As the uncertainty of the branching : Branching ratio (in 10 −8 ) for three different q 2 bins as a function of λ B for the decay mode B − → µ − µ + e −ν e . We include the variation from µ h,hc and from δr LP as inner (red) and outer (blue) bands, respectively. ratio is dominated by the error on r LP , we consider two options; δr LP = 0.2 and δr LP = 0.1. The latter option shrinks the error by half. These predictions include our model for the long-distance resonance contributions as described above, for which we do not add an uncertainty. We note that the sensitivity to λ B is best for the small q 2 bins, while it is significantly reduced for higher q 2 bins. Comparing to B → γ ν [20], we conclude that the sensitivity to λ B for B → ν ¯ in the low-q 2 bin is comparable (compare to Fig. 5 in [5]) for λ B < 200 MeV, but less when it is larger. However, in this bin the resonance contribution is sizeable and there is unquantified model dependence related to its interference with the factorization contribution. For the [1,3] GeV 2 bin, the resonance contribution is less pronounced and thus this bin could still provide information on λ B despite its smaller sensitivity.

Dependence on the shape of the B LCDA
Up to now, we used the exponential model (4.2) to present our main results. However, it is known that for B → γ ν [20] the shape of the B meson LCDA has a significant effect through the dependence of radiative corrections on the logarithmic inverse moments, and of the power-suppressed form factor ξ through its dependence on the shape of the LCDA in the sum rule calculation. In four-lepton decay the generalized inverse moments λ ± B (n − q) introduce further dependence on the shape of the LCDA.
To study this dependence, we consider three two-parameter models [20] φ I where U (α, β, z) is the confluent hypergeometric function of the second kind. Given (ω 0 , a) one determines λ B and the dimensionless shape parameter σ 1 , related to the first inverselogarithmic moment. The range of a is chosen such that the range −0.693147 < σ 1 < 0.693147 is covered, where σ 1 = 0 in the exponential model, see [20] for more details. To study the influence of the shape of the LCDA, we are then interested in the envelope of theoretical predictions of all three models spanned by the variation of a for given λ B . For simplicity, we assume that these forms of the B meson LCDA hold at the scale µ hc = 1.5 GeV, so that no renormalization group evolution to the hard-collinear scale needs to be performed. We obtain φ − (ω) using the Wandzura-Wilczek (WW) relation [18] φ − (ω) = ∞ ω dω ω φ + (ω ) . (4.5) The n − q dependent moments λ ± B (n − q) are then obtained using (3.12) (and equivalently for λ − B (n − q)). Again we define λ B ≡ λ + B (n − q = 0), such that ω 0 can be related to λ B via In Figures 8 and 9, respectively, we show the q 2 dependence of 1/λ + B (n − q) and 1/λ − B (n − q) for fixed n + q = 4 GeV for the three B LCDA models by varying the parameter a within the ranges indicated in (4.4), fixing λ B = 350 MeV. The black solid line represents the exponential model. There is a significant dependence of 1/λ ± B (n − q) on the B-meson LCDA shape -this is expected, as for instance, the q 2 dependence of the imaginary part 1/λ ± B (n − q) is directly related to the ω-dependence of φ ± (ω).
Finally, we compute the effect of the B meson LCDA shape on the branching ratio. In Figure 10, we show the dependence of the branching ratio in the [4m 2 µ , 0.96] GeV 2 q 2 bin on λ B and the shape parameter a for the three models. For given λ B on the horizontal  axis, the bands are obtained by varying a in its allowed range. In black, we also show the exponential model. Comparing with our previous results, we observe that the dependence on the shape is about as large as the dependence on the r LP variation from the powersuppressed form factor, see Figure 7. The conclusion is thus similar to the case of B → γ ν [20]. Once sufficient data is available, a correlated determination of λ B together with the shape parameter σ 1 (and, perhaps, others) should be performed. The important point is that the predicted branching fractions are highly sensitive to B meson LCDA input, even if not necessarily λ B alone.

Conclusion
Motivated by the first search and upper limit [8] for the rare charged-current B decay to a four-lepton final state ν ( )¯ ( ) , this work considered the calculation of the decay amplitude with factorization methods. Combining methods previously applied to B − → −ν γ [5] and B s → µ + µ − γ [11], we obtain the B → γ * form factors, which depend on the invariant masses of the two lepton pairs, in QCD factorization at next-to-leading order in α s and leading power in an expansion in Λ QCD /m b , and to leading order in α s at next-to-leading power. To this we added a simple Breit-Wigner parametrization of the ρ, ω intermediate resonances.
Although suppressed beyond next-to-leading power, the resonances dominate the spectrum in the ( )¯ ( ) invariant mass q 2 locally, making the predictions more uncertain in this region than at large invariant mass or for B − → −ν γ. Quite generally it must be noted that the parametric counting that justifies the heavy-quark expansion is not well respected, as is evidenced by the large contribution of the longitudinal polarization state of the intermediate virtual photon.
Our calculations predict branching fractions of a few times 10 −8 in the q 2 bin up to 1 GeV 2 , which are accessible to the LHC experiments. The branching fraction rapidly drops with increasing q 2 , reaching 10 −9 in the bin [1.5, 6] GeV 2 .
Confronting these results to measurements checks our understanding of Standard Model dynamics in these rare decays. An important further motivation for this investigation has been to explore the sensitivity of the decay rate to the inverse moment λ B of the leadingtwist B meson light-cone distribution amplitude. For non-vanishing q 2 the access to λ B is less direct than in B − → −ν γ, and requires some knowledge of the shape of the LCDA as well. At large q 2 , the sensitivity disappears. We find these expectations confirmed in Fig. 7, which shows that λ B is best determined from the small-q 2 bin. In this bin the sensitivity to λ B is almost comparable to B − → −ν γ when λ B < 200 MeV, but less when it is larger. However, one should be aware that in this bin the resonance contribution is sizeable and there is unquantified model dependence related to its interference with the factorization contribution. The q 2 bin above 1 GeV 2 can still yield useful bounds on λ B , despite its weaker sensitivity. As for the case of B → γ ν [20], once sufficient data is available, a correlated determination of λ B together with B meson LCDA shape parameters should be performed. Overall, we conclude that the four-lepton final state cannot fully replace the B − → −ν γ mode to measure λ B . However, given the current state of knowledge, any complementary experimental result on λ B is worthwhile pursuing.

Note added
While this paper was being finalized, Ref. [31] appeared. We note the following important differences: (1) The third independent form factor, related to F A , is missed, see App. A.
(2) The q 2 distribution is computed without a cut on n + q, hence includes significant phasespace regions where the adopted QCD factorization treatment is not applicable. (3) The residual scale dependence of the leading-power form factor at NLO in the strong coupling is much larger than ours. Presumably this is because it is assumed (incorrectly) that the form of the exponential model is preserved by renormalization group evolution. (4) Only the region of small q 2 < 1 GeV 2 is discussed. In this region our results are dominated locally by the Breit-Wigner parameterization of the ρ and ω resonances, whereas [31] adopts the QCD sum rule expression [32] for the power-suppressed form factor ξ, but in the timelike region. (5) For the case of identical lepton flavours, we present partially integrated branching fractions that correspond to experimental observables.

A Decomposition of the hadronic tensor
Using Lorentz covariance only, the most general decomposition of the hadronic tensor T µν (p, q) defined in (2.2) contains six independent scalar form factors F 1...6 = F 1...6 (k 2 , q 2 ): T µν = F 1 g µν + F 2 µναβ k α q β + F 3 k µ q ν + F 4 q µ k ν + F 5 k µ k ν + F 6 q µ q ν . (A.1) For γ * →¯ and W * → ν , the q µ (k ν ) terms do not contribute to the decay amplitude if ( ) is massless. The Ward identity q µ T µν = f B (k + q) ν implies the relations and hence reduces the number of independent form factors to four. We choose to eliminate F 3 and F 5 , and write T µν = F 1 g µν + F 2 µναβ k α q β + f B − F 1 − q 2 F 6 k · q k µ q ν + F 4 q µ k ν + f B − q 2 F 4 k · q k µ k ν + F 6 q µ q ν . (A.3) Since we consider massless leptons we now can drop all terms in the second line, which leaves three independent form factors. The number of independent form factors can be associated with the number of independent polarization states of the virtual photon. Note that dropping F 4,5,6 in (A.1) before applying the Ward identity would lead to the omission of the F 6 term in the coefficient of the k µ q ν term in (A.3), and to the wrong conclusion (since F 6 does not vanish) that there are only two independent form factors F 1,2 for massless leptons. 6 It is straightforward to work out the relations between the form factors F 1,2,6 and F A ⊥ ,V,A used in the main text.