Effective field theory approach to lepton number violating decays $K^\pm\rightarrow \pi^\mp l^{\pm}l^{\pm}$: short-distance contribution

This is the first paper of our systematic efforts on lepton number violating (LNV) hadronic decays in the effective field theory approach. These decays provide information complementary to popular nuclear neutrinoless double-$\beta$ ($0\nu\beta\beta$) decay in that they can probe LNV interactions involving heavier quarks and charged leptons. We may call them hadronic $0\nu\beta\beta$ decays in short, though $\beta$ refers to all charged leptons. In this work we investigate the decays $K^\pm\rightarrow\pi^\mp l^{\pm}l^{\pm}$ that arise from short-distance or contact interactions involving four quark fields and two charged lepton fields, which have canonical dimension nine (dim-9) at leading order in low energy effective field theory (LEFT). We make a complete analysis on the basis of all dim-9 operators that violate lepton number by two units, and compute their one-loop QCD renormalization effects. We match these effective interactions in LEFT to those in chiral perturbation theory ($\chi$PT) for pseudoscalar mesons, and determine the resulting hadronic low energy constants (LECs) by chiral symmetry and lattice results in the literature. The obtained decay rate is general in that all physics at and above the electroweak scale is completely parameterized by the relevant Wilson coefficients in LEFT and hadronic LECs in $\chi$PT. Assuming the standard model effective field theory (SMEFT) is the appropriate effective field theory between some new physics scale and the electroweak scale, we match our LEFT results to SMEFT whose leading effective interactions arise from LNV dim-7 operators. This connection to SMEFT simplifies significantly the interaction structures entering in the kaon decays, and we employ the current experimental bounds to set constraints on the relevant Wilson coefficients in SMEFT.


Introduction
Neutrino oscillation experiments have confirmed that neutrinos have mass and their weak interactions with charged leptons mix lepton flavors. However, the origin of mass and the nature of neutrinos are still unclear. Being neutral, neutrinos could be Majorana particles and would thus violate lepton number conservation that arises as an accidental symmetry in the standard model (SM). This issue can be addressed both at high energy colliders where a typical signal for lepton number violation would be pair production of like-sign charged leptons from new heavy particles, and in high intensity experiments where the most popular so far is the nuclear neutrinoless double β (0νβ β ) decay appearing as a low energy manifestation of new physics at a high scale. The two types of experiments are necessary and complementary in searching for imprints of the Majorana nature of neutrinos.
Starting from this work we will make a systematic study on lepton number violating (LNV) decays of the mesons and the τ lepton. The motivations for the efforts are evident. There are rich data on the LNV decays of the charged mesons such as K ± , D ± , D ± s , B ± and the τ lepton from experiments such as LHCb, Babar, Belle, etc [1,2,3,4,5,6,7,8,9,10,11,12,13,14], and the bounds on some of the decays are expected to be considerably improved in proposed or upgraded experiments. This is particularly relevant considering the null results in current experiments on nuclear 0νβ β decays; for reviews, see, e.g., Refs. [15,16]. From the theoretical point of view the above decays involve heavier quarks and charged leptons that do not appear in nuclear 0νβ β decays, and thus can at least provide new information on lepton number violation that cannot be extracted from nuclear 0νβ β decays. Since heavy quarks and leptons may have enhanced interactions with new particles compared to light quarks and the electron, this might also be the case with LNV processes at low energy, although we are aware that the data samples in meson decays are generally much less than in nuclear 0νβ β decays. Another advantage is that one avoids uncertainties associated with nuclear physics in nuclear 0νβ β decays and that for kaon and B mesons and the τ lepton we can employ well-established chiral perturbation theory or heavy quark effective theory whose errors can be systematically estimated.
In this work we will investigate the decays K ± → π ∓ l ± l ± (l = e, µ) in the framework of effective field theory (EFT). Our approach is pictorially explained in Fig. 1 which shows the series of EFTs relevant to the decays. We start with the low energy effective field theory (LEFT) for quarks (excluding the top) and leptons that enjoys the QCD and QED gauge symmetries SU(3) C × U(1) EM [17,18]. There are then two types of contributions to the decays at the quark level: one is long-distance and the other is short-distance. We will concentrate in what follows on the short-distance or contact contribution, and reserve the long-distance part for a future publication. At the leading order in LEFT, the short-distance contribution arises from effective interactions of dim-9 LNV operators that involve four quark fields and two charged lepton fields. To calculate the meson decay rate we match at the chiral symmetry breaking scale Λ χ ≈ 4πF π the effective interactions of quarks to those of mesons that can be organized in chiral perturbation theory (χPT) [19,20]. The result thus obtained is general in the sense that it depends only on the Wilson coefficients of effective interactions in LEFT based on the QCD and QED gauge symmetries, and on the hadronic low energy constants (LECs) in χPT parameterizing nonperturbative QCD physics. These LECs may be extracted by chiral symmetry from other measured processes or computed in lattice theory. The merit in such an approach is that the uncertainties incurred in the result may be estimated systematically. This is in contrast to the studies in the literature [21,22,23,24,25,26,27] where hadronic models or approximations such as vacuum insertion are appealed to estimate the hadronic matrix elements.
To translate the experimental constraints at low energy to those on new physics at a high scale, we have to climb up the ladder of energy scales in Fig. 1. If there are no new particles with a mass at or below the electroweak scale Figure 1: Flow chart for a general EFT study is exemplified by the decay K − → π + l − l − from SMEFT, through LEFT, and to χPT, with a sequence of matching calculation and renormalization group equations (RGEs).
Λ EW ≈ m W , SM serves as a good starting point for an effective field theory, namely, the standard model effective field theory (SMEFT) [28,29,30,31,32,33,34,35,36], between some new physics scale Λ NP and the scale Λ EW . SMEFT enjoys the complete SM gauge symmetries SU(3) C × SU(2) L × U(1) Y but does not assume other symmetries such as lepton or baryon number conservation. When we match SMEFT and LEFT at Λ EW by integrating out heavy SM particles, i.e., the Higgs boson h, the weak gauge bosons W ± , Z, and the top quark t, it turns out that the interaction structures entering the above general result for the decay rate simplify significantly. This simplification would disappear generically if one assumes a different EFT above the scale Λ EW , such as the νSMEFT with relatively light sterile neutrinos [37,38,39,40] or the Higgs-Electroweak Chiral Lagrangian (EWCHL ) [41].
The outline of this paper is as follows. In Sec. 2, we first find out in LEFT a complete and independent basis for all dim-9 operators that violate lepton number L by two units (|∆ L| = 2) relevant to the type of processes under consideration. Other seemingly independent operators are removed as redundant in App. A. Then, in Sec. 3 we calculate the one-loop QCD contribution to the anomalous dimension matrix for the basis operators, and we correct as a byproduct the submatrix relevant to nuclear 0νβ β decay obtained in the literature [42]. Our analytical and numerical solutions to the complete set of RGEs are delegated to App. B. In Sec. 4 we match the above operators downwards the scale to those in χPT, and formulate a master formula for the decay width of K − → π + l − l − . And in Sec. 5 we do the opposite by matching the dim-9 operators in LEFT upwards the scale to dim-7 operators in SMEFT, and obtain the bounds on the relevant Wilson coefficients in SMEFT. Our conclusion is summarized in Sec. 6, while App. C and App. D reproduce respectively the dim-7 basis operators in SMEFT and renormalization group equations relevant to the decay under consideration.

Basis for dim-9 LNV operators in LEFT
Low energy effective field theory (LEFT) is an effective field theory for quarks and leptons that respects the QCD and QED gauge symmetries SU(3) C × U(1) EM . If SM is considered as a fundamental theory or the leading interactions of an effective field theory (SMEFT), LEFT is the effective field theory below the electroweak scale Λ EW in which spontaneous symmetry breaking has already taken place and heavy particles like the Higgs and weak gauge bosons and the top quark have been integrated out. It has a natural expansion parameter of momentum (derivative) over Λ EW , and the relative importance of effective interactions is judged by the canonical dimension of their operators.
For a consistent physics analysis it is important to establish first a basis of complete and independent operators.
The short-distance contribution to the meson decays under consideration originates at leading order in LEFT from dim-9 LNV operators that involve four quark fields and two lepton fields. The basis of dim-9 operators relevant to nuclear 0νβ β decay has been studied by several groups [42,43,44,45]. For kaon decays, Ref. [27] suggested a set of operators by assigning a free flavor index to quark and lepton fields in the operators given in Ref. [42]. As we will show below, this set is incomplete. We will make a thorough analysis on the basis of operators and verify our result by an independent approach based on the Hilbert series.
Since the operators with L = +2 are Hermitian conjugates of those with L = −2, it is sufficient to consider the latter which correspond to the decay M − 1 → M + 2 l − α l − β with M ± i being charged mesons. By Fierz identities, we can always reshuffle the fields to reach a quark-lepton separated form of operators: Here u p X (d p X ) is an up-type (a down-type) quark field of flavor p and chirality X, l α a charged lepton field of flavor α and l C α is its charge conjugate. Γ i refers to the standard sixteen Dirac matrices, while the two brackets ( , ) and [ , ] indicate the two pairs of color contraction in the products of quark fields. We denote the lepton bilinears as The bilinears in Eqs. (2) and (3) are respectively symmetric and antisymmetric under the interchange of the lepton flavors α and β . Thus the latter disappear for identical leptons, and for simplicity we will focus on the decay into identical leptons K − → π + l − l − in our later phenomenological analysis.
By making judicious applications of the properties of Dirac matrices and Fierz identities derived in Refs. [32,46], we find all dim-9, L = −2 operators shown in Tab Operator Notation  Removing the redundant operators due to flavor symmetries, there are in total 5886 independent operators with five quarks and three charged leptons. A detailed discussion on removing redundancy is delegated to App. A, in which some apparently independent operators are shown to be actually redundant and decomposed into a linear combination of the basis operators in Tab. 1. We have verified our results by the Hilbert series approach in Ref. [35] which counts the number of independent operators for a specified set of fields but does not spell out their forms.
Let us compare our results with those in the literature. First, the operators with a tensor lepton bilinear are new to Ref. [27], and moreover there are color exchanged operators that are not included in that reference. Second, if we restrict ourselves to the subspace of operators relevant to nuclear 0νβ β decay, our results match the ones in Refs. [42] and [45] when Fierz identities are applied. The relations among the three bases are where the first two notations in each line refer respectively to Refs. [42] and [45] and the last one is ours, and N = 3 is color number. Since the operators in Ref. [45] do not include the lepton bilinear, we have implicitly striped the lepton bilinears from the other two bases in order to write down the above relations. We have also dropped a factor of 4 in the definition of operators in Ref. [42], which will not affect comparison of RGE results in Sec. 3. These relations further confirm our results against Ref. [27], and will be used in Sec. 3 to clarify differences in one-loop QCD running to Ref. [42].
Our basis of operators shown in Tab. 1 is responsible for all leading-order short-distance mechanisms of low energy processes that violate lepton number by two units: where the subscript i covers all indices appearing in the operators. These processes include in particular the popular nuclear 0νβ β decays and the three-and four-body decays of the charged mesons K, D, D s , B and the τ lepton. The current experimental upper bounds on the three-body decays of the charged meson and τ lepton are summarized in Tab. 2. Using the above effective Lagrangian we can match downwards the scale to heavy quark effective theory or χPT that is appropriate to the process under consideration, and the upper bounds then translate into constraints on the Wilson coefficients C i . We can further match L |∆L|=2 LEFT upwards the scale to an EFT such as SMEFT defined between the electroweak scale Λ EW and some new physics scale Λ NP , so that we can set constraints on Λ NP . In this manner all low energy data are connected through a sequence of EFTs to potential new physics defined at a high scale.
Before finishing this section we record here the 36 operators (involving 26 four-quark combinations) that contribute where (i, j) = (d, s), (s, d) and the operators in the second line are the parity partners in the first.

QCD RGEs for dim-9 LNV operators in LEFT
To improve convergence in fixed-order perturbation theory we have to resum the large logarithms generated between two well-separated energy scales, i.e., Λ EW and Λ χ for the processes under consideration. In this section we compute the one-loop QCD anomalous dimension matrix for all the LEFT operators shown in Tab. 1 and Eq. (16). The Feynman diagrams for the quark factors (u of the basis operators are displayed in Fig. 2. We perform the calculation with dimensional regularization and in the general R ξ gauge with gauge parameter ξ 3 , and work with the MS scheme. The disappearance of ξ 3 in the final results will serve as a check of our calculation.
Since QCD conserves parity, it is sufficient to consider the operators listed in the left part of Tab. 1. Computing the Feynman diagrams in Fig. 2 and including the field wavefunction renormalization, we get the following one-loop QCD RGEs for the Wilson coefficients: where α s = g 2 3 /4π is the strong coupling constant and C F = (N 2 − 1)/2N = 4/3 the second Casimir invariant of quarks. We see that there remain no ξ 3 -dependent terms. The general solutions to the above RGEs are given in App. B, together with numerical estimates running from the electroweak scale Λ EW to the chiral symmetry breaking scale Λ χ .
In Ref. [42], the one-loop QCD RGEs were computed for the subset of operators contributing to nuclear 0νβ β decay. If we restrict ourselves to the same subset, our RGEs reduce to µ d dµ For comparison, we now recast our results in terms of the operator basis in Ref. [42]. Consider two bases of operators  (15). To summarize the result of comparison, we confirmed the anomalous dimension matricesγ XY 31 ,γ XY 12 , andγ XX 3 in the notation of Ref. [42], but found differences for the other two matrices, for which we obtained Ourγ XX 45 is half the result in Ref. [42], whileγ XY 45 is completely different. We also computed those two matrices in the basis of that reference and confirmed our result.

Matching onto chiral perturbation theory
Our discussion on LEFT in the previous two sections applies generally to the case with five quarks (and all leptons).
To study the specific hadronic process K − → π + l − l − , we restrict ourselves to the LEFT with the three light quarks u, d, s and match it at the scale Λ χ to χPT for the Nambu-Goldstone bosons. The only guide we have for this matching calculation is the spontaneous chiral symmetry breaking SU(3) L ⊗ SU(3) R → SU(3) V , whose consequences can be systematically worked out by χPT [19,20]; see Refs. [44,45,47,48] for discussions in the case of lepton number violation. We will follow closely the technique clearly demonstrated in Ref. [45].
In the matching to χPT the lepton bilinear of a dim-9 operator in LEFT behaves as a fixed external source, thus we only have to cope with the quark factor of the operator. Suppose the latter has been decomposed into a sum of irreducible representations (irreps) of the chiral group. A general irrep takes the form, where the set of pure numbers T ab cd depends on the irrep under consideration. T is promoted as a spurion field that transforms properly together with chiral transformations of quarks, where L ∈ SU(3) L and R ∈ SU(3) R , so that O looks like a chiral invariant. On the χPT side, we introduce the standard matrix for the Nambu-Goldstone bosons (NGBs), where F 0 is the decay constant in the chiral limit and ξ 2 ≡ Σ. Under chiral transformations we have where U ∈ SU(3) V depends on the NGB fields. To form leading-order (LO) operators, matching is accomplished by the substitutions, where the free indices are to be contracted when forming an operator with T ab cd . For hadronic operators appearing at the next-to-leading (NLO) or next-to-next-to-leading order (NNLO), covariant derivatives and quark masses will be involved in matching: where the quark mass M = diag(m u , m d , m s ) as a spurion transforms like M → R † ML under the chiral group, and covariant derivatives transform as D µ ξ → UD µ ξ R † , D µ ξ † → UD µ ξ † L † , with D µ = ∂ µ + (ξ † ∂ µ ξ + ξ ∂ µ ξ † )/2. They are related to the ordinary derivative of the Σ field by ξ D µ ξ † = (Σ∂ µ Σ † )/2 and ξ † (D µ ξ ) = (Σ † ∂ µ Σ)/2. The identities (D µ ξ )ξ † = −ξ (D µ ξ ) † and ξ † (D µ ξ ) = −(D µ ξ ) † ξ are useful when reducing redundant operators. 4 Now we turn to the decay K − → π + l − l − . According to the procedure outlined above, we first decompose the quark factors in operators of Eq. (17) into irreps under the chiral group, and then we obtain their leading nonvanishing hadronic counterparts according to Eqs. (36) and (38). Our results are shown in Tab. 3. The matching coefficients are denoted by the LECs, g X . These constants are difficult to compute because of strong dynamics, but some of them may be related by chiral symmetry to other constants that have been experimentally measured or computed by lattice simulations. In particular, the matrix elements of K − → π + associated with g 27×1 , g i 6×6 , and g i 8×8 can be related to those of π − → π + which are responsible for the short-range mechanism of nuclear 0νβ β decay and computed in [48], and to those of K + → π + π 0 [47] and K 0 →K 0 [48] which have been computed by lattice simulations in [49,50,51,52,53] and [54,55], respectively. Here we adopt the results from [56]; in our notations: We make some clarifications concerning the results in Tab. 3. An operator with a symbol (P) is the parity partner of a corresponding operator with a ( ). Parity invariance of QCD implies that such a pair of operators shares the same LEC: This leaves us with 13 LECs. We are now ready to write down the effective Lagrangian contributing to the decay K − → π + l − l − according to Eqs. (16)-(17) and Tab. 3 at the leading order of each operator: where the parameters c i are  Table 3: Quark factors of dim-9 operators in LEFT are matched at the leading nonvanishing order to their hadronic counterparts in χPT that are relevant to K − → π + l − l − . Nonvanishing lepton bilinears j = (ll C ), j 5 = (lγ 5 l C ), j µ5 = (lγ µ γ 5 l C ) for an identical lepton pair l = e, µ are not shown in hadronic factors for brevity.
The spin-summed squared matrix element for the decay where m K,π,l are the masses of the K − , π + , and l respectively, and s = (q 1 + q 2 ) 2 , t = (p + q 2 ) 2 . The decay width is calculated as where the integration domains are with We have been working so far at the first nonvanishing order of each operator in χPT. To estimate the errors incurred by ignoring higher-order terms we calculate the one-loop chiral logarithms contributing to the K − → π + transition amplitude as shown by Feynman diagrams in Fig. 3. These nonanalytic chiral logarithms cannot be cancelled by the counterterms to the one-loop diagrams, and thus can be used as a rough estimate of higher order terms. For this, we need the K − π − |P| 2 terms from expansion of the hadronic operators shown in Tab. 3 with P being any NGB, and the where B 0 is the order parameter for chiral symmetry breaking through the quark condensate, −3F 2 0 B 0 = 0|qq|0 . We will use the Gell-Mann-Okubo formula for NGBs, 3m 2 η = 4m 2 K − m 2 π . For the purpose of illustrating the relative size of chiral logarithms to the leading terms, we consider the interactions from g 27×1 , g a/b 8×8 , and g a/b 6×6 in Tab. 3. As will be clear in the next section, only the first two types of couplings can be obtained from matching to the SMEFT dim-7 operators. Since this discussion concerns only the accuracy of χPT at LO, we consider only the hadronic factors of those operators with a ( ) in Tab. 3, instead of a complete effective interaction, which is a sum involving various unknown Wilson coefficients and different lepton bilinears. Dropping the momentum squared of the leptonic system, (k − p) 2 = 0, we find where L P = m 2 P /(4πF 0 ) 2 ln(µ 2 /m 2 P ) and µ is the renormalization scale. We have taken into account in the same approximation the renormalization of the decay constants [20] and the wave function renormalization constants The relative corrections in the three amplitudes have a magnitude of about 49.8%, 41.6%, 28.3% (34%, 32.5%, 31.3%) respectively, at the renormalization scale µ = Λ χ (µ = m K ). This roughly fits the usual expectation on the accuracy of SU(3) χPT.

Matching to SMEFT
Our previous results on short-distance contributions to the decay K − → π + l − l − are general in that they are not specific to physics above the electroweak scale Λ EW but rely only on well established symmetries at low energy. But to study the impact of low energy measurements on new physics at a high scale, we should connect our results to the effective field theory above Λ EW . For this we make the minimal assumption that there are no new particles of mass of order Λ EW or smaller, so that SM appears as the leading terms in an EFT, i.e., SMEFT. SMEFT is defined between Λ EW and some new physics scale Λ NP where new heavy particles have been integrated out. Its Lagrangian consists of a tower of effective operators built by the SM fields and satisfying the SM gauge symmetry SU(3) C × SU(2) L ×U(1) Y [28,29,30,31,32,33,34,35,36]: where O d i is the i-th operator of dimension d with the corresponding Wilson coefficient C d i . Lepton number violation first appears at dim-5 through the Weinberg operator that yields the Majorana neutrino mass upon electroweak symmetry breaking [28]. This operator however is not directly relevant to our purpose of calculating short-distance contributions to the LNV K ± decay. The leading contributions then arise from dim-7 operators in SMEFT. The effects of dim-7 operators in nuclear 0νβ β decay have been analyzed recently by many groups, see, e.g., Refs. [56,57,33] and references cited therein, and a few Wilson coefficients have been constrained from the experimental bounds [58,59].
The Wilson coefficients to be constrained by the K ± decay below are actually left free in nuclear 0νβ β decay.
where both sides are evaluated at the scale µ = Λ EW , and G F is the Fermi constant, V ud and V us are the CKM matrix elements. On the right-hand side we have used the basis and notation of dim-7 operators in Ref. [32]; e.g., C 21ll † duLLD refers to the s, u quarks and the lepton l = e, µ.
The above matching between LEFT and SMEFT simplifies considerably our subsequent phenomenological analysis. First of all, the QCD RGEs for the above coefficients decouple from other coefficients: where (i, j) = (d, s), (s, d). The solutions are C LLLL,S/P udus where b = −11 + 2n f /3 with n f being the number of active quark flavors between scales µ 1 and µ 2 ; in particular, C LLLL,S/P udus The coefficients c i in Eqs. (42)-(47) for the K decay simplify to and the squared matrix element becomes a compact form Normalizing the LNV K decay width to its total width [60] yields the branching ratios for the decays to two identical leptons: To get some feel about the bound on the relevant energy scale, we assume naively that the above Wilson coefficients in SMEFT scale as Λ −3 , then the experiments upper bounds in Tab. 3 translate into a loose bound Λ > O(10 GeV). This bound is indeed much weaker than that from nuclear 0νβ β decay, Λ > O(10 TeV) [33], but it concerns the quarks and leptons of the second generation. This relative weakness arises largely from much smaller data samples in the K decays.

Conclusion
The Majorana nature of neutrinos has so far been intensively explored in nuclear 0νβ β decay both experimentally and theoretically. Considering the null result in current experimental searches it is important that we seek other potential signals that could reveal the Majorana nature of neutrinos. The LNV decays of the charged mesons and τ lepton may play a role in probing interactions of heavier quarks and leptons to which nuclear 0νβ β decay is not sensitive.
In this work we have for the first time investigated the LNV decay K − → π + l − l − in the framework of effective field theory. We established the basis of |∆L| = 2 dim-9 operators in LEFT that are responsible for the leading order shortdistance contributions to LNV processes including, e.g., the decays of the mesons and τ lepton and nuclear 0νβ β First of all, from the properties of Dirac matrices, we have the following Fierz identities i 2 ε µνρσ σ ρσ = σ µν γ 5 (ε 0123 = +1), (A.1) where j αβ 5µν = l α γ 5 σ µν l β ,C and P ± = (1 ± γ 5 )/2. Thus we can discard the lepton bilinear j αβ 5µν together with the tensor ε µνρσ for contraction between the quark and lepton factors and include only lepton bilinears listed in Tab. 1.

Appendix B Solutions to RGEs for dim-9 LNV operators in LEFT
In this Appendix we solve the complete set of one-loop QCD RGEs (18)- (25). Our results in previous sections on nuclear 0νβ β decay and the decays K ± → π ∓ → l α l β form subsets of the results recorded below.
We denote the diagonal matrices formed with the eigenvalues of anomalous dimension matrices in Eqs. where ζ 2/1 = α s (µ 2 )/α s (µ 1 ) and b = −11 + 2n f /3 with n f being the number of active quarks between the scales µ 1 and µ 2 . The corresponding diagonalization matrices are found to be, has been worked out for the subset of operators violating both baryon and lepton numbers in Ref. [32] and for the subset violating only lepton number in Ref. [33]. where α i = g 2 i /(4π) with g i being the gauge couplings for the SM gauge group SU(3) C × SU(2) L × U(1) Y , α t = y 2 t /(4π) with y t being the Yukawa coupling of the top quark (in the convention, m t = y t v/ √ 2 with v ≈ 246 GeV), and α λ = λ /(4π) with λ being the Higgs self-coupling (in the convention, m 2 h = 2λ v 2 ). The above equations ignore much smaller Yukawa couplings of other fermions.