Radiative lepton flavor violating B, D, and K decays

We argue that radiative lepton flavor violating (RLFV) decays $P \to \gamma \ell_1 \overline{\ell}_2$ of $P =B^0_q$, $\bar{D}^0$, and $K^0$ meson states are robust probes of new physics models. In particular, they could be used to put constraints on the Wilson coefficients of effective operators describing lepton flavor-changing neutral current interactions at low energy scales. We set up a generic framework for describing these transitions and review new physics constraints from $P \to \ell_1 \bar \ell_2$ decays. There is discussion of how RLFV transitions provide access to the operators that cannot be constrained in two-body decays and we in turn motivate further experimental searches via these channels.


I. INTRODUCTION
Currently operating and future B-factories, such as LHCb and Belle-II, will be accumulating significant amounts of beauty and charm decay data. These large data sets will be quite useful in studies of extremely small decay rates of B and D mesons, which could probe new physics (NP) at unprecedentedly high energy scales. In particular, studies of pseudoscalar meson decays P = B 0 q ,D 0 , and K 0 into the final states containing charged leptons of different flavors such as P → 1 2 and P → γ 1 2 could be performed. Such decays are induced by the operators that generate flavor-changing neutral currents (FCNC) in the lepton sector, which provide a fruitful approach to probing beyond the standard model (BSM) physics, assuming of course that such flavor-violating interactions are allowed in the BSM models. There are indeed many well-established new physics models (see, e.g. [1][2][3][4]) that meet this opportunity and predict charged lepton flavor violating (CLFV) transition rates that are significantly larger than the standard model (SM) rates [1].
A convenient way to describe CLFV transitions in low energy experiments is by introducing an effective Lagrangian, L eff . Such a Lagrangian is a convenient parameterization of all new physics models that include lepton flavor violation with the details of the models encoded in the Wilson coefficients (WCs) of L eff , which are obtained by matching the effective Lagrangian to a given BSM model at the new physics scale Λ [5]. This Lagrangian is required to be invariant under the unbroken symmetry groups SU (3) c × U (1) em below the electroweak symmetry breaking scale. At the low scale for which a given process occurs the effective operators would exhibit the relevant standard model (SM) degrees of freedom with the effective operators written completely using quarks (q i = b, c, s, u, and d) and leptons ( i = τ, µ, and e). In what follows, we assume that top quarks are integrated out of the theory, and we do not consider neutrinos. The effective Lagrangian L eff that involves CLFV can be written as where L D is a dipole part, and L q is the part that contains four-fermion interactions. Since here we are interested in the decays of electrically-neutral pseudoscalar B 0 q ,D 0 , and K 0 mesons to flavor-off-diagonal lepton pairs and other particles, the transitions involve FCNC interactions on both quark and lepton sides.
Note that it is known that the quark FCNC transitions, at least in the decays of downtype quarks, are dominated by the SM contributions. For instance, the dipole operator describing q 1 → q 2 γ can be written as [7] L peng = Here λ P q = V qq 2 V * qq 1 denotes the appropriate Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, m q 1 is the heavier quark, and C 7γ is the corresponding Wilson coefficient [7].
Here m q H is the mass of the heavier quark (m q H = max[m q 1 , m q 2 ]) and P R,L = (1±γ 5 )/2 is the right (left) chiral projection operator. In general the Wilson coefficients would be different for different lepton flavors i and quark flavors q i . Note that, contrary to some previous studies, we include tensor operator in Eq. (4) (see [8] for motivation). CP-conservation is assumed so all the Wilson coefficients in Eq. (4) should be viewed as real numbers.
In this paper we discuss the possibility of the Wilson coefficients of the effective Lagrangian in Eq. (1) for different i and q i be determined from experimental data on leptonic and radiative leptonic CLFV decays of B 0 q ,D 0 , and K 0 states. We review two-body decays P → 1 2 in Sect. III. We will note that restricted kinematics of the two-body transitions would allow us to select operators with particular quantum numbers significantly reducing the reliance on the single operator dominance assumption [8]. The main part of the paper, Sect. III, will be devoted to discussion of radiative lepton-flavor violating (RLFV) decays P → γ 1 2 . We will summarize our results in Sect. IV and conclude in Sect. V.
Note that here we only consider short distance effects in Kaon decays. In the SM long distance effects on decays such as K 0 L(S) → γ ¯ dominate the dynamics [9]. In light of this, our kaon results may be modified by long distance effects.
In what follows, we will use the convention that the subscript of "1" will denote the lighter lepton and the subscript "2" will denote the heavier lepton. Unless otherwise specified when studying the branching ratios we assume for a meson, P , that B (P → (γ) 1 2 ) = B P → (γ) 1 2 + B P → (γ) 1 2 . Finally, it is important to note that some of the twobody and all of the three-body transitions have yet to be experimentally studied. Numerical constraints on some Wilson coefficients of the effective Lagrangian, L ef f , from these unstudied decays are not available.

II. TWO-BODY DECAYS P → 1 2
Many studies have focused on rare leptonic decays of B 0 q mesons, B q → ¯ , as both precision tests of the SM and as an opportunity to search for new physics (e.g. [10][11][12][13][14]).
The abundance of produced B 0 q andD 0 states at the LHCb, Belle II, and BESIII experiments also allows for studies of lepton-flavor violating decays at these experiments [15,16]. Such decays were discussed at length previously, mainly in the context of particular models. Here we shall review these transitions emphasizing the possibility to constrain Wilson coefficients of the axial and pseudoscalar operators of the effective Lagrangian in Eq. (1). These decays would provide information about C q 1 q 2 1 2 P L (C q 1 q 2 1 2 P R ) and/or C q 1 q 2 1 2 AL (C q 1 q 2 1 2 AR ) in Eq. (4).
One can write the most general expression for the P → 1 2 decay amplitude as [8] A(P → 1 2 ) = u(p 1 , s 1 ) E q 1 q 2 1 2 with E q 1 q 2 1 2 P and F q 1 q 2 1 2 P being dimensionless constants which depend on the Wilson coefficients of operators in Eq. (1) and various decay constants.
The amplitude of Eq. (5) leads to the branching ratio for flavor off-diagonal leptonic decays of pseudoscalar mesons:  [16][17][18][19][20]. Center dots signify that no experimental data are available; "FPS" means that the transition is forbidden by phase space.  [21,22], total decay widths, and meson masses [16] used in the calculation of branching ratios B(P → 1 2 ). Here Γ P is the total width of the pseudoscalar state. We have once again neglected the mass of the lighter lepton and set y = m 2 /m P . Calculating E q 1 q 2 1 2 P and F q 1 q 2 1 2 P for P = B 0 d (q 1 q 2 = db), B 0 s (q 1 q 2 = sb),D 0 (q 1 q 2 = cu), and , K 0 L (q 1 q 2 = ds), the coefficients are The hadronic matrix element in Eq. (7) is defined as [23] 0|q Here p is the momentum of the meson. The constant κ P is 1 for B 0 q ,D 0 , and K 0 ; and 1/ √ 2 for K 0 L(S) . The experimental limits and numerical values of the pseudo-scalar decay constants used in the calculations can be found in Tables I and II. The resulting constraints on the Wilson coefficients are found in Table III. Particle masses and other input parameters are from [16][17][18][19][20].

Leptons Initial state
Wilson coefficient Similarly to the B 0 s → µ + µ − γ transition [24][25][26][27][28], addition of a photon to the 1 2 final state allows one to probe operators of the effective Lagrangian that do not contribute to P → 1 2 transition. This was pointed out for the LFV decays in [8], and, more importantly in [29] (for a calculation of B 0 s → 1 2 γ in the model of [30]). In addition, P → 1¯ 2 decays suffer from chiral suppression (see Eq. (7)), which three-body radiative decays do not neccessarily exhibit. Thus, it is possible that RLFV decays might have larger branching ratios than two-body LFV transitions (see [24][25][26][27][28] for similar effects in lepton flavor conserving decays).
Here we evaluate radiative lepton-flavor violating decays of the pseudoscalar mesons with the model-independent effective Lagrangian of Eq. (1).
It might be theoretically easier to deal with a three-body final state that contains no strongly-interacting composite particles. Still, the calculation of the P → 1 2 γ decay is more complicated than P → 1 2 , where all nonperturbative effects are summarized in one decay constant f P . Further, because of the electromagnetic gauge invariance, it is important to have a good understanding of what kind of constraints the kinematic structure of the decay amplitude imposes on the dynamics of these transitions. Let us now derive the most general amplitude for P → 1 2 γ.
A. General amplitude and differential decay rate for P → 1 2 γ The most general expression for the P (p) → γ(k) 1 (p 1 ) 2 (p 2 ) decay amplitude can be obtained using the Bardeen-Tung formalism [31]. The decay amplitude can be written as where u(p 1 , s 1 ) and v(p 2 , s 2 ) are spinors for 1 and¯ 2 , q = 1 2 (p 1 − p 2 ), and ε * µ (k) is the polarization vector of the photon. The function M µ (p, k, q), which we seek to parameterize, transforms as a tensor under Lorentz transformations. This function should only contain dynamical singularities, so particular care should be taking by writing it in such a way that it does not contain kinematical ones. The most general expression for the M µ (p, k, q) from Eq. (9) can be written by expanding it into simpler Lorentz structures µ i (p, q, k) multiplied by the invariant functions M P 1 2 i , which only depend on Lorentz invariants, The most general parameterization of Eq. (10) contains twelve form-factors, + p µ M q 1 2 9 + / kM q 1 2 10 + iγ 5 p µ M P 1 2 11 + / kM P 1 2
In writing of Eq. (11) we used the equation of motion for the lepton spinors, and rewrote terms containing σ µν in terms of components, e.g. iσ µν q ν = q µ − γ µ / q. Note that terms proportional to / q can be expressed as terms proportional to / k using momentum conservation and equations of motion. Next, terms proportional to the µναβ tensor, such as µναβ γ ν p α k β , can be written in terms of the existing form factors of Eq. (11) using the relation and the equations of motion. Finally, all possible terms in Eq. (11) proportional to k µ trivially vanish by gauge invariance.
The set of Eq. (11) is still not minimal, as the condition of gauge invariance k µ M µ (p, k, q) = 0 implies that some of the M P 1 2 i in Eq. (11) are not independent. An elegant way of finding the minimal set of gauge-invariant Lorentz structures has been given in [31], which we shall apply to our analysis. To get the minimal set, it is most convenient to apply a projection operator to M µ (p, k, q). Since P µν M ν = M µ and k µ P µν = 0, P µν does indeed project out gaugeinvariant structures in M µ (p, k, q). Applying P µν to Eq. (11) we learn that terms proportional to p µ do not give contributions to the minimal set and should be dropped, leaving the number of independent amplitudes at eight. Applying the condition k µ µ i = 0 and eliminating kinematical singularities we write the Lorentz structures L µ i for the set of amplitudes as which are defined in a manner that removes all kinematical singularities. The A P 1 2 i (p 2 , ...) are new scalar form factors, while L µ i are This implies that the decay amplitude can be written as Using this general amplitude for a three-body pseudoscalar decay, P → γ 1 2 , we calculate a general differential decay rate, which depends on the same scalar functions Here the Mandelstam variables have the usual definitions: m 2 12 = (p 1 + p 2 ) 2 , m 2 13 = (p 1 + k) 2 , m 2 23 = (p 2 + k) 2 , where p 1,2 is the 1,2 lepton momentum, k is the γ photon momentum, and they are related to the pseudoscalar momentum, p, by p = p 1 + p 2 + k. The mass m P is the pseudo-scalar mass, m 2 is the heavier lepton mass, and y = m 2 /m P . The superscript of P 1 2 on the scalar functions A P 1 2 i (p 2 , ...) is dropped for brevity in Eq. (17). We introduce a photon mass, m γ , to regulate the infrared divergences that will appear via bremsstrahlung diagrams. We use a value of m γ = 60 MeV as our cut-off, which is near the final state invariant mass resolution of experiments [29].  decays, such as P + → γ +ν or P 0 → γ ¯ . These form factors are defined as [25][26][27]29] Here Q = p − k and the tensor form factors are defined for an off-shell photon. The tensor form factors f P T 1,2,3 [k 2 1 , k 2 2 ] are functions of two variables: k 1 , which is the momentum flowing from a vertex associated with the tensor current, and k 2 , which is the momentum of the photon emitted from the valence quark of the meson. Note that for the on-shell photon , so the tensor matrix element simplifies to [25] Using Eqs. (18), (19), and (21) we can calculate the scalar function contributions of the axial, vector, and tensor operators from the Lagrangian in Eq.
where q 1 = q 2 , which are found in Fig. (1). The contributions of these diagrams to the Note that in this section (e.g. in writing Eq. (22)) we suppressed the previously used superscript of P 1 2 in favor of a superscript related to the associated diagrams, which consists of the figure number and sub-figure letters (i.e. 1ab). We only show the odd subscript scalar function equations. The even subscript equations can be found from the odd subscript equations by replacing the left-handed WCs by their negative magnitudes (i.e. ) and multiplying the odd subscript scalar function by the imaginary constant i. This may be used to find A 2 from A 1 , A 4 from A 3 , A 6 from A 5 , and A 8 from A 7 and is true throughout this section.
There is contribution in Fig. (1) from the pseudo-scalar operators of the Lagrangian in Eq. (4). This can be seen by taking a matrix element of the divergence of axial current to relate the axial and pseudo-scalar matrix elements, and using Eq. (18) to get A similar argument can be made to prove that the scalar operators also do not give form factor contributions.
The bremsstrahlung diagrams in Fig. (2) are calculated similarly to the two-body decays of Sect. III using the matrix element of Eq. (8). We have given the photon a small mass, m γ , to regulate the infrared divergences. This divergence only appears in the quark flavor changing axial and pseudoscalar operator terms of the scalar functions, Eq. (25), so the photon mass is set to zero for the non-divergent terms. The same is true for the differential The black circles represent the four-fermion LFV vertex defined in L ef f of Eq. (4).
decay rate in Eq. (17). The axial and pseudoscalar operator scalar function terms are defined here as The dipole operator diagrams of Eq. (2) found in Fig. (3) contain contributions from the SM dipole penguin operator, Eq. (3). This is directly related to both the on and offshell tensor matrix elements in Eqs. (20) and (21) from which we need to write matrix elements of the form γ(k)|q 1 σ µν (1 ± γ 5 )q 2 |P (p) . These can be found by using the relation σ µν γ 5 = − i 2 µναβ σ αβ , which yields: The on-shell matrix element in Eq. Eq. (27). Using these matrix elements we find the dipole operator components of the scalar  2)). Note that the contributions of these diagrams are severely constrained by already available data on 1 → 2 γ decays.

functions which are
where we have used the shorthand notations f P T,I and f P T,II that we define as So far we have not addressed the contributions of the diagrams in Fig. (4). These diagrams contain contributions from the axial, vector, and tensor operators from the Lagrangian in Eq. (4) of type 1 2 qq, where the quarks are both the same flavor. As was the case for the T i +f q 2 T i (e.g. see [32]). For convenience we will use a definition with the quark charge explicitly included in the formula, f T i = Q q 1 f q 1 T i + Q q 2 f q 2 T i . This is important because in the case of Fig. 4(a) we only have contributions from f q 1 T i and in Fig. 4(b) we only have f q 2 T i .
Applying this information to the decays of B 0 q ,D 0 , and K 0 mesons shown in Figs. (1)-(4), we find that each scalar function A P 1 2 i is written as which are functions of model independent form factors and decay constants.

A. Spectra
Inputting the scalar functions of Eq. (31) in the differential decay rate, Eq. (17), and integrating over the Mandelstam variables m 2 23 and m 2 12 , we calculate the differential decay rate, dΓ/dm 2 12 , and total decay rate, Γ P → γ 1 2 . Using these results we may predict the differential decay spectra for individual operators, (1/Γ) (dΓ/dE γ ). Where we make the variable change from m 2 12 to E γ , the photon energy in the meson rest frame, and normalize to the total decay rate. This analysis requires the practical assumption of single operator dominance so that the unknown WCs of individual operators will cancel between the differential and total decay rates.
Here we have suppressed the superscripts of the WCs for brevity (e.g. C q 1 q 2 1 2 We drop terms higher in order than y 2 , which is a good approximation in most cases as the ratio y is small. The vector and tensor operators with flavor change on both the quark and lepton side are of particular importance to our analysis. They cannot be constrained via two-body decays and so the three-body decay channels present us with a unique opportunity to place limits on the associated WCs. The vector operators also have an advantage over the tensor operators because they are not chirally suppressed by quark and lepton masses. Assuming WCs are of similar size, this means the vector operators would give a larger contribution to the overall decay rate and conversely are better constrained by experimental   Table IV for B 0 q decays and in Table V forD 0 and K 0 L decays. K 0 L is used in lieu of K 0 for the limits on the branching ratios due to a lack of experimental information on the total decay rate of K 0 . The normalized differential decay plots of K 0 are the same as K 0 L because the normalization to the total decay rate cancels out the numerical differences (i.e. a factor of 1/ √ 2).
The predicted upper limits of the four-fermion axial and pseudo-scalar operators for radiative pseudo-scalar decays P → γ 1 2 in Tables IV and V demonstrate that these operators Wilson coefficient constraints using form factors for four-fermion axial and pseudo-scalar operators of type O ∼ ( 1 2 )(q 1 q 2 ) where q 1 = q 2 . Note the K 0 L results are for short distance (SD) interactions.

Wilson
Upper limits   Table I we see they are one to two orders of magnitude smaller. Therefore the three-body decays could still provide complimentary access to these operators.  [16,33,34]. In our previous work we were able to provide complimentary access via two-body vector quarkonium decays V → γ 1 2 [8].

The tensor form factors in Appendix
Using the WC constraints obtained from the radiative lepton decays 2 → 1 γ in [8], we predict the dipole operator decay upper limits for P → γ 1 2 in Table VI. Here the predicted upper limits range from 10 −21 -10 −38 , which is much lower than we would expect to be within experimental reach during the foreseeable future. Despite showing that P → γ 1 2 is not a useful means to constrain the dipole operators, the results in Table VI  (1). We argued that those Wilson coefficients can be constrained through the studies of radiative B 0 q ,D 0 , and K 0 decays to two different flavored leptons.
It is clear that studies of two-body P → 1 2 decays allowed for the quantum number selection of a smaller subset of the effective operators, which reduced our reliance on single operator dominance. Yet, the radiative three-body decays to γ 1 2 allowed access to the effective operators in Eq. (1) which cannot be probed via any two-body meson decays.
In addition to probing new operators, the three-body radiative transitions also allowed for complimentary access to four-fermion operators constrained by two-body decays without the need to include a composite strongly-interacting meson to the final state. Finally, we provide evidence that the dipole operators are so well constrained by radiative LFV transitions As more data is produced by Belle II and the LHCb experiment, we emphatically encourage our experimental colleagues to produce experimental limits on both LFV and radiative LFV decays of the B 0 q ,D 0 , and K 0 mesons discussed in this work.

Acknowledgments
We would like to thank Dmitri Melikhov  To estimate differential decay rates and the upper limits of the total decay rates of the radiative decays in Section IV, we must apply the form factors of Eqs. (18)- (20) and the numerical constants of Tables VII and VIII. Numerical inputs for the CKM matrix elements are found in [16]. Before we can apply these form factors, we must relate them to those calculated in the literature, which are defined as [25][26][27][28][29] These form factors are functions of two momenta, k 1 , which is emitted from the q 1 → q 2 weak transition current, and k 2 , which is emitted from one of the valence quarks of the meson   [25]. P . Here the photon is off-shell, but the on-shell definitions may be found by assuming k 2 2 = 0 and applying the momentum conservation relation p = k 1 + k 2 .
Assuming k 2 = 0 and making the appropriate substitutions of Q = p − k and k for k 1 and k 2 we find the necessary relations between the form factors in Eqs. (18)- (20) and Eq.
(A1) as To make use of these relations we employ the parameterizations of [25] for the F V , F A , F T V , and F T A form factors. For the B 0 q → γ form factor parameterization when the photon γ is emitted from the valence quarks (k 1 = Q, k 2 = k) we use where E γ is the photon energy in the P -meson rest-frame. The constants β and ∆ are numerical parameters which can be found in Table IX.
For the parameterization of theD 0 , K 0 → γ form factors when the photon γ is emitted  [29,36]. The K 0 tensor form factors will be calculated elsewhere. Parameter Here Q d(s) = − 1 3 , Q u(c) = 2 3 , and the remaining parameters are found in Table X [29]. The form factors F P T V , T A [0, Q 2 ] for B 0 q andD 0 decays are parameterized using vector meson dominance in [27,28], which gives The vector meson dominance input parameter values are found in Table XI Tables VII and VIII we are able to plot the normalized differential decay rates and estimate the upper limits for the radiative branching ratios assuming single operator dominance in Section IV.   Leptons Quark

Appendix B: Quark Model
When the necessary form factors are unavailable to take a model independent approach to the calculation of the four-fermion operator contributions of the diagrams in Fig. (4), we may choose a model dependent approach. We apply a constituent quark model to calculate the contributions of four-fermion vector, axial, and tensor operators of the type ( 1 2 )(qq).
We constrained both the vector and tensor Wilson coefficients for these operators previously in [8]. The results are reproduced here in Table XII and can be used to find a predicted upper bound on the branching ratio of B P → γ 1 2 for individual operators using the single operator dominance assumption.

Consituent Quark Model
The amplitude for the diagrams in Fig. (4) using this model is This amplitude is dependent on matrix elements of the form 0| q 1 Γq 2 |P with the ma- In modeling the quark anti-quark distribution, we chose to follow [38][39][40], where we can write the wave function of the ground state, P (p), as The variable x is the momentum fraction of one of the quarks and I c is the identity matrix of color space. We have assigned the momenta in Fig. (4) such that the valence quark q 1 has momentum xP and the valence quark q 2 has momentum (1 − x)P .
Here m q L is the mass of the light quark and the normalization is related to the decay constant f P . By taking the trace and integrating over the momentum fraction we find the matrix element 0| q 1 Γ µ q 2 |P = 1 0 Tr[Γ µ ψ P ] dx. (B6)