Study of $B_s\to\ell^+\ell^-\gamma$ decays in covariant quark model

We study the rare radiative leptonic decays $B_s\to\ell^+\ell^-\gamma$ ($\ell=e,\mu,\tau$) within the Standard Model, considering both the structure-dependent amplitude and bremsstrahlung. In the framework of the covariant confined quark model developed by us, we calculate the form factors characterizing the $B_s\to\gamma$ transition in the full kinematical region of the dilepton momentum squared and discuss their behavior. We provide the analytic formula for the differential decay distribution and give predictions for the branching fractions in both cases: with and without long-distance contributions. Finally, we compare our results with those obtained in other approaches.


I. INTRODUCTION
The rare radiative decays B q → ℓ + ℓ − γ with ℓ = e, µ, τ and q = d, s are of great interest for several reasons. First, they are complementary to the well-known decays B → K ( * ) ℓ + ℓ − and therefore, provide us with extra tests of the Standard Model (SM) predictions for processes which proceed at loop level. Second, this process is not helicity suppressed as compared with the pure leptonic decays B q → ℓ + ℓ − due to the appearance of a photon in the final state.
Theoretical estimates of the decay branching fractions have shown that B (B s → µ + µ − γ) may be an order of magnitude larger than B (B s → µ + µ − ).
There is a number of theoretical calculations of the branching fractions B(B q → ℓ + ℓ − γ) performed in different approaches. Among them one can mention the early studies in the framework of: a constituent quark model [1], light-cone QCD sum rules [2,3], and the lightfront model [4]. The structure-dependent amplitude of the decays B s → ℓ + ℓ − γ was analyzed in Ref. [5] by taking a universal form for the form factors, which is motivated by QCD and related to the light-cone wave function of the B s meson. In Ref. [6] it was shown that efficient constraints on the behavior of the form factors can be obtained from the gauge-invariance requirement of the B q → ℓ + ℓ − γ amplitude, as well as from the resonance structure of the form factors and their relations at large photon energies. Universality of nonperturbative QCD effects in radiative B decays was studied in Ref. [7]. In Ref. [8] long-distance QCD effects in the B d/s → ℓ + ℓ − γ decays were analyzed. It was shown that the contribution of light vector-meson resonances related to the virtual photon emission from valence quarks of the B meson gives sizable impact on the dilepton differential distribution. In Ref. [9] the B d/s → γ transition form factors were calculated within the relativistic dispersion approach based on the constituent quark picture. A detailed analysis of the charm-loop contributions to the radiative leptonic decays was also performed. Very recently, a novel strategy to search for the decays B s → µ + µ − γ in the event sample selected for B s → µ + µ − searches was presented [10].
In this paper we calculate the matrix elements and the differential decay rates of the decays B s → ℓ + ℓ − γ in the framework of the covariant confined quark model previously developed by us (see, e.g., Ref. [11]). This is a quantum-field theoretical model based on relativistic Lagrangians which effectively describe the interaction of hadrons with their constituent quarks. The quark confinement is realized by cutting the integration variable, which is called the proper time, at the upper limit. The interaction with the electromagnetic field is introduced by gauging the interaction Lagrangian in such a way to keep gauge invariance of the matrix elements at all calculation steps. This model has been successfully applied for the description of the matrix elements and form factors in the full kinematical region in semileptonic and rare decays of heavy mesons as well as baryons (see, e.g., Refs. [12][13][14][15][16][17]).
The rest of the paper is organized as follows. In Sec. II we give a necessary brief sketch of our approach. The introduction of electromagnetic interactions in the model is described in Sec. III. Sec. IV is devoted to the calculation of the decay matrix elements. We also briefly discuss their gauge invariance. In Sec. V we recalculate the formula for the twofold decay distribution in terms of the Mandelstam variables (t, s). Then we integrate out the t variable analytically and present the expression for the dilepton differential distribution. In Sec. VI we provide numerical results for the form factors, the differential decay widths, and the branching fractions. A comparison with existing results in the literature is included.
Finally, we conclude in Sec. VII.

II. BRIEF SKETCH OF THE COVARIANT CONFINED QUARK MODEL
The covariant confined quark model (CCQM) has been developed by our group in a series of papers. In this section, we mention several key elements of the model only for completeness. For a more detailed description of the model, as well as the calculation techniques used for the quark-loop evaluation, we refer to Refs. [11][12][13][14][15][16][17][18] and references therein.
In the CCQM, the interaction Lagrangian of the B s meson with its constituent quarks is constructed from the hadron field B s (x) and the interpolating quark current J Bs (x): where the latter is given by The hadron-quark coupling g Bs is obtained with the help of the compositeness condition, which requires the wave function renormalization constant of the hadron to be equal to zero Z H = 0. Here, F Bs (x; x 1 , x 2 ) is the vertex function whose form is chosen such as to reflect the intuitive expectations about the relative quark-hadron positions where we require w 1 + w 2 = 1. We actually adopt the most natural choice in which the barycenter of the hadron is identified with that of the quark system. The interaction strength Φ Bs (x 1 − x 2 ) 2 is assumed to have a Gaussian form which is, in the momentum representation, written as Here, Λ Bs is a hadron-related size parameter, regarded as an adjustable parameter of the model. Additional free parameters include the constituent quark masses and a universal cutoff parameter λ cutoff related to the embedded infrared confinement of constituent quarks inside hadrons. In this paper we use the results of the updated fit performed in Refs. [16,19,20]. The model parameters involved in this paper are given by (in GeV) Once these parameters are fixed, one can employ the CCQM as a frame-independent tool for hadronic calculation. The parameters are determined from a least-squares fit to available experimental data and some lattice calculations. We have observed that the errors of the fitted parameters are within 10%. Besides, based on a comparison of our predictions with experimental data and lattice QCD in previous studies, we estimate the theoretical error of the CCQM to be of the order of 10%.

III. ELECTROMAGNETIC INTERACTIONS
Within the CCQM framework, interactions with electromagnetic fields are introduced as follows. First, one gauges the free-quark Lagrangian in the standard manner by using minimal substitution that gives the quark-photon interaction Lagrangian In order to guarantee local invariance of the strong interaction Lagrangian, one multiplies each quark field q(x i ) in L str int with a gauge field exponential. One then has where The path P connects the endpoints of the path integral.
It is readily seen that the full Lagrangian is invariant under the transformations One then expands the gauge exponential up to the required power of e q A µ needed in the perturbative series. This will give rise to a second term in the nonlocal electromagnetic interaction Lagrangian L em−nonloc int . At first glance, it seems that the results will depend on the path P taken to connect the endpoints of the path integral in Eq. (10). However, one needs to know only the derivatives of the path integral expressions when calculating the perturbative series. Therefore, we use the formalism suggested in Refs. [21,22], which is based on the path-independent definition of the derivative of I(x, y, P ) As a result of this rule, the Lagrangian describing the nonlocal interaction of the B s meson, the quarks, and electromagnetic fields reads (to the first order in the electromagnetic charge) where p = w 2 p 1 + w 1 p 2 .

IV. MATRIX ELEMENTS OF THE DECAYS
The decays B s → ℓ + ℓ − γ are described by three sets of diagrams shown in Figs. 1-3.
Diagrams from the first set ( Fig. 1) correspond to the case when the real photon is emitted from the quarks or the meson-quark vertex. The effective Hamiltonian describing the b → sℓ + ℓ − weak transition is written as where λ t = V tb V * ts , andm b is the QCD quark mass which is different from the constituent quark mass m b used in our model. Here and in the following we denote the QCD quark masses with a tilde to distinguish them from the constituent quark masses used in the model (see Eq. (6)). The Wilson coefficients C eff 7 = C 7 − C 5 /3 − C 6 and C 10 depend on the scale parameter µ. The Wilson coefficient C eff 9 effectively takes into account, first, the contributions from the four-quark operators O i (i = 1, ..., 6) and, second, nonperturbative effects coming from the cc-resonance contributions which are as usual parametrized by the Breit-Wigner ansatz [23]: where iπ.
We will use the value of µ =m b pole for the renormalization scale. The numerical values for the model-independent input parameters and the corresponding Wilson coefficients are given in Table I. The SM Wilson coefficients are taken from Ref. [24]. They were computed at the matching scale µ 0 = 2M W , and run down to the hadronic scale µ b = 4.8 GeV. The evolution of couplings and current quark masses proceeds analogously.  [24]. Bare quark massesm q are given in GeV.
Diagrams from the first two sets contribute to the structure-dependent (SD) part of the decay amplitude. They can be parametrized by a set of invariant form factors. In order to define the form factors, we specify our choice for the momenta in the decays as follows: where p 1 = p 2 + k + + k − , p 1 − p 2 = k + + k − ≡ q, with p 2 1 = M 2 Bs , p 2 2 = 0, ǫ † 2 · p 2 = 0, and k 2 + = k 2 − = m 2 ℓ . We will use the definition of the B s → γ transition form factors given, for instance, in Ref. [8] Here we use the short notations σ µq ≡ σ µβ q β and ε µαp 1 p 2 ≡ ε µαβδ p 1β p 2δ . The The process with the virtual photon emitted from the light s quark is described by the diagram in Fig. 2 (left panel). This diagram has a branch point at q 2 = 4m 2 s if one does not cut the upper limit of the integration over proper time. The infrared cutoff eliminates this singularity and thereby guarantees the quark confinement. However, the physical region for Bs , which is too far from the value q 2 = 4m 2 s . In this case, the form factorsF s(lℓ)s T V /T A become very large and unreliable. This is a drawback of our model. We are planning to modify our approach in order to avoid such unphysical growths, and then include the intermediate φ resonance into this diagram. For the time being we drop the form factorsF s(lℓ)s T V /T A from the consideration. Finally, the SD part of the amplitude is written in terms of the form factors as follows: where T µα 1 ≡ (g µα p 1 p 2 − p α 1 p µ 2 ). The structure-independent part of the amplitude (Bremsstrahlung) is described by the diagrams in Fig. 3. Only the operator O 10 contributes to this process, and effectively gives the leptonic decay constant f Bs . One has Here, t = ( where the bounds t ± are given by One can see that t ± = m 2 ℓ at minimum recoil s = q 2 = M 2 Bs , which leads to the infrared pole in Eq. (22). A cut in the photon energy is required. In the B s center-of-mass system one has Note that there are also weak annihilation diagrams with u(c) anomalous triangle in addition to the above diagrams. However, the contribution from these diagrams is much smaller than that from other diagrams as was shown in Ref. [8]. Therefore, we will drop this type of diagrams in what follows.
, the amplitude has also the non-gauge-invariant pieces g µα and p µ 1 p α 1 . We have checked numerically that the form factors corresponding to the non-gauge-invariant part vanish for arbitrary momentum transfer squared q 2 .

V. DIFFERENTIAL DECAY RATE
The twofold decay distribution is written as where pol denotes the summation over polarizations of both the photon and leptons. The physical region was discussed in the previous section, which reads 4m 2 ℓ ≤ s ≡ q 2 ≤ M 2 Bs , and t − ≤ t ≤ t + . It is more convenient to write the final result for the twofold decay distribution in terms of dimensionless momenta and masses: The decay distribution is written as a sum of the structure dependent (SD) part, the Bremsstrahlung (BR), and the interference (IN) ones as follows: One has We have checked that all the expressions above are in agreement with those given in Ref. [9].
After integrating out the variablet we obtain the following analytic expressions for the differential decay rate where β = 1 − 4m 2 ℓ /ŝ.

VI. NUMERICAL RESULTS
In Fig. 4 we present the q 2 dependence of the calculated form factors in the full kinematical region 0 ≤ q 2 ≤ M 2 Bs . The results of our numerical calculations can be approximated with high accuracy by the double-pole parametrization with the relative error less than 1%. The parameters F (0), a, and b are listed in Table II.
For completeness we also list here the values of the form factors at zero recoil (q 2 max ). As mentioned in Sec. II, the estimated theoretical error on the form factors is about 10%.  factors calculated in Ref. [9]. Using the definitions in Eqs. (19) and (20) we can relate our form factors F i (q 2 ) to the KMN form factors F i (q 2 , 0) as follows (see Ref. [9] for more detail): One can see that in the low-q 2 region (q 2 20 GeV 2 ) the corresponding form factors from the two sets are very close. In the high-q 2 region, the KMN form factors steeply increase and largely exceed our form factors.
In order to have a better picture of the behavior of the form factors, in Fig. 6 we plot all of them together and compare with those from Ref. [9]. It is very interesting to note that our form factors share with the corresponding KMN ones not only similar shapes (especially  In Fig. 7 we plot the differential branching fractions 10 9 dB (B s → γℓ + ℓ − ) /dŝ as functions of the dimensionless variableŝ = q 2 /M 2 1 . We also plot here the ratio which is a promising observable for testing lepton flavor universality (LFU) in these channels [25]. The ratio r γ is very close to unity in low-q 2 region but far above unity at large q 2 due to Bremsstrahlung. As was pointed out in Ref. [9], in the high-q 2 region (q 2 15 GeV 2 ), the ratio r γ is mainly described by the form factors F A (q 2 ) and F V (q 2 ). Therefore, knowledge of their behavior at large q 2 plays an important role in testing LFU. In Fig. 8 we plot the ratio r γ at large q 2 in comparison with Ref. [9]. The ratios are very close in the range 16 q 2 20 GeV 2 (i.e. 0.55 ŝ 0.69), which is a result of the similarity between the form factors discussed above. The authors of Ref. [25] suggested a useful observable with the optimal choiceŝ 1 = 0.55 andŝ 2 = 0.8 (corresponding to q 2 1 = 15.8 GeV 2 and q 2 2 = 23.0 GeV 2 , respectively). In this range, the ratio is dominated by the form factors F A (q 2 ) and F V (q 2 ). We provide our prediction for this ratio R γ (0.55, 0.8) = 1.46, about 20% larger than the prediction R γ (0.55, 0.8) = 1.115 ± 0.030 given in Ref. [25]. Note that from the results of Ref. [9], one obtains R γ (0.55, 0.8) = 1.32.
In Table III we give the values of the branching fractions calculated without and with long distance contributions. In the calculation with long distance contributions, the region of two low lying charmonia is excluded by assuming 0.33 ≤ŝ ≤ 0.55, as usually done in experimental data analysis. It is seen that the Bremsstrahlung contribution to the electron mode is negligible, while for the tau mode it becomes the main part.  Finally, in Table IV we compare our results for the branching fractions with those obtained in other approaches. In most cases, our predictions for the light leptons are smaller than those from other studies. This is due mainly to the fact that in this paper we drop the contribution from the diagram with the virtual photon emited from the s quark, which should has a sizable impact on the final results. In the case of the tau mode, the contribution from Bremsstrahlung dominates the decay branching fractions and the SD amplitude becomes less important. As a result, our prediction for B(B s → γτ + τ − ) agrees well with other studies.

VII. SUMMARY AND CONCLUSIONS
We have studied the rare radiative leptonic decays B s → γℓ + ℓ − (ℓ = e, µ, τ ) in the framework of the covariant confined quark model. The relevant transition form factors have been obtained in the full kinematical range of dilepton momentum transfer squared. We found a very good agreement between our form factors and those from Ref. [9], especially in the region q 2 20 GeV 2 . We have provided predictions for the decay branching fractions and their ratio. The branching fractions for the light leptons are lower than most existing values in the literature. In our opinion, the main reason is that we have not considered the contribution from the diagram with the virtual photon emitted from the light s quark of the B s meson, due to the involved unphysical growths. However, we note that this diagram has no effects on the form factors F V (q 2 ) and F A (q 2 ) calculated in the present study. Therefore, the two form factors may be useful to study LFU in these rare channels, especially in the high-q 2 region. In our future works, we plan to modify our approach to treat this problem and include the diagram contribution into the process.