Already defined observables in FlavorKit

Many observables are already implemented in FlavorKit and can be used out-of-the-box

Lepton flavor observables

Lepton flavor violation in the SM or MSSM without neutrino masses vanishes exactly. Even adding Dirac neutrino masses to the SM predicts LFV rates which are far beyond the experimental reach. However, many extensions of the SM can introduce new sources for LFV of a size which is testable nowadays. The best-known examples are SUSY and non-SUSY models with high- or low-scale seesaw mechanism, models with vector-like leptons and SUSY models with R-parity violation, see for instance Refs. .

We discuss in the following the implementation of the most important LFV observables in SARAH and SPheno using the previously defined operators which are calculated by SPheno.

Radiative decays of leptons (ℓ_α → ℓ_βγ)

The decay width is given by 

\Gamma \left( \ell_\alpha \to \ell_\beta \gamma \right) = \frac{\alpha m_{\ell_\alpha}^5}{4} \left( |K_2^L|^2 + |K_2^R|^2 \right) \, ,

where α is the fine structure constant and the dipole Wilson coefficients K₂^L, R are defined in Eq..

NameProcess = "LLpGamma";
NameObservables = {{muEgamma, 701, "BR(mu->e gamma)"},
                   {tauEgamma, 702, "BR(tau->e gamma)"},
                   {tauMuGamma, 703, "BR(tau->mu gamma)"}};

NeededOperators = {K2L, K2R};

Body = "LLpGamma.f90";

Real(dp) :: width
Integer :: i1, gt1, gt2

! ----------------------------------------------------------------
! l -> l' gamma
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on J. Hisano et al, PRD 53 (1996) 2442 [hep-ph/9510309]
! ----------------------------------------------------------------

Do i1=1,3

If (i1.eq.1) Then         ! mu -> e gamma
 gt1 = 2
 gt2 = 1
Elseif (i1.eq.2) Then     !tau -> e gamma
 gt1 = 3
 gt2 = 1
Else                      ! tau -> mu gamma
 gt1 = 3
 gt2 = 2
End if

width=0.25_dp*mf_l(gt1)**5*(Abs(K2L(gt1,gt2))**2 &
           & +Abs(K2R(gt1,gt2))**2)*Alpha

If (i1.eq.1) Then
muEgamma = width/(width+GammaMu)
Elseif (i1.eq.2) Then
tauEgamma = width/(width+GammaTau)
Else
tauMuGamma = width/(width+GammaTau)
End if

End do

Decays into three leptons (ℓ_α → 3****ℓ_β)

The decay width is given by

\Gamma \left( \ell_\alpha \to 3 \ell_\beta \right) = \frac{m_{\ell_\alpha}^5}{512 \pi^3} \left[ e^4 \, \left( \left| K_2^L \right|^2 + \left| K_2^R \right|^2 \right) \left( \frac{16}{3} \log{\frac{m_{\ell_\alpha}}{m_{\ell_\beta}}} - \frac{22}{3} \right) \right. \\ + \frac{1}{24} \left( \left| A_{LL}^S \right|^2 + \left| A_{RR}^S \right|^2 \right) + \frac{1}{12} \left( \left| A_{LR}^S \right|^2 + \left| A_{RL}^S \right|^2 \right) \nonumber \\ + \frac{2}{3} \left( \left| \hat A_{LL}^V \right|^2 + \left| \hat A_{RR}^V \right|^2 \right) + \frac{1}{3} \left( \left| \hat A_{LR}^V \right|^2 + \left| \hat A_{RL}^V \right|^2 \right) + 6 \left( \left| A_{LL}^T \right|^2 + \left| A_{RT}^T \right|^2 \right) \nonumber \\ + \frac{e^2}{3} \left( K_2^L A_{RL}^{S \ast} + K_2^R A_{LR}^{S \ast} + c.c. \right) - \frac{2 e^2}{3} \left( K_2^L \hat A_{RL}^{V \ast} + K_2^R \hat A_{LR}^{V \ast} + c.c. \right) \nonumber \\ - \frac{4 e^2}{3} \left( K_2^L \hat A_{RR}^{V \ast} + K_2^R \hat A_{LL}^{V \ast} + c.c. \right) \nonumber \\ - \left. \frac{1}{2} \left( A_{LL}^S A_{LL}^{T \ast} + A_{RR}^S A_{RR}^{T \ast} + c.c. \right) - \frac{1}{6} \left( A_{LR}^S \hat A_{LR}^{V \ast} + A_{RL}^S \hat A_{RL}^{V \ast} + c.c. \right) \right] \nonumber \, .

Here we have defined

Â_X**Y^V = A_X**Y^V + e²K₁^X  (X,Y=L,R) .

The mass of the leptons in the final state has been neglected in this formula, with the exception of the dipole terms K₂^L, R, where an infrared divergence would otherwise occur due to the massless photon propagator. Eq. is in agreement with , but also includes the coefficients A_L**R^S and A_R**L^S.

NameProcess = "Lto3Lp";
NameObservables = {{BRmuTo3e, 901, "BR(mu->3e)"},
                   {BRtauTo3e, 902, "BR(tau->3e)"},
                   {BRtauTo3mu, 903, "BR(tau->3mu)"}
                  };

ExternalStates =  {Electron};
NeededOperators = {K1L, K1R, K2L, K2R,
 O4lSLL, O4lSRR, O4lSRL, O4lSLR ,
 O4lVRR, O4lVLL, O4lVRL, O4lVLR ,
 O4lTLL, O4lTRR };

Body = "Lto3Lp.f90";

Complex(dp) :: cK1L, cK1R, cK2L, cK2R
Complex(dp) :: CSLL, CSRR, CSLR, CSRL, CVLL, &
                    & CVRR, CVLR, CVRL, CTLL, CTRR
Real(dp) :: BRdipole, BRscalar, BRvector, BRtensor
Real(dp) :: BRmix1, BRmix2, BRmix3, BRmix4, GammaLFV
Real(dp) :: e2, e4
Integer :: i1, gt1, gt2, gt3, gt4

! ----------------------------------------------------------------
! l -> 3 l'
! Observable implemented by W. Porod, F. Staub and A. Vicente
! ----------------------------------------------------------------

e2 = (4._dp*Pi*Alpha_MZ)
e4 = e2**2

Do i1=1,3

If (i1.eq.1) Then
 gt1 = 2
 gt2 = 1
Elseif (i1.eq.2) Then
 gt1 = 3
 gt2 = 1
Else
 gt1 = 3
 gt2 = 2
End if
gt3 = gt2
gt4 = gt2

cK1L = K1L(gt1,gt2)
cK1R = K1R(gt1,gt2)

cK2L = K2L(gt1,gt2)
cK2R = K2R(gt1,gt2)

CSLL = O4lSLL(gt1,gt2,gt3,gt4)
CSRR = O4lSRR(gt1,gt2,gt3,gt4)
CSLR = O4lSLR(gt1,gt2,gt3,gt4)
CSRL = O4lSRL(gt1,gt2,gt3,gt4)

CVLL = O4lVLL(gt1,gt2,gt3,gt4)
CVRR = O4lVRR(gt1,gt2,gt3,gt4)
CVLR = O4lVLR(gt1,gt2,gt3,gt4)
CVRL = O4lVRL(gt1,gt2,gt3,gt4)

CVLL = CVLL + e2*cK1L
CVRR = CVRR + e2*cK1R
CVLR = CVLR + e2*cK1L
CVRL = CVRL + e2*cK1R

CTLL = O4lTLL(gt1,gt2,gt3,gt4)
CTRR = O4lTRR(gt1,gt2,gt3,gt4)

! Photonic dipole contributions
BRdipole = (Abs(cK2L)**2+Abs(cK2R)**2)&
&*(16._dp*Log(mf_l(gt1)/mf_l(gt2))-22._dp)/3._dp

! Scalar contributions
BRscalar = (Abs(CSLL)**2+Abs(CSRR)**2)/24._dp&
&+(Abs(CSLR)**2+Abs(CSRL)**2)/12._dp

! Vector contributions
BRvector = 2._dp*(Abs(CVLL)**2+Abs(CVRR)**2)/3._dp&
&+(Abs(CVLR)**2+Abs(CVRL)**2)/3._dp

! Tensor contributions
BRtensor = 6._dp*(Abs(CTLL)**2+Abs(CTRR)**2)

! Mix: dipole x scalar
BRmix1 = 2._dp/3._dp*Real(cK2L*Conjg(CSRL) + cK2R*Conjg(CSLR),dp)

! Mix: dipole x vector
BRmix2 = -4._dp/3._dp*Real(cK2L*Conjg(CVRL) + cK2R*Conjg(CVLR),dp) &
     & -8._dp/3._dp*Real(cK2L*Conjg(CVRR) + cK2R*Conjg(CVLL),dp)

! Mix: scalar x vector
BRmix3 = -1._dp/3._dp*Real(CSLR*Conjg(CVLR) + CSRL*Conjg(CVRL),dp)

! Mix: scalar x tensor
BRmix4 = -1._dp*Real(CSLL*Conjg(CTLL) + CSRR*Conjg(CTRR),dp)

GammaLFV = oo512pi3*mf_l(gt1)**5* &
     & (e4*BRdipole + BRscalar + BRvector + BRtensor &
     & + e2*BRmix1 + e2*BRmix2 + BRmix3 + BRmix4)

!----------------------------------------------------------------------
!taking alpha(Q=0) instead of alpha(m_Z) as this contains most of the
!running of the Wilson coefficients
!----------------------------------------------------------------------

If (i1.Eq.1) Then
 BRmuTo3e=GammaLFV/GammaMu
Else If (i1.Eq.2) Then
 BRtauTo3e=GammaLFV/GammaTau
Else
 BRtauTo3mu=GammaLFV/GammaTau
End If
End do

Coherent muon-electron conversion in nuclei

The conversion rate, relative to the the muon capture rate, can be expressed as

{\rm CR} (\mu- e, {\rm Nucleus}) = \frac{p_e \, E_e \, m_\mu^3 \, G_F^2 \, \alpha^3 \, Z_{\rm eff}^4 \, F_p^2}{8 \, \pi^2 \, Z} \nonumber \\ \times \left\{ \left| (Z + N) \left( g_{LV}^{(0)} + g_{LS}^{(0)} \right) + (Z - N) \left( g_{LV}^{(1)} + g_{LS}^{(1)} \right) \right|^2 + \right. \nonumber \\ \ \ \ \ \left. \,\, \left| (Z + N) \left( g_{RV}^{(0)} + g_{RS}^{(0)} \right) + (Z - N) \left( g_{RV}^{(1)} + g_{RS}^{(1)} \right) \right|^2 \right\} \frac{1}{\Gamma_{\rm capt}}\,.

Z and N are the number of protons and neutrons in the nucleus and Z_{\rm eff} is the effective atomic charge . Similarly, G_F is the Fermi constant, F_p is the nuclear matrix element and \Gamma_{\rm capt} represents the total muon capture rate. α is the fine structure constant, p_e and E_e ( ≃m_μ in the numerical evaluation) are the momentum and energy of the electron and m_μ is the muon mass. In the above, g_X**K⁽⁰⁾ and g_X**K⁽¹⁾ (with X = L, R and K = S, V) can be written in terms of effective couplings at the quark level as

g_{XK}^{(0)} = \frac{1}{2} \sum_{q = u,d,s} \left( g_{XK(q)} G_K^{(q,p)} + g_{XK(q)} G_K^{(q,n)} \right)\,, \nonumber \\ g_{XK}^{(1)} = \frac{1}{2} \sum_{q = u,d,s} \left( g_{XK(q)} G_K^{(q,p)} - g_{XK(q)} G_K^{(q,n)} \right)\,.

For coherent μ − e conversion in nuclei, only scalar (S) and vector (V) couplings contribute. Furthermore, sizable contributions are expected only from the u, d, s quark flavors. The numerical values of the relevant G_K factors are 

G_V^{(u, p)}\, =\, G_V^{(d, n)\,} =\, 2 \,;\, \ \ \ \ G_V^{(d, p)}\, =\, G_V^{(u, n)}\, = 1\,; \nonumber \\ G_S^{(u, p)}\, =\, G_S^{(d, n)}\, =\, 5.1\,;\, \ \ G_S^{(d, p)}\, =\, G_S^{(u, n)}\, = \,4.3 \,;\, \nonumber \\ G_S^{(s, p)}\,=\, G_S^{(s, n)}\, = \,2.5\,.

Finally, the g_X**K(q) coefficients can be written in terms of the Wilson coefficients in Eqs., and as

g_{LV(q)} = \frac{\sqrt{2}}{G_F} \left[ e^2 Q_q \left( K_1^L - K_2^R \right)- \frac{1}{2} \left( C_{\ell\ell qq}^{VLL} + C_{\ell\ell qq}^{VLR} \right) \right] \\ g_{RV(q)} = \left. g_{LV(q)} \right|_{L \to R} \\ g_{LS(q)} = - \frac{\sqrt{2}}{G_F} \frac{1}{2} \left( C_{\ell\ell qq}^{SLL} + C_{\ell\ell qq}^{SLR} \right) \\ g_{RS(q)} = \left. g_{LS(q)} \right|_{L \to R} \, .

Here Q_q is the quark electric charge (Q_d = −1/3, Q_u = 2/3) and C_ℓℓq**q^IXK = B_X**Y^K (C_X**Y^K) for d-quarks (u-quarks), with X = L, R and K = S, V.

NameProcess = "MuEconversion";
NameObservables = {{CRmuEAl, 800, "CR(mu-e, Al)"},
                   {CRmuETi, 801, "CR(mu-e, Ti)"},
                   {CRmuESr, 802, "CR(mu-e, Sr)"},
                   {CRmuESb, 803, "CR(mu-e, Sb)"},
                   {CRmuEAu, 804, "CR(mu-e, Au)"},
                   {CRmuEPb, 805, "CR(mu-e, Pb)"}
                  };

NeededOperators = {K1L, K1R, K2L, K2R,
  OllddSLL, OllddSRR, OllddSRL, OllddSLR, OllddVRR, OllddVLL,
  OllddVRL, OllddVLR, OllddTLL, OllddTLR, OllddTRL, OllddTRR,
  OlluuSLL, OlluuSRR, OlluuSRL, OlluuSLR, OlluuVRR, OlluuVLL,
  OlluuVRL, OlluuVLR, OlluuTLL, OlluuTLR, OlluuTRL, OlluuTRR
};

Body = "MuEconversion.f90";

Complex(dp) :: gPLV(3), gPRV(3)
Complex(dp),Parameter :: mat0(3,3)=0._dp
Real(dp) :: Znuc,Nnuc, Zeff, Fp, GammaCapt, GSp(3), GSn(3), &
     & GVp(3), GVn(3), e2
Complex(dp) :: Lcont,Rcont,gLS(3),gRS(3),gLV(3),gRV(3),g0LS,g0RS, &
     & g0LV,g0RV,g1LS,g1RS,g1LV,g1RV
Integer :: i1, i2

! ----------------------------------------------------------------
! Coherent mu-e conversion in nuclei
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on Y. Kuno, Y. Okada, Rev. Mod. Phys. 73 (2001) 151 [hep-ph/9909265]
! and E. Arganda et al, JHEP 0710 (2007) 104 [arXiv:0707.2955]
! ----------------------------------------------------------------

e2 = 4._dp*Pi*Alpha_MZ

! 1: uu
! 2: dd
! 3: ss

! vector couplings

gLV(1) = 0.5_dp*(OlluuVLL(2,1,1,1) + OlluuVLR(2,1,1,1))
gRV(1) = 0.5_dp*(OlluuVRL(2,1,1,1) + OlluuVRR(2,1,1,1))
gLV(2) = 0.5_dp*(OllddVLL(2,1,1,1) + OllddVLR(2,1,1,1))
gRV(2) = 0.5_dp*(OllddVRL(2,1,1,1) + OllddVRR(2,1,1,1))
gLV(3) = 0.5_dp*(OllddVLL(2,1,2,2) + OllddVLR(2,1,2,2))
gRV(3) = 0.5_dp*(OllddVRL(2,1,2,2) + OllddVRR(2,1,2,2))

gLV = -gLV*Sqrt(2._dp)/G_F
gRV = -gRV*Sqrt(2._dp)/G_F

gPLV(1) = (K1L(2,1)-K2R(2,1))*(2._dp/3._dp)
gPRV(1) = (K1R(2,1)-K2L(2,1))*(2._dp/3._dp)
gPLV(2) = (K1L(2,1)-K2R(2,1))*(-1._dp/3._dp)
gPRV(2) = (K1R(2,1)-K2L(2,1))*(-1._dp/3._dp)
gPLV(3) = (K1L(2,1)-K2R(2,1))*(-1._dp/3._dp)
gPRV(3) = (K1R(2,1)-K2L(2,1))*(-1._dp/3._dp)
gPLV = gPLV*Sqrt(2._dp)/G_F*e2
gPRV = gPRV*Sqrt(2._dp)/G_F*e2

gLV=gPLV+gLV
gRV=gPRV+gRV


! scalar couplings

gLS(1) = 0.5_dp*(OlluuSLL(2,1,1,1)+OlluuSLR(2,1,1,1))
gRS(1) = 0.5_dp*(OlluuSRL(2,1,1,1)+OlluuSRR(2,1,1,1))
gLS(2) = 0.5_dp*(OllddSLL(2,1,1,1)+OllddSLR(2,1,1,1))
gRS(2) = 0.5_dp*(OllddSRL(2,1,1,1)+OllddSRR(2,1,1,1))
gLS(3) = 0.5_dp*(OllddSLL(2,1,2,2)+OllddSLR(2,1,2,2))
gRS(3) = 0.5_dp*(OllddSRL(2,1,2,2)+OllddSRR(2,1,2,2))

gLS = -gLS*Sqrt(2._dp)/G_F
gRS = -gRS*Sqrt(2._dp)/G_F


Do i1=1,6
 If(i1.eq.1) Then
Znuc=13._dp
Nnuc=14._dp
Zeff=11.5_dp
Fp=0.64_dp
GammaCapt=4.64079e-19_dp
Else If(i1.eq.2) Then
Znuc=22._dp
Nnuc=26._dp
Zeff=17.6_dp
Fp=0.54_dp
GammaCapt=1.70422e-18_dp
Else If(i1.eq.3) Then
Znuc=38._dp
Nnuc=42._dp
Zeff=25.0_dp
Fp=0.39_dp
GammaCapt=4.61842e-18_dp
Else If(i1.eq.4) Then
Znuc=51._dp
Nnuc=70._dp
Zeff=29.0_dp
Fp=0.32_dp
GammaCapt=6.71711e-18_dp
Else If(i1.eq.5) Then
Znuc=79._dp
Nnuc=118._dp
Zeff=33.5_dp
Fp=0.16_dp
GammaCapt=8.59868e-18_dp
Else If(i1.eq.6) Then
Znuc=82._dp
Nnuc=125._dp
Zeff=34.0_dp
Fp=0.15_dp
GammaCapt=8.84868e-18_dp
End If

! numerical values
! based on Y. Kuno, Y. Okada, Rev. Mod. Phys. 73 (2001) 151 [hep-ph/9909265]
! and T. S. Kosmas et al, PLB 511 (2001) 203 [hep-ph/0102101]
GSp=(/5.1,4.3,2.5/)
GSn=(/4.3,5.1,2.5/)
GVp=(/2.0,1.0,0.0/)
GVn=(/1.0,2.0,0.0/)

g0LS=0._dp
g0RS=0._dp
g0LV=0._dp
g0RV=0._dp
g1LS=0._dp
g1RS=0._dp
g1LV=0._dp
g1RV=0._dp
Do i2=1,3
g0LS=g0LS+0.5_dp*gLS(i2)*(GSp(i2)+GSn(i2))
g0RS=g0RS+0.5_dp*gRS(i2)*(GSp(i2)+GSn(i2))
g0LV=g0LV+0.5_dp*gLV(i2)*(GVp(i2)+GVn(i2))
g0RV=g0RV+0.5_dp*gRV(i2)*(GVp(i2)+GVn(i2))
g1LS=g1LS+0.5_dp*gLS(i2)*(GSp(i2)-GSn(i2))
g1RS=g1RS+0.5_dp*gRS(i2)*(GSp(i2)-GSn(i2))
g1LV=g1LV+0.5_dp*gLV(i2)*(GVp(i2)-GVn(i2))
g1RV=g1RV+0.5_dp*gRV(i2)*(GVp(i2)-GVn(i2))
End Do
Lcont=(Znuc+Nnuc)*(g0LV+g0LS)+(Znuc-Nnuc)*(g1LV-g1LS)
Rcont=(Znuc+Nnuc)*(g0RV+g0RS)+(Znuc-Nnuc)*(g1RV-g1RS)

! Conversion rate
If (i1.eq.1) Then
 CRMuEAl =oo8pi2*mf_l(2)**5*G_F**2*Alpha**3*Zeff**4*Fp**2/Znuc*&
   & (Abs(Lcont)**2+Abs(Rcont)**2)/GammaCapt
Else if (i1.eq.2) Then
 CRMuETi =oo8pi2*mf_l(2)**5*G_F**2*Alpha**3*Zeff**4*Fp**2/Znuc*&
   & (Abs(Lcont)**2+Abs(Rcont)**2)/GammaCapt
Else if (i1.eq.3) Then
 CRMuESr =oo8pi2*mf_l(2)**5*G_F**2*Alpha**3*Zeff**4*Fp**2/Znuc*&
   & (Abs(Lcont)**2+Abs(Rcont)**2)/GammaCapt
Else if (i1.eq.4) Then
 CRMuESb =oo8pi2*mf_l(2)**5*G_F**2*Alpha**3*Zeff**4*Fp**2/Znuc*&
   & (Abs(Lcont)**2+Abs(Rcont)**2)/GammaCapt
Else if (i1.eq.5) Then
 CRMuEAu =oo8pi2*mf_l(2)**5*G_F**2*Alpha**3*Zeff**4*Fp**2/Znuc*&
   & (Abs(Lcont)**2+Abs(Rcont)**2)/GammaCapt
Else if (i1.eq.6) Then
 CRMuEPb =oo8pi2*mf_l(2)**5*G_F**2*Alpha**3*Zeff**4*Fp**2/Znuc*&
   & (Abs(Lcont)**2+Abs(Rcont)**2)/GammaCapt
End if
End do

Tau decays into mesons (τ → P****ℓ)

Our analytical expressions for τ → Pℓ, where ℓ = e, μ and P is a pseudoscalar meson, generalize the results in . The decay width is given by

\Gamma \left( \tau \to \ell P \right) = \frac{1}{4 \pi} \frac{\lambda^{1/2}(m_\tau^2,m_\ell^2,m_P^2)}{m_\tau^2} \frac{1}{2} \sum_{i,f} |\mathcal{M}_{\tau \ell P}|^2 \, ,

where the averaged squared amplitude can be written as

\frac{1}{2} \sum_{i,f} |\mathcal{M}_{\tau \ell P}|^2 = \frac{1}{4 m_\tau} \sum_{I,J = S,V} \left[ 2 m_\tau m_\ell \left( a_P^I a_P^{J \, \ast} - b_P^I b_P^{J \, \ast} \right) + (m_\tau^2 + m_\ell^2 - m_P^2) \left( a_P^I a_P^{J \, \ast} + b_P^I b_P^{J \, \ast} \right) \right] \, .

The coefficients a_P^S, V and b_P^S, V can be expressed in terms of the Wilson coefficients in Eqs. and as

a_P^S = \frac{1}{2} f_\pi \, \sum_{X = L,R} \left[ \frac{D_X^d(P)}{m_d} \left( B^S_{LX} + B^S_{RX} \right) + \frac{D_X^u(P)}{m_u} \left( C^S_{LX} + C^S_{RX} \right) \right] \\ b_P^S = \frac{1}{2} f_\pi \, \sum_{X = L,R} \left[ \frac{D_X^d(P)}{m_d} \left( B^S_{RX} - B^S_{LX} \right) + \frac{D_X^u(P)}{m_u} \left( C^S_{RX} - C^S_{LX} \right) \right] \\ a_P^V = \frac{1}{4} f_\pi \, C(P) (m_\tau - m_\ell) \left[ - B_{LL}^V + B_{LR}^V - B_{RL}^V + B_{RR}^V \right. \nonumber \\ \left. + C_{LL}^V - C_{LR}^V + C_{RL}^V - C_{RR}^V \right] \\ b_P^V = \frac{1}{4} f_\pi \, C(P) (m_\tau + m_\ell) \left[ - B_{LL}^V + B_{LR}^V + B_{RL}^V - B_{RR}^V \right. \nonumber \\ \left. + C_{LL}^V - C_{LR}^V - C_{RL}^V + C_{RR}^V \right] \, .

In these expressions m_d and m_u are the down- and up-quark masses, respectively, f_π is the pion decay constant and the coefficients C(P),D_L, R^d, u(P) take different forms for each pseudoscalar meson P . For P = π one has

C(\pi) = 1 \\ D_L^d(\pi) = - \frac{m_\pi^2}{4} \\ D_L^u(\pi) = \frac{m_\pi^2}{4} \, ,

for P = η

C(\eta) = \frac{1}{\sqrt{6}} \left( \sin \theta_\eta + \sqrt{2} \cos \theta_\eta \right) \\ D_L^d(\eta) = \frac{1}{4 \sqrt{3}} \left[ (3 m_\pi^2 - 4 m_K^2) \cos \theta_\eta - 2 \sqrt{2} m_K^2 \sin \theta_\eta \right] \\ D_L^u(\eta) = \frac{1}{4 \sqrt{3}} m_\pi^2 \left( \cos \theta_\eta - \sqrt{2} \sin \theta_\eta \right) \, ,

and for P = η^′

C(\eta^\prime) = \frac{1}{\sqrt{6}} \left( \sqrt{2} \sin \theta_\eta - \cos \theta_\eta \right) \\ D_L^d(\eta^\prime) = \frac{1}{4 \sqrt{3}} \left[ (3 m_\pi^2 - 4 m_K^2) \sin \theta_\eta + 2 \sqrt{2} m_K^2 \cos \theta_\eta \right] \\ D_L^u(\eta^\prime) = \frac{1}{4 \sqrt{3}} m_\pi^2 \left( \sin \theta_\eta + \sqrt{2} \cos \theta_\eta \right) \, .

Here m_π and m_K are the masses of the neutral pion and Kaon, respectively, and θ_η is the η − η^′ mixing angle. In addition, D_R^d, u(P)= − (D_L^d, u(P))^*.

Notice that the Wilson coefficients in Eq. include all pseudoscalar and axial contributions to τ → ℓP. Therefore, this goes beyond some well-known results in the literature, see for example , where box contributions were neglected.

NameProcess = "TauLMeson";
NameObservables = {{BrTautoEPi, 2001, "BR(tau->e pi)"},
                   {BrTautoEEta, 2002, "BR(tau->e eta)"},
                   {BrTautoEEtap, 2003, "BR(tau->e eta')"},
           {BrTautoMuPi, 2004, "BR(tau->mu pi)"},
                   {BrTautoMuEta, 2005, "BR(tau->mu eta)"},
                   {BrTautoMuEtap, 2006, "BR(tau->mu eta')"}};

NeededOperators = {OllddSLL, OllddSRR, OllddSRL, OllddSLR,
  OllddVRR, OllddVLL, OllddVRL, OllddVLR,
  OlluuSLL, OlluuSRR, OlluuSRL, OlluuSLR,
  OlluuVRR, OlluuVLL, OlluuVRL, OlluuVLR
};

Body = "TauLMeson.f90";

Real(dp) :: Fpi, thetaEta, mPi, mK, mEta, mEtap, meson_abs_T2, cont, &
     & mP, CP, factor, BR
Complex(dp) :: BSLL, BSLR, BSRL, BSRR, BVLL, BVLR, BVRL, BVRR, &
     & CSLL, CSLR, CSRL, CSRR, CVLL, CVLR, CVRL, CVRR, aP(2), bP(2), &
     & DLdP, DRdP, DLuP, DRuP
Integer :: i1, i2, out, k1, k2

! ----------------------------------------------------------------
! tau -> l meson
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Generalizes the analytical expressions in
! E. Arganda et al, JHEP 0806 (2008) 079 [arXiv:0803.2039]
! ----------------------------------------------------------------

Fpi=0.0924_dp! Pion decay constant in GeV
thetaEta=-Pi/10._dp! eta-eta' mixing angle
mPi=0.13497_dp! Pion mass in GeV
mK=0.49761_dp! Kaon mass in GeV
mEta=0.548_dp! Eta mass in GeV
mEtap=0.958_dp! Eta' mass in GeV

!Mesons:
!1:Pi0
!2:Eta
!3:Eta'
Do i1=1,3
   If(i1.eq.1) Then !1:Pi0
      mP = mPi
      CP = 1._dp
      DLdP = - mPi**2/4._dp
      DRdP = - Conjg(DLdP)
      DLuP = mPi**2/4._dp
      DRuP = - Conjg(DLuP)
   Else If(i1.eq.2) Then !2:Eta
      mP = mEta
      CP = (Sin(thetaEta)+Sqrt(2._dp)*Cos(thetaEta))/Sqrt(6._dp)
      DLdP = 1._dp/(4._dp*Sqrt(3._dp))*((3._dp*mPi**2-4._dp*mK**2) &
        & *Cos(thetaEta)-2._dp*Sqrt(2._dp)*mK**2*Sin(thetaEta))
      DRdP = - Conjg(DLdP)
      DLuP = 1._dp/(4._dp*Sqrt(3._dp))*mPi**2*(Cos(thetaEta)       &
        & -Sqrt(2._dp)*Sin(thetaEta))
      DRuP = - Conjg(DLuP)
   Else If(i1.eq.3) Then !3:Eta'
      mP = mEtap
      CP = (Sqrt(2._dp)*Sin(thetaEta)-Cos(thetaEta))/Sqrt(6._dp)
      DLdP = 1._dp/(4._dp*Sqrt(3._dp))*((3._dp*mPi**2-4._dp*mK**2) &
        & *Sin(thetaEta)+2._dp*Sqrt(2._dp)*mK**2*Cos(thetaEta))
      DRdP = - Conjg(DLdP)
      DLuP = 1._dp/(4._dp*Sqrt(3._dp))*mPi**2*(Sin(thetaEta)+      &
        & Sqrt(2._dp)*Cos(thetaEta))
      DRuP = - Conjg(DLuP)
   End If

!Leptons:
!1:e
!2:mu
Do i2=1,2
If (i2.eq.1) Then         ! tau -> e P
 out = 1
Elseif (i2.eq.2) Then     ! tau -> mu P
 out = 2
End if

! d-quark coefficients

BSLL = OllddSLL(3,out,1,1)
BSLR = OllddSLR(3,out,1,1)
BSRL = OllddSRL(3,out,1,1)
BSRR = OllddSRR(3,out,1,1)
BVLL = OllddVLL(3,out,1,1)
BVLR = OllddVLR(3,out,1,1)
BVRL = OllddVRL(3,out,1,1)
BVRR = OllddVRR(3,out,1,1)

! u-quark coefficients

CSLL = OlluuSLL(3,out,1,1)
CSLR = OlluuSLR(3,out,1,1)
CSRL = OlluuSRL(3,out,1,1)
CSRR = OlluuSRR(3,out,1,1)
CVLL = OlluuVLL(3,out,1,1)
CVLR = OlluuVLR(3,out,1,1)
CVRL = OlluuVRL(3,out,1,1)
CVRR = OlluuVRR(3,out,1,1)

! aP, bP scalar
aP(1) = Fpi/2._dp*(DLdP/mf_d(1)*(BSLL+BSRL) + DRdP/mf_d(1)*(BSLR+BSRR) &
         & + DLuP/mf_u(1)*(CSLL+CSRL) + DRuP/mf_u(1)*(CSLR+CSRR))
bP(1) = Fpi/2._dp*(DLdP/mf_d(1)*(BSRL-BSLL) + DRdP/mf_d(1)*(BSRR-BSLR) &
         & + DLuP/mf_u(1)*(CSRL-CSLL) + DRuP/mf_u(1)*(CSRR-CSLR))

! aP, bP vector
aP(2) = Fpi/4._dp*CP*(mf_l(3)-mf_l(out))*(-BVLL+BVLR-BVRL+BVRR+        &
         & CVLL-CVLR+CVRL-CVRR)
bP(2) = Fpi/4._dp*CP*(mf_l(3)+mf_l(out))*(-BVLL+BVLR+BVRL-BVRR+        &
         & CVLL-CVLR-CVRL+CVRR)

! averaged squared amplitude
meson_abs_T2=0._dp
Do k1=1,2
   Do k2=1,2
      cont=2._dp*mf_l(out)*mf_l(3)*(aP(k1)*conjg(aP(k2))               &
           & -bP(k1)*conjg(bP(k2)))+                                   &
           & (mf_l(3)**2+mf_l(out)**2-mP**2)*(aP(k1)*conjg(aP(k2))+    &
           & bP(k1)*conjg(bP(k2)))
      meson_abs_T2=meson_abs_T2+cont
   End Do
End Do
meson_abs_T2=meson_abs_T2/(2._dp*mf_l(3))

! branching ratio
factor=oo4pi*Sqrt(lamb(mf_l(3)**2,mf_l(out)**2,mP**2))                 &
            & /(mf_l(3)**2*GammaTau)*0.5_dp
BR=factor*meson_abs_T2
If (i1.eq.1) Then !pi
   If(i2.eq.1) Then
      BrTautoEPi = BR
   Else
      BrTautoMuPi = BR
   End If
Elseif (i1.eq.2) Then !eta
   If(i2.eq.1) Then
      BrTautoEEta = BR
   Else
      BrTautoMuEta = BR
   End If
Else !eta'
   If(i2.eq.1) Then
      BrTautoEEtap = BR
   Else
      BrTautoMuEtap = BR
   End If
End if

End Do
End Do

Contains

Real(dp) Function lamb(x,y,z)
Real(dp),Intent(in)::x,y,z
 lamb=(x+y-z)**2-4._dp*x*y
End Function lamb

Higgs two-body decays (h → ℓ_αℓ_β)

The decay width is given by 

\Gamma \left( h \to \ell_\alpha \ell_\beta \right) \equiv \Gamma \left( h \to \ell_\alpha \bar \ell_\beta \right) + \Gamma \left( h \to \bar \ell_\alpha \ell_\beta \right) = \\ \frac{1}{16 \pi m_h} \left[ \left(1-\left(\frac{m_{\ell_\alpha} + m_{\ell_\beta}}{m_h}\right)^2\right)\left(1-\left(\frac{m_{\ell_\alpha} - m_{\ell_\beta}}{m_h}\right)^2\right)\right]^{1/2} \nonumber \\ \times \left[ \left( m_h^2 - m_{\ell_\alpha}^2 - m_{\ell_\beta}^2 \right) \left( |S_L|^2 + |S_R|^2 \right)_{\alpha \beta} - 4 m_{\ell_\alpha} m_{\ell_\beta} \text{Re}(S_L S_R^\ast)_{\alpha \beta} \right] \nonumber \\ + (\alpha \leftrightarrow \beta) \nonumber

NameProcess = "hLLp";
NameObservables = {{BrhtoMuE, 1101, "BR(h->e mu)"},
                   {BrhtoTauE, 1102, "BR(h->e tau)"},
                   {BrhtoTauMu, 1103, "BR(h->mu tau)"}};

NeededOperators = {OH2lSL, OH2lSR};

Body = "hLLp.f90";

Real(dp) :: width1, width2, width, mh, gamh, kinfactor
Complex(dp) :: SL1, SR1, SL2, SR2
Integer :: i1, gt1, gt2, hLoc

! ----------------------------------------------------------------
! h -> l l'
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on E. Arganda et al, PRD 71 (2005) 035011 [hep-ph/0407302]
! ----------------------------------------------------------------

!! NEXT LINE HAVE TO BE PARSED BY SARAH
! Checking if there are several generations of Scalars and what is the SM-like doublet
@ If[getGen[HiggsBoson]>1, "hLoc = MaxLoc(Abs("<>ToString[HiggsMixingMatrix]<>"(2,:)),1)", "hLoc = 1"]

@ "mh = "<>ToString[SPhenoMass[HiggsBoson]]<>If[getGen[HiggsBoson]>1, "(hLoc)", ""]

@ "gamh ="<>ToString[SPhenoWidth[HiggsBoson]]<>If[getGen[HiggsBoson]>1, "(hLoc)", ""]

If (.not.L_BR) gamh = 4.5E-3_dp  ! Decays not calculated; using SM value

Do i1=1,3

If (i1.eq.1) Then         ! h -> e mu
 gt1 = 1
 gt2 = 2
Elseif (i1.eq.2) Then     ! h -> e tau
 gt1 = 1
 gt2 = 3
Else                      ! h -> mu tau
 gt1 = 2
 gt2 = 3
End if

! width = Gamma(h -> \bar{l1} l2) + Gamma(h -> l1 \bar{l2})

SL1 = OH2lSL(gt1,gt2,hLoc)
SR1 = OH2lSR(gt1,gt2,hLoc)
SL2 = OH2lSL(gt2,gt1,hLoc)
SR2 = OH2lSR(gt2,gt1,hLoc)

kinfactor = (1-(mf_l(gt1)+mf_l(gt2)/mh)**2)*&
       & (1-(mf_l(gt1)-mf_l(gt2)/mh)**2)

width1 = (mh**2-mf_l(gt1)**2-mf_l(gt2)**2)*(Abs(SL1)**2+Abs(SR1)**2) &
     & - 4._dp*mf_l(gt1)*mf_l(gt2)*Real(SL1*Conjg(SR1),dp)
width2 = (mh**2-mf_l(gt1)**2-mf_l(gt2)**2)*(Abs(SL2)**2+Abs(SR2)**2) &
     & - 4._dp*mf_l(gt1)*mf_l(gt2)*Real(SL2*Conjg(SR2),dp)

! decay width
width = oo16pi/mh * sqrt(kinfactor) * (width1+width2)

If (i1.eq.1) Then
BrhtoMuE = width/(width+gamh)
Elseif (i1.eq.2) Then
BrhtoTauE = width/(width+gamh)
Else
BrhtoTauMu = width/(width+gamh)
End if

End do

Z two-body decays (Z → ℓ_αℓ_β)

The decay width is given by 

\Gamma \left( Z \to \ell_\alpha \ell_\beta \right) \equiv \Gamma \left( Z \to \ell_\alpha \bar \ell_\beta \right) + \Gamma \left( Z \to \bar \ell_\alpha \ell_\beta \right) = \\ \frac{m_Z}{48 \pi} \left[ 2 \left( |R_1^L|^2 + |R_1^R|^2 \right) + \frac{m_Z^2}{4} \left( |R_2^L|^2 + |R_2^R|^2 \right) \right] \, , \nonumber

where the charged lepton masses have been neglected.

NameProcess = "ZLLp";
NameObservables = {{BrZtoMuE, 1001, "BR(Z->e mu)"},
                   {BrZtoTauE, 1002, "BR(Z->e tau)"},
                   {BrZtoTauMu, 1003, "BR(Z->mu tau)"}};

NeededOperators = {OZ2lSL, OZ2lSR,OZ2lVL,OZ2lVR};

Body = "ZLLp.f90";

Real(dp) :: width
Integer :: i1, gt1, gt2

! ----------------------------------------------------------------
! Z -> l l'
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on X. -J. Bi et al, PRD 63 (2001) 096008 [hep-ph/0010270]
! ----------------------------------------------------------------

Do i1=1,3

If (i1.eq.1) Then         ! Z -> e mu
 gt1 = 1
 gt2 = 2
Elseif (i1.eq.2) Then     !Z -> e tau
 gt1 = 1
 gt2 = 3
Else                      ! Z -> mu tau
 gt1 = 2
 gt2 = 3
End if

! decay width
width = oo48pi*(2*(Abs(OZ2lVL(gt1,gt2))**2 +            &
  &                    Abs(OZ2lVR(gt1,gt2))**2)*mZ      &
  & + (Abs(OZ2lSL(gt1,gt2))**2+Abs(OZ2lSR(gt1,gt2))**2) &
  & * mZ * mZ2 * 0.25_dp)

If (i1.eq.1) Then
BrZtoMuE = width/(width+gamZ)
Elseif (i1.eq.2) Then
BrZtoTauE = width/(width+gamZ)
Else
BrZtoTauMu = width/(width+gamZ)
End if

End do

Quark flavor observables

QFV has been observed and its description in the SM due to the CKM matrix is well established. However, the large majority of BSM models causes additional contributions which have to be studied carefully, see for instance Refs. .

We give also here a description of the implementation of the different observables using the operators present in the SPheno output of SARAH.

Neutral B-meson decays in charged leptons (B_s**,** d⁰ → ℓ⁺ℓ⁻)

Our analytical results for B_s, d⁰ → ℓ⁺ℓ⁻ follow . The B⁰ ≡ B_s, d⁰ decay width to a pair of charged leptons can be written as

\Gamma \left(B^0 \to \ell_\alpha^+ \ell_\beta^- \right) = \frac{|\mathcal{M_{B\ell\ell}}|^2}{16 \pi M_{B}} \left[ \left(1-\left(\frac{m_{\ell_\alpha} + m_{\ell_\beta}}{m_B}\right)^2\right)\left(1-\left(\frac{m_{\ell_\alpha} - m_{\ell_\beta}}{m_B}\right)^2\right)\right]^{1/2} \, .

Here

|\mathcal{M_{B\ell\ell}}|^2 = 2 |F_S|^2 \left[ m_B^2 - \left(m_{\ell_\alpha} + m_{\ell_\beta}\right)^2 \right] + 2 |F_P|^2 \left[ m_B^2 - \left(m_{\ell_\alpha} - m_{\ell_\beta}\right)^2 \right] \nonumber \\ + 2 |F_V|^2 \left[ m_B^2 \left(m_{\ell_\alpha} - m_{\ell_\beta}\right)^2 - \left(m_{\ell_\alpha}^2 - m_{\ell_\beta}^2\right)^2 \right] \nonumber \\ + 2 |F_A|^2 \left[ m_B^2 \left(m_{\ell_\alpha} + m_{\ell_\beta}\right)^2 - \left(m_{\ell_\alpha}^2 - m_{\ell_\beta}^2\right)^2 \right] \nonumber \\ + 4 \, \text{Re}(F_S F_V^\ast) \left(m_{\ell_\alpha} - m_{\ell_\beta}\right) \left[ m_B^2 + \left(m_{\ell_\alpha} + m_{\ell_\beta}\right)^2 \right] \nonumber \\ + 4 \, \text{Re}(F_P F_A^\ast) \left(m_{\ell_\alpha} + m_{\ell_\beta}\right) \left[ m_B^2 - \left(m_{\ell_\alpha} - m_{\ell_\beta}\right)^2 \right] \, ,

and the F_X coefficients are defined in terms of our Wilson coefficients as[1]

F_S = \frac{i}{4} \frac{m_B^2 f_B}{m_d + m_{d^\prime}} \left( E_{LL}^S + E_{LR}^S - E_{RR}^S - E_{RL}^S \right) \\ F_P = \frac{i}{4} \frac{m_B^2 f_B}{m_d + m_{d^\prime}} \left( - E_{LL}^S + E_{LR}^S - E_{RR}^S + E_{RL}^S \right) \\ F_V = - \frac{i}{4} f_B \left( E_{LL}^V + E_{LR}^V - E_{RR}^V - E_{RL}^V \right) \\ F_A = - \frac{i}{4} f_B \left( - E_{LL}^V + E_{LR}^V - E_{RR}^V + E_{RL}^V \right) \, ,

where f_B ≡ f_Bd, s⁰ is the B_d, s⁰ decay constant and m_d, d^′ are the masses of the quarks contained in the B meson, B_d⁰ ≡ b̄**d and B_s⁰ ≡ b̄**s. In the lepton flavor conserving case, α = β, the F_V contribution vanishes. In this case, the results in are in agreement with previous computations .

NameProcess = "B0toLL";
NameObservables = {{BrB0dEE, 4000, "BR(B^0_d->e e)"},
                   {ratioB0dEE, 4001, "BR(B^0_d->e e)/BR(B^0_d->e e)_SM"},
                   {BrB0sEE, 4002, "BR(B^0_s->e e)"},
                   {ratioB0sEE, 4003, "BR(B^0_s->e e)/BR(B^0_s->e e)_SM"},
                   {BrB0dMuMu, 4004, "BR(B^0_d->mu mu)"},
                   {ratioB0dMuMu, 4005, "BR(B^0_d->mu mu)/BR(B^0_d->mu mu)_SM"},
                   {BrB0sMuMu, 4006, "BR(B^0_s->mu mu)"},
                   {ratioB0sMuMu, 4007, "BR(B^0_s->mu mu)/BR(B^0_s->mu mu)_SM"},
                   {BrB0dTauTau, 4008, "BR(B^0_d->tau tau)"},
                   {ratioB0dTauTau, 4009, "BR(B^0_d->tau tau)/BR(B^0_d->tau tau)_SM"},
                   {BrB0sTauTau, 4010, "BR(B^0_s->tau tau)"},
                   {ratioB0sTauTau, 4011, "BR(B^0_s->tau tau)/BR(B^0_s->tau tau)_SM"} };


NeededOperators = {OddllSLL, OddllSRR, OddllSRL, OddllSLR,
                   OddllVRR, OddllVLL, OddllVRL, OddllVLR,
                   OddllSLLSM, OddllSRRSM, OddllSRLSM, OddllSLRSM,
                   OddllVRRSM, OddllVLLSM, OddllVRLSM, OddllVLRSM};

Body = "B0ll.f90";

Real(dp) :: AmpSquared,AmpSquared2, AmpSquared_SM, AmpSquared2_SM, &
              & width_SM, width
Real(dp) :: MassB0s, MassB0d, fBs, fBd, TauB0s, TauB0d
Real(dp) :: hbar=6.58211899E-25_dp
Real(dp) :: MassB0,MassB02,fB0,GammaB0
Complex(dp) :: CS(4), CV(4), CT(4)
Complex(dp) :: FS=0._dp, FP=0._dp, FV=0._dp, FA=0._dp
Integer :: i1, gt1, gt2, gt3, gt4

! ----------------------------------------------------------------
! B0 -> l l
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on A. Dedes et al, PRD 79 (2009) 055006 [arXiv:0812.4320]
! ----------------------------------------------------------------

! Using global hadronic data
fBd = f_B0d_CONST
fBs = f_B0s_CONST
TauB0d = tau_B0d
TauB0s = tau_B0s
MassB0d = mass_B0d
MassB0s = mass_B0s

Do i1=1,6
gt1 = 3
If (i1.eq.1) Then ! B0d -> e+ e-
  MassB0 = MassB0d
  MassB02 = MassB0d**2
  fB0 = fBd
  GammaB0 = (hbar)/(TauB0d)
  gt2 = 1
  gt3 = 1
  gt4 = 1
Else if (i1.eq.2) Then ! B0s -> e+ e-
  MassB0 = MassB0s
  MassB02 = MassB0s**2
  fB0 = fBs
  GammaB0 = (hbar)/(TauB0s)
  gt2 = 2
  gt3 = 1
  gt4 = 1
Else if (i1.eq.3) Then ! B0d -> mu+ mu-
  MassB0 = MassB0d
  MassB02 = MassB0d**2
  fB0 = fBd
  GammaB0 = (hbar)/(TauB0d)
  gt2 = 1
  gt3 = 2
  gt4 = 2
Else if (i1.eq.4) Then ! B0s -> mu+ mu-
  MassB0 = MassB0s
  MassB02 = MassB0s**2
  fB0 = fBs
  GammaB0 = (hbar)/(TauB0s)
  gt2 = 2
  gt3 = 2
  gt4 = 2
Else if (i1.eq.5) Then ! B0d -> tau+ tau-
  MassB0 = MassB0d
  MassB02 = MassB0d**2
  fB0 = fBd
  GammaB0 = (hbar)/(TauB0d)
  gt2 = 1
  gt3 = 3
  gt4 = 3
Else if (i1.eq.6) Then ! B0s -> tau+ tau-
  MassB0 = MassB0s
  MassB02 = MassB0s**2
  fB0 = fBs
  GammaB0 = (hbar)/(TauB0s)
  gt2 = 2
  gt3 = 3
  gt4 = 3
End if

! BSM contributions

CS(1) = OddllSRR(gt1,gt2,gt3,gt4)
CS(2) = OddllSRL(gt1,gt2,gt3,gt4)
CS(3) = OddllSLL(gt1,gt2,gt3,gt4)
CS(4) = OddllSLR(gt1,gt2,gt3,gt4)

CV(1) = OddllVLL(gt1,gt2,gt3,gt4)
CV(2) = OddllVLR(gt1,gt2,gt3,gt4)
CV(3) = OddllVRR(gt1,gt2,gt3,gt4)
CV(4) = OddllVRL(gt1,gt2,gt3,gt4)

FS= 0.25_dp*MassB02*fB0/(MFd(gt1)+MFd(gt2))*( CS(1)+CS(2)-CS(3)-CS(4))
FP= 0.25_dp*MassB02*fB0/(MFd(gt1)+MFd(gt2))*(-CS(1)+CS(2)-CS(3)+CS(4))
FV= -0.25_dp*fB0*( CV(1)+CV(2)-CV(3)-CV(4))
FA= -0.25_dp*fB0*(-CV(1)+CV(2)-CV(3)+CV(4))

AmpSquared = 2 * abs(FS)**2 * (MassB02 - (mf_l(gt3)+mf_l(gt4))**2) &
     & + 2 *abs(FP)**2 * (MassB02 - (mf_l(gt3)-mf_l(gt4))**2) &
     & + 2 *abs(FV)**2 * (MassB02*(mf_l(gt4)-mf_l(gt3))**2    &
             & - (mf_l2(gt4)-mf_l2(gt3))**2) &
     & + 2 *abs(FA)**2 * (MassB02*(mf_l(gt4)+mf_l(gt3))**2 -  &
             & (mf_l2(gt4)-mf_l2(gt3))**2) &
     & + 4 *REAL(FS*conjg(FV)) *(mf_l(gt3)-mf_l(gt4)) *(MassB02 &
             & + (mf_l(gt3)+mf_l(gt4))**2) &
     & + 4 *REAL(FP*conjg(FA)) *(mf_l(gt3)+mf_l(gt4)) *(MassB02 &
             & - (mf_l(gt3)-mf_l(gt4))**2)

width = oo16pi * AmpSquared / MassB0 * &
    & sqrt(1-((mf_l(gt4)+mf_l(gt3))/MassB0)**2) &
    & * sqrt(1-((mf_l(gt4)-mf_l(gt3))/MassB0)**2)*(Alpha/Alpha_160)**4


! SM contributions

CS(1) = OddllSRRSM(gt1,gt2,gt3,gt4)
CS(2) = OddllSRLSM(gt1,gt2,gt3,gt4)
CS(3) = OddllSLLSM(gt1,gt2,gt3,gt4)
CS(4) = OddllSLRSM(gt1,gt2,gt3,gt4)

CV(1) = OddllVLLSM(gt1,gt2,gt3,gt4)
CV(2) = OddllVLRSM(gt1,gt2,gt3,gt4)
CV(3) = OddllVRRSM(gt1,gt2,gt3,gt4)
CV(4) = OddllVRLSM(gt1,gt2,gt3,gt4)

FS= 0.25_dp*MassB02*fB0/(MFd(gt1)+MFd(gt2))*( CS(1)+CS(2)-CS(3)-CS(4))
FP= 0.25_dp*MassB02*fB0/(MFd(gt1)+MFd(gt2))*(-CS(1)+CS(2)-CS(3)+CS(4))
FV= -0.25_dp*fB0*( CV(1)+CV(2)-CV(3)-CV(4))
FA= -0.25_dp*fB0*(-CV(1)+CV(2)-CV(3)+CV(4))

AmpSquared = 2 * abs(FS)**2 * (MassB02 - (mf_l(gt3)+mf_l(gt4))**2) &
     & + 2 *abs(FP)**2 * (MassB02 - (mf_l(gt3)-mf_l(gt4))**2) &
     & + 2 *abs(FV)**2 * (MassB02*(mf_l(gt4)-mf_l(gt3))**2 -  &
        &  (mf_l2(gt4)-mf_l2(gt3))**2) &
     & + 2 *abs(FA)**2 * (MassB02*(mf_l(gt4)+mf_l(gt3))**2 -  &
        &  (mf_l2(gt4)-mf_l2(gt3))**2) &
     & + 4 *REAL(FS*conjg(FV)) *(mf_l(gt3)-mf_l(gt4)) *(MassB02 &
        & + (mf_l(gt3)+mf_l(gt4))**2) &
     & + 4 *REAL(FP*conjg(FA)) *(mf_l(gt3)+mf_l(gt4)) *(MassB02 &
        & - (mf_l(gt3)-mf_l(gt4))**2)

width_SM = oo16pi * AmpSquared / MassB0 * sqrt(1-((mf_l(gt4)+   &
     & mf_l(gt3))/MassB0)**2) &
     & * sqrt(1-((mf_l(gt4)-mf_l(gt3))/MassB0)**2)*(Alpha/Alpha_160)**4


If (i1.Eq.1) Then
  BrB0dEE= width / GammaB0
  ratioB0dEE= width / width_SM
Else If (i1.Eq.2) Then
  BrB0sEE= width / GammaB0
  ratioB0sEE= width / width_SM
Else If (i1.Eq.3) Then
  BrB0dMuMu= width / GammaB0
  ratioB0dMuMu= width / width_SM
Else If (i1.Eq.4) Then
  BrB0sMuMu= width / GammaB0
  ratioB0sMuMu= width / width_SM
Else If (i1.Eq.5) Then
  BrB0dTauTau= width / GammaB0
  ratioB0dTauTau= width / width_SM
Else If (i1.Eq.6) Then
  BrB0sTauTau= width / GammaB0
  ratioB0sTauTau= width / width_SM
End If

End do

Radiative B-meson decay (B̄ → X_sγ)

The branching ratio for B̄ → X_sγ, with a cut E_γ > 1.6 GeV in the B̄ rest frame, can be obtained as 

\BR \hspace\*{-0.1cm} \left( \bar B \to X_s \gamma \right)_{E_\gamma gt; 1.6 \text{GeV}} = 10^{-4} \bigg\[ a_{SM} + a_{77} \left( |\delta C_7^{(0)}|^2 + |\delta C_7^{\prime (0)}|^2 \right) + a_{88} \left( |\delta C_8^{(0)}|^2 + |\delta C_8^{\prime (0)}|^2 \right) \nonumber \\ + \text{Re} \left( a_7 \, \delta C_7^{(0)} + a_8 \, \delta C_8^{(0)} + a_{78} \left( \delta C_7^{(0)} \, \delta C_8^{(0) \ast} + \delta C_7^{\prime (0)} \, \delta C_8^{\prime (0) \, \ast} \right) \right) \bigg\] \, ,

where a_S**M = 3.15 is the NNLO SM prediction , the other a coefficients in Eq. are found to be

a_{77} = 4.743 \nonumber \\ a_{88} = 0.789 \nonumber \\ a_7 = -7.184 + 0.612 \, i \nonumber \\ a_8 = -2.225 - 0.557 \, i \nonumber \\ a_{78} = 2.454 - 0.884 \, i

and we have defined δ**C_i⁽⁰⁾ = C_i⁽⁰⁾ − C_i^(0) SM. Finally, the C_i⁽⁰⁾ coefficients can be written in terms of Q_1, 2^L, R in Eqs. and as

C_7^{(0)} = n_{CQ} \, Q_1^R \\ C_7^{\prime (0)} = n_{CQ} \, Q_1^L \\ C_8^{(0)} = n_{CQ} \, Q_2^R \\ C_8^{\prime (0)} = n_{CQ} \, Q_2^L

where n_{CQ}^{-1} = - \frac{G_F}{4 \sqrt{2} \pi^2} V_{tb} V_{ts}^\ast and V is the Cabibbo-Kobayashi-Maskawa (CKM) matrix.

NameProcess = "bsGamma";
NameObservables = {{BrBsGamma, 200, "BR(B->X_s gamma)"},
                   {ratioBsGamma, 201, "BR(B->X_s gamma)/BR(B->X_s gamma)_SM"}};

NeededOperators = {CC7, CC7p, CC8, CC8p,
                   CC7SM, CC7pSM, CC8SM, CC8pSM};

Body = "bsGamma.f90";

Integer :: gt1, gt2
Complex(dp) :: norm, delta_C7_0, delta_C7p_0, delta_C8_0, delta_C8p_0
Real(dp) :: NNLO_SM

! ----------------------------------------------------------------
! \bar{B} -> X_s gamma (Egamma > 1.6 GeV)
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on E. Lunghi, J. Matias, JHEP 0704 (2007) 058 [hep-ph/0612166]
! ----------------------------------------------------------------

gt1=3 !b
gt2=2 !s

! normalization of our Wilson coefficients
! relative to the ones used in hep-ph/0612166
norm = -CKM_160(3,3)*Conjg(CKM_160(gt1,gt2))*Alpha_160/ &
         & (8._dp*Pi*sinW2_160*mW2)

! Wilson coefficients
delta_C7_0 =(CC7(gt1,gt2)-CC7SM(gt1,gt2))/norm
delta_C7p_0=(CC7p(gt1,gt2)-CC7pSM(gt1,gt2))/norm
delta_C8_0 =(CC8(gt1,gt2)-CC8SM(gt1,gt2))/norm
delta_C8p_0=(CC8p(gt1,gt2)-CC8pSM(gt1,gt2))/norm

! NNLO SM prediction
! as obtained in M. Misiak et al, PRL 98 (2007) 022002
! and M. Misiak and M. Steinhauser, NPB 764 (2007) 62
NNLO_SM=3.15_dp

BrBsGamma=NNLO_SM+4.743_dp*(Abs(delta_C7_0)**2+Abs(delta_C7p_0)**2)&
&+0.789_dp*(Abs(delta_C8_0)**2+Abs(delta_C8p_0)**2)&
&+Real((-7.184_dp,0.612_dp)*delta_C7_0&
&+(-2.225_dp,-0.557_dp)*delta_C8_0+(2.454_dp,-0.884_dp)*&
&(delta_C7_0*conjg(delta_C8_0)+delta_C7p_0*conjg(delta_C8p_0)),dp)

! ratio BSM/SM
ratioBsGamma = BrBsGamma/NNLO_SM

! branching ratio
BrBsGamma=1E-4_dp*BrBsGamma

B-meson decay into a S-meson and two charged leptons (B̄ → X_sℓ⁺ℓ⁻)

Our results for B̄ → X_sℓ⁺ℓ⁻ are based on , expanded with the addition of prime operators contributions . The branching ratios for the ℓ = e case can be written as

10^7 \, BR \left( \bar B \to X_s e^+ e^- \right) = 2.3148 - 0.001658 Im (R_{10}) + 0.0005 Im(R_{10} R_8^\ast + R_{10}^\prime R_8^{\prime \, \ast}) \nonumber \\ + 0.0523 Im(R_7) + 0.02266 Im(R_7 R_8^\ast + R_7^\prime R_8^{\prime \, \ast}) + 0.00496 Im(R_7 R_9^\ast + R_7^\prime R_9^{\prime \, \ast}) \nonumber \\ + 0.00518 Im(R_8) + 0.0261 Im(R_8 R_9^\ast + R_8^\prime R_9^{\prime \, \ast}) - 0.00621 Im(R_9) - 0.5420 Re(R_{10}) \nonumber \\ - 0.03340 Re(R_7) + 0.0153 Re(R_7 R_{10}^\ast + R_7^\prime R_{10}^{\prime \, \ast}) + 0.0673 Re(R_7 R_8^\ast + R_7^\prime R_8^{\prime \, \ast}) \nonumber \\ - 0.86916 Re(R_7 R_9^\ast + R_7^\prime R_9^{\prime \, \ast}) - 0.0135 Re(R_8) + 0.00185 Re(R_8 R_{10} + R_8^\prime R_{10}^{\prime \, \ast}) \nonumber \\ - 0.09921 Re(R_8 R_9^\ast + R_8^\prime R_9^{\prime \, \ast}) + 2.833 Re(R_9) - 0.10698 Re(R_9 R_{10}^\ast + R_9^\prime R_{10}^{\prime \, \ast}) \nonumber \\ + 11.0348 \left( |R_{10}|^2 + |R_{10}^\prime|^2 \right) + 0.2804 \left( |R_7|^2 + |R_7^\prime|^2 \right) \nonumber \\ + 0.003763 \left( |R_8|^2 + |R_8^\prime|^2 \right) + 1.527 \left( |R_9|^2 + |R_9^\prime|^2 \right) \, ,

whereas for the ℓ = μ case one gets

10^7 \, BR \left( \bar B \to X_s \mu^+ \mu^- \right) = 2.1774 - 0.001658 Im (R_{10}) + 0.0005 Im(R_{10} R_8^\ast + R_{10}^\prime R_8^{\prime \, \ast}) \nonumber \\ + 0.0534 Im(R_7) + 0.02266 Im(R_7 R_8^\ast + R_7^\prime R_8^{\prime \, \ast}) + 0.00496 Im(R_7 R_9^\ast + R_7^\prime R_9^{\prime \, \ast}) \nonumber \\ + 0.00527 Im(R_8) + 0.0261 Im(R_8 R_9^\ast + R_8^\prime R_9^{\prime \, \ast}) - 0.0115 Im(R_9) - 0.5420 Re(R_{10}) \nonumber \\ + 0.0208 Re(R_7) + 0.0153 Re(R_7 R_{10}^\ast + R_7^\prime R_{10}^{\prime \, \ast}) + 0.0648 Re(R_7 R_8^\ast + R_7^\prime R_8^{\prime \, \ast}) \nonumber \\ - 0.8545 Re(R_7 R_9^\ast + R_7^\prime R_9^{\prime \, \ast}) - 0.00938 Re(R_8) + 0.00185 Re(R_8 R_{10} + R_8^\prime R_{10}^{\prime \, \ast}) \nonumber \\ - 0.0981 Re(R_8 R_9^\ast + R_8^\prime R_9^{\prime \, \ast}) + 2.6917 Re(R_9) - 0.10698 Re(R_9 R_{10}^\ast + R_9^\prime R_{10}^{\prime \, \ast}) \nonumber \\ + 10.7652 \left( |R_{10}|^2 + |R_{10}^\prime|^2 \right) + 0.2880 \left( |R_7|^2 + |R_7^\prime|^2 \right) \nonumber \\ + 0.003763 \left( |R_8|^2 + |R_8^\prime|^2 \right) + 1.527 \left( |R_9|^2 + |R_9^\prime|^2 \right) \, .

Here we have defined the ratios of Wilson coefficients

R_{7,8} = \frac{Q_{1,2}^R}{Q_{1,2}^{R, \text{SM}}} \qquad , \qquad R_{7,8}^\prime = \frac{Q_{1,2}^L}{Q_{1,2}^{L, \text{SM}}}

as well as

R_{9,10} = \frac{E_{LL}^V \pm E_{LR}^V}{E_{LL}^{V, \text{SM}} \pm E_{LR}^{V, \text{SM}}} \qquad , \qquad R_{9,10}^\prime = \frac{E_{RR}^V \pm E_{RL}^V}{E_{RR}^{V, \text{SM}} \pm E_{RL}^{V, \text{SM}}} \, .

NameProcess = "BtoSLL";
NameObservables = {{BrBtoSEE, 5000, "BR(B-> s e e)"},
                   {ratioBtoSEE, 5001, "BR(B-> s e e)/BR(B-> s e e)_SM"},
                   {BrBtoSMuMu, 5002, "BR(B-> s mu mu)"} ,
                   {ratioBtoSMuMu, 5003, "BR(B-> s mu mu)/BR(B-> s mu mu)_SM"}};

NeededOperators = {OddllVRR, OddllVLL, OddllVRL, OddllVLR,
                   CC7, CC7p, CC8, CC8p,
                   OddllVRRSM, OddllVLLSM, OddllVRLSM, OddllVLRSM,
                   CC7SM, CC7pSM, CC8SM, CC8pSM
};

Body = "BtoSLL.f90";

Complex(dp) :: c7(2), c7p(2), c8(2), c8p(2), r7, r7p, r8, r8p, norm, &
     & r9(2), r9p(2), r10(2), r10p(2),                               &
     & c9ee(2), c9pee(2), c10ee(2), c10pee(2),                       &
     & c9_cee(2), c9p_cee(2), c10_cee(2), c10p_cee(2),               &
     & c9mm(2), c9pmm(2), c10mm(2), c10pmm(2), c9_cmm(2),            &
     & c9p_cmm(2), c10_cmm(2), c10p_cmm(2)

! ----------------------------------------------------------------
! \bar{B} -> X_s l+ l-
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on T. Huber et al, NPB 740 (2006) 105, [hep-ph/0512066]
! Prime operators added after private communication with E. Lunghi
! ----------------------------------------------------------------

! Wilson coefficients

c7(1) = CC7(3,2)
c7(2) =  CC7SM(3,2)
c7p(1) = CC7p(3,2)
c7p(2) = CC7pSM(3,2)

c8(1) = CC8(3,2)
c8(2) = CC8SM(3,2)
c8p(1) = CC8p(3,2)
c8p(2) = CC8pSM(3,2)

c9ee(1) = OddllVLL(3,2,1,1)+OddllVLR(3,2,1,1)
c9ee(2) = (OddllVLLSM(3,2,1,1)+OddllVLRSM(3,2,1,1))
c9mm(1) = OddllVLL(3,2,2,2)+OddllVLR(3,2,2,2)
c9mm(2) =  (OddllVLLSM(3,2,2,2)+OddllVLRSM(3,2,2,2))
c9pee(1) = OddllVRR(3,2,1,1)+OddllVRL(3,2,1,1)
c9pee(2) =  (OddllVRRSM(3,2,1,1)+OddllVRLSM(3,2,1,1))
c9pmm(1) = OddllVRR(3,2,2,2)+OddllVRL(3,2,2,2)
c9pmm(2) =  (OddllVRRSM(3,2,2,2)+OddllVRLSM(3,2,2,2))

c10ee(1) = OddllVLL(3,2,1,1)-OddllVLR(3,2,1,1)
c10ee(2) =  (OddllVLLSM(3,2,1,1)-OddllVLRSM(3,2,1,1))
c10mm(1) = OddllVLL(3,2,2,2)-OddllVLR(3,2,2,2)
c10mm(2) =  (OddllVLLSM(3,2,2,2)-OddllVLRSM(3,2,2,2))
c10pee(1) = OddllVRR(3,2,1,1)-OddllVRL(3,2,1,1)
c10pee(2) =  (OddllVRRSM(3,2,1,1)-OddllVRLSM(3,2,1,1))
c10pmm(1) = OddllVRR(3,2,2,2)-OddllVRL(3,2,2,2)
c10pmm(2) = (OddllVRRSM(3,2,2,2)-OddllVRLSM(3,2,2,2))

! ratios

r7 = c7(1) / c7(2)
r7p = c7p(1) / c7(2)
r8 = c8(1) / c8(2)
r8p = c8p(1) / c8(2)

r9(1) = c9ee(1)/c9ee(2)
r9(2) = c9mm(1)/c9mm(2)
r9p(1) = c9pee(1)/c9ee(2)
r9p(2) = c9pmm(1)/c9mm(2)

r10(1) = c10ee(1)/c10ee(2)
r10(2) = c10mm(1)/c10mm(2)
r10p(1) = c10pee(1)/c10ee(2)
r10p(2) = c10pmm(1)/c10mm(2)

BrBtoSEE = (2.3148_dp - 1.658e-3_dp * Aimag(R10(1))                   &
 & + 5.e-4_dp * Aimag(r10(1)*Conjg(r8) + r10p(1)*Conjg(r8p) )         &
 & + 5.23e-2_dp * Aimag(r7) + 5.18e-3_dp * Aimag(r8)                  &
 & + 2.266e-2_dp * Aimag(r7 * Conjg(r8) + r7p * Conjg(r8p) )          &
 & + 4.96e-3_dp * Aimag(r7 * Conjg(r9(1)) + r7p * Conjg(r9p(1)) )     &
 & + 2.61e-2_dp * Aimag(r8 * Conjg(r9(1)) + r8p * Conjg(r9p(1)) )     &
 & - 6.21e-3_dp * Aimag(r9(1)) - 0.5420_dp * Real( r10(1), dp)        &
 & - 3.340e-2_dp * Real(r7,dp) - 1.35e-2_dp * Real(r8,dp)             &
 & + 1.53e-2_dp * Real(r7*Conjg(r10(1)) + r7p*Conjg(r10p(1)), dp )    &
 & + 6.73e-2_dp * Real(r7 * Conjg(r8) + r7p * Conjg(r8p), dp )        &
 & - 0.86916_dp * Real(r7*Conjg(r9(1)) + r7p*Conjg(r9p(1)), dp )      &
 & + 1.85e-3_dp * Real(r8*Conjg(r10(1)) + r8p*Conjg(r10p(1)), dp )    &
 & - 9.921e-2_dp * Real(r8*Conjg(r9(1)) + r8p*Conjg(r9p(1)), dp )     &
 & + 2.833_dp* Real(r9(1),dp) + 0.2804_dp * (Abs(r7)**2 + Abs(r7p)**2)&
 & - 0.10698_dp * Real( r9(1) * Conjg(r10(1))                         &
 &                    + r9p(1) * Conjg(r10p(1)), dp)                  &
 & + 11.0348_dp * (Abs(r10(1))**2 + Abs(r10p(1))**2 )                 &
 & + 1.527_dp * (Abs(r9(1))**2 + Abs(r9p(1))**2 )                     &
 & + 3.763e-3_dp * (Abs(r8)**2 + Abs(r8p)**2 ) )

  ! ratio BR(B -> Xs mu+ mu-)/BR(B -> Xs e+ e-)_SM
  ratioBtoSee = BrBtoSEE/16.5529_dp

  ! branching ratio B -> Xs e+ e-
  BrBtoSEE = BrBtoSEE* 1.e-7_dp

BrBtoSMuMu = (2.1774_dp - 1.658e-3_dp * Aimag(R10(2))                 &
  & + 5.e-4_dp * Aimag(r10(2)*Conjg(r8) + r10p(2)*Conjg(r8p) )        &
  & + 5.34e-2_dp * Aimag(r7) + 5.27e-3_dp * Aimag(r8)                 &
  & + 2.266e-2_dp * Aimag(r7 * Conjg(r8) + r7p * Conjg(r8p) )         &
  & + 4.96e-3_dp * Aimag(r7 * Conjg(r9(2)) + r7p * Conjg(r9p(2)) )    &
  & + 2.61e-2_dp * Aimag(r8 * Conjg(r9(2)) + r8p * Conjg(r9p(2)) )    &
  & - 1.15e-2_dp * Aimag(r9(2)) - 0.5420_dp * Real( r10(2), dp)       &
  & + 2.08e-2_dp * Real(r7,dp) - 9.38e-3_dp * Real(r8,dp)             &
  & + 1.53e-2_dp * Real(r7*Conjg(r10(2)) + r7p*Conjg(r10p(2)), dp )   &
  & + 6.848e-2_dp * Real(r7 * Conjg(r8) + r7p * Conjg(r8p), dp )      &
  & - 0.8545_dp * Real(r7*Conjg(r9(2)) + r7p*Conjg(r9p(2)), dp )      &
  & + 1.85e-3_dp * Real(r8*Conjg(r10(2)) + r8p*Conjg(r10p(2)), dp )   &
  & - 9.81e-2_dp * Real(r8*Conjg(r9(2)) + r8p*Conjg(r9p(2)), dp )     &
  & + 2.6917_dp * Real(r9(2),dp) + 0.2880_dp*(Abs(r7)**2+Abs(r7p)**2) &
  & - 0.10698_dp * Real( r9(2) * Conjg(r10(2))                        &
  &                    + r9p(2) * Conjg(r10p(2)), dp)                 &
  & + 10.7652_dp * (Abs(r10(2))**2 + Abs(r10p(2))**2 )                &
  & + 1.4884_dp * (Abs(r9(2))**2 + Abs(r9p(2))**2 )                   &
  & + 3.81e-3_dp * (Abs(r8)**2 + Abs(r8p)**2 ) )

  ! ratio BR(B -> Xs mu+ mu-)/BR(B -> Xs mu+ mu-)_SM
  ratioBtoSMuMu = BrBtoSMuMu/16.0479_dp

  ! branching ratio B -> Xs mu+ mu-
  BrBtoSMuMu = BrBtoSMuMu* 1.e-7_dp

Charged B-meson decay into charged Kaons and two charged leptons (B⁺ → K⁺ℓ⁺ℓ⁻)

Our results for B⁺ → K⁺ℓ⁺ℓ⁻ are based on the expressions given in . The branching ratio for B⁺ → K⁺μ⁺μ⁻ in the high-q² region, q² being the dilepton invariant mass squared, can be written as

\BR \left( B^+ \to K^+ \mu^+ \mu^- \right)_{q^2 \in \[14.18,22\] \text{GeV}^2} \simeq 1.11 + 0.22 \left( C_7^{\text{NP}} + C_7^\prime \right) + 0.27 \left( C_9^{\text{NP}} + C_9^\prime \right) - 0.27 \left( C_{10}^{\text{NP}} + C_{10}^\prime \right) \, .

The coefficients in Eq. can be related to the ones in our generic Lagrangian as

C_7^{\text{NP}} = n_{CQ} \, \left( Q_1^R - Q_1^{R,\text{SM}} \right) \\ C_7^\prime = n_{CQ} \, Q_1^L \\ C_9^{\text{NP}} = n_{CQ} \, \left[ \left( E_{LL}^V + E_{LR}^V \right) - \left( E_{LL}^{V,\text{SM}} + E_{LR}^{V,\text{SM}} \right) \right] \\ C_9^\prime = n_{CQ} \, \left( E_{RR}^V + E_{RL}^V \right) \\ C_{10}^{\text{NP}} = n_{CQ} \, \left[ \left( E_{LL}^V - E_{LR}^V \right) - \left( E_{LL}^{V,\text{SM}} - E_{LR}^{V,\text{SM}} \right) \right] \\ C_{10}^\prime = n_{CQ} \, \left( E_{RR}^V - E_{RL}^V \right)

where the normalization factor n_C**Q was already defined after Eq. .

NameProcess = "BtoKLL";
NameObservables = {{BrBtoKmumu, 6000, "BR(B -> K mu mu)"},
                   {ratioBtoKmumu, 6001, "BR(B -> K mu mu)/BR(B -> K mu mu)_SM"}};

NeededOperators = {OddllVRR, OddllVLL, OddllVRL, OddllVLR, CC7, CC7p,
                   OddllVRRSM, OddllVLLSM, OddllVRLSM, OddllVLRSM, CC7SM, CC7pSM
};

Body = "BtoKLL.f90";

Complex(dp) :: c7NP, c7p, c9NP, c9p, c10NP, c10p, norm
Real(dp) ::  GF

! ----------------------------------------------------------------
! B^+ -> K^+ l+ l- (14.18 GeV^2 < q^2 < 22 GeV^2)
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on W. Altmannshofer, D. M. Straub, EPJ C 73 (2013) 2646
! [arXiv:1308.1501]
! ----------------------------------------------------------------

c7NP = (CC7(3,2) - CC7SM(3,2))
c7p = CC7p(3,2)
c9NP = (OddllVLL(3,2,1,1)+OddllVLR(3,2,1,1) - &
            & (OddllVLLSM(3,2,1,1)+OddllVLRSM(3,2,1,1)))
c9p =  (OddllVRR(3,2,1,1)+OddllVRL(3,2,1,1))
c10NP = (OddllVLL(3,2,1,1)-OddllVLR(3,2,1,1) - &
            &  (OddllVLLSM(3,2,1,1)-OddllVLRSM(3,2,1,1)))
c10p = (OddllVRR(3,2,1,1)-OddllVRL(3,2,1,1))


! running GF
GF = (Alpha_160*4._dp*Pi/sinW2_160)/mw**2*sqrt2/8._dp

! normalization of our Wilson coefficients
! relative to the ones used in arXiv:1308.1501
norm = - oo16pi2*4._dp*GF/sqrt2*CKM_160(3,3)*Conjg(CKM_160(3,2))

! Branching ratio in the high-q^2 region
! q^2 in [14.18,22] GeV^2
BrBtoKmumu = (1.11_dp + 0.22_dp*(c7Np+c7p)/norm + &
   & 0.27_dp*(c9NP+c9p)/norm - 0.27_dp*(c10NP+c10p)/norm)

! ratio relative to SM
ratioBtoKmumu = BrBtoKmumu/1.11_dp

! total BR
BrBtoKmumu = BrBtoKmumu*1.0E-7_dp

B-meson decays into neutrino pairs (B̄ → X_d**,** sν****ν̄)

The branching ratio for B̄ → X_qν**ν̄, with q = d, s, is given by 

\BR \left( \bar B \to X_q \nu \bar \nu \right) = \frac{\alpha^2}{4 \pi^2 \sin^4 \theta_W} \frac{|V_{tb} V_{tq}^\ast|^2}{|V_{cb}|^2} \frac{\BR \left( \bar B \to X_c e \bar \nu_e \right) \kappa(0)}{f(\hat m_c) \kappa(\hat m_c)} \\ \times \sum_f \left[ \left( |c_L|^2 + |c_R|^2 \right) f(\hat m_q) - 4 \RE \left(c_L c_R^\ast \right) \hat m_q \tilde f(\hat m_q) \right] \nonumber \, .

The sum runs over the three neutrinos and m̂_i ≡ m_i/m_b. The functions f(m̂_c) and κ(m̂_c) represent the phase-space and the 1-loop QCD corrections, respectively. In case of κ(m̂_c), one needs the numerical values κ(0)=0.83 and κ(m̂_c)=0.88. The functions f(x) and f̃(x) take the form

f(x) = 1 - 8x^2 + 8 x^6-x^8 -24 x^4 \log x \\ \tilde f(x) = 1 + 9x^2 - 9x^4-x^6+12x^2(1+x^2)\log x \, .

Finally, \BR \left( \bar B \to X_c e \bar \nu_e \right)_{\text{exp}} = 0.101 and the coefficients c_L and c_R are given by

c_L = n^q_{BX \nu\nu} \, F_{LL}^V \\ c_R = n^q_{BX \nu\nu} \, F_{RL}^V \, ,

where \left(n^q_{BX \nu\nu}\right)^{-1} = \frac{4 G_F}{\sqrt{2}} \frac{\alpha}{2 \pi \sin^2 \theta_W} V_{tb}^\ast V_{tq} is the relative factor between our Wilson coefficients and the ones in .

NameProcess = "BtoQnunu";
NameObservables = {{BrBtoSnunu, 7000, "BR(B->s nu nu)"},
                   {ratioBtoSnunu, 7001, "BR(B->s nu nu)/BR(B->s nu nu)_SM"},
                   {BrBtoDnunu, 7002, "BR(B->D nu nu)"},
                   {ratioBtoDnunu, 7003, "BR(B->D nu nu)/BR(B->D nu nu)_SM"}};

NeededOperators = {OddvvVRR, OddvvVLL, OddvvVRL, OddvvVLR,
                   OddvvVRRSM, OddvvVLLSM, OddvvVRLSM, OddvvVLRSM};

Body = "BtoQnunu.f90";

Complex(dp) :: cL, cR, br, br_SM, cL_SM, cR_SM, norm
Real(dp) :: f_mq, tf_mq, kappa_0, kappa_c, f_mc, BrBXeNu, sw2, mq
Real(dp) :: prefactor, factor1, factor2, GF
Integer :: out, i1, i2

! ----------------------------------------------------------------
! \bar{B} -> X_{d,s} nu nu
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on C. Bobeth et al, NPB 630 (2002) 87 [hep-ph/0112305]
! ----------------------------------------------------------------

kappa_0 = 0.830_dp
kappa_c = 0.88_dp
f_mc = 0.53_dp
BrBXeNu = 0.101_dp ! PDG central value

sw2 = sinw2_160
GF = (Alpha_160*4._dp*Pi/sinW2_160)/mw**2*sqrt2/8._dp

Do out = 1,2
If (out.eq.1) Then ! B -> X_d nu nu
 mq = mf_d(1)/mf_d(3)
 norm = Alpha_160*4._dp*GF/sqrt2/(2._dp*pi*sinw2_160)* &
             & Conjg(CKM_160(3,3)*Conjg( CKM_160(3,1) ))
Else ! B -> X_s nu nu
 mq = mf_d(2)/mf_d(3)
 norm = Alpha_160*4._dp*GF/sqrt2/(2._dp*pi*sinw2_160)* &
             & Conjg(CKM_160(3,3)*Conjg( CKM_160(3,2) ))
End if

! f and tilde f functions
f_mq = 1._dp - 8._dp*mq**2 + 8._dp*mq**6 - &
            & mq**8 -24._dp*mq**4*Log(mq)
tf_mq = 1._dp + 9._dp*mq**2 - 9._dp*mq**4 - mq**6 + &
            & 12._dp*mq**2*(1._dp + mq**2)*Log(mq)

prefactor =  Alpha_mz**2/(4._dp*pi**2*sw2**2)*Abs(CKM_160(3,3)/ &
               & CKM_160(2,3))**2*BrBXeNu/(f_mc*kappa_c)*kappa_0
factor1 = f_mq
factor2 = - 4._dp*mq*tf_mq

br = 0._dp
br_SM = 0._dp

 Do i1= 1,3
  Do i2 = 1,3

   ! BSM
   cL = OddvvVLL(3,out,i1,i2)/norm
   cR = OddvvVRL(3,out,i1,i2)/norm
   br = br + factor1*(Abs(cL)**2 + Abs(cR)**2) +    &
               &  factor2*Real(cL*Conjg(cR),dp)

   ! SM
   cL = OddvvVLLSM(3,out,i1,i2)/norm
   cR = OddvvVRLSM(3,out,i1,i2)/norm
   br_SM = br_SM + factor1*(Abs(cL)**2 + Abs(cR)**2) + &
            & factor2*Real(cL*Conjg(cR),dp)

  End Do
 End do
If (out.eq.1) Then ! B -> X_d nu nu
  BrBtoDnunu = prefactor*br*Abs(CKM_160(3,1))**2
  ratioBtoDnunu = br/br_SM
Else ! B -> X_s nu nu
  BrBtoSnunu = prefactor*br*Abs(CKM_160(3,2))**2
  ratioBtoSnunu = br/br_SM
End if
End Do

Kaon decays into Pions and pairs of neutrinos (K → πνν̄)

Following , the branching ratios for rare Kaon decays involving neutrinos in the final state can be written as

\BR \left( K^+ \to \pi^+ \nu \bar \nu \right) = 2 r_1 \, r_2 \, r_{K^+} \sum_f \left[ \left( \text{Im} \lambda_t X_f \right)^2 + \left( \text{Re} \lambda_c X_{NL} + \text{Re} \lambda_t X_f \right)^2 \right] \\ \BR \left( K_L \to \pi^0 \nu \bar \nu \right) = 2 r_1 \, r_{K_L} \sum_f \left( \text{Im} \lambda_t X_f \right)^2 \, ,

where the sums are over the three neutrino species, X_N**L = 9.78 ⋅ 10⁻⁴ is the SM NLO charm correction , λ_t = V_t**s^*V_t**d and λ_c = V_c**s^*V_c**d, the coefficients r₁, r₂, r_K⁺ and r_KL take the numerical values

r_1 = 1.17 \cdot 10^{-4} \nonumber \\ r_2 = 0.24 \nonumber \\ r_{K^+} = 0.901 \nonumber \\ r_{K_L} = 0.944

and X_f contains the Wilson coefficients contributing to the processes, F_L**L^V and F_R**L^V, as

X_f = n_Kπν**ν (F_L**L^V+F_R**L^V) .

Here n_{K \pi \nu\nu}^{-1} = \frac{4 G_F}{\sqrt{2}} \frac{\alpha}{2 \pi \sin^2 \theta_W} V_{ts}^\ast V_{td}.

NameProcess = "KtoPInunu";
NameObservables = {{BrKptoPipnunu, 8000, "BR(K^+ -> pi^+ nu nu)"},
                   {ratioKptoPipnunu, 8001, "BR(K^+ -> pi^+ nu nu)/BR(K^+ -> pi^+ nu nu)_SM"},
                   {BrKltoPinunu, 8002, "BR(K_L -> pi^0 nu nu)"},
                   {ratioKltoPinunu, 8003, "BR(K_L -> pi^0 nu nu)/BR(K_L -> pi^0 nu nu)_SM"}};

NeededOperators = {OddvvVRR, OddvvVLL, OddvvVRL, OddvvVLR,
                   OddvvVRRSM, OddvvVLLSM, OddvvVRLSM, OddvvVLRSM};

Body = "KtoPInunu.f90";

Complex(dp) :: br, r1, r2, rKp, rKl, Xx, XNL, Lt, Lc
Complex(dp) :: Xx_SM, br_SM, norm
Real(dp) :: GF
Integer :: out, i1, i2

! ----------------------------------------------------------------
! K -> pi nu nu
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on C. Bobeth et al, NPB 630 (2002) 87 [hep-ph/0112305]
! ----------------------------------------------------------------

GF = (Alpha_160*4._dp*Pi/sinW2_160)/mw**2*sqrt2/8._dp
norm = Alpha_160*4._dp*GF/sqrt2/(2._dp*pi*sinw2_160) &
          & *Conjg(CKM_160(3,2))*CKM_160(3,1)

r1 = 1.17E-4_dp
r2 = 0.24_dp
rKp =  0.901
rKl = 0.944

! SM NLO charm correction
! See G. Buchalla and A. Buras, NPB 412 (1994) 106 and NPB 548 (1999) 309
XNL = 9.78E-4_dp

! out = 1 : K^+ -> pi^+ nu nu
! out = 2 : K_L -> pi^0 nu nu

Do out = 1,2
br = 0._dp
br_SM = 0._dp
 Do i1= 1,3
  Do i2 = 1,3
   Xx = ((OddvvVLL(2,1,i1,i2)+OddvvVRL(2,1,i1,i2))/norm)
   Xx_SM = ((OddvvVLLSM(2,1,i1,i2)+OddvvVRLSM(2,1,i1,i2))/norm)
   Lt = Conjg(CKM_160(3,2))*CKM_160(3,1)
   Lc = Conjg(CKM_160(2,2))*CKM_160(2,1)
   If (out.eq.1) Then
     br = br + Aimag(Xx*Lt)**2 + (Real(Lc*XNL,dp) + Real(Xx*Lt,dp))**2
     br_SM = br_SM + Aimag(Xx_SM*Lt)**2 + &
                 & (Real(Lc*XNL,dp) + Real(Xx_SM*Lt,dp))**2
   Else
    br = br + Abs(Aimag(Xx*Lt))**2
    br_SM = br_SM + Abs(Aimag(Xx_SM*Lt))**2
   End if
  End Do
 End do
If (out.eq.1) Then ! K^+ -> pi^+ nu nu
 BrKptoPipnunu = 2._dp*r1*r2*rKp*br
 RatioKptoPipnunu = br/br_SM
 ! SM expectation: (7.2 +/- 2.1)*10^-11 (hep-ph/0112135)
Else ! K_L -> pi^0 nu nu
 BrKltoPinunu = 2._dp*r1*rKl*br
 RatioKltoPinunu = br/br_SM
 ! SM expectation: (3.1 +/- 1.0)*10^-11 (hep-ph/0408142)
End if
End Do

Mass splitting of neutral B-mesons (Δ****M_Bs**,** d)

The B_q⁰ − B̄_q⁰ mass difference can be written as

\Delta M_{B_q} = \frac{G_F^2 m_W^2}{6 \pi^2} m_{B_q} \eta_B f_{B_q}^2 \hat B_{B_q} |V_{tq}^{\text{eff}}|^2 |F_{tt}^q| \, ,

where q = s, d, m_Bq and f_Bq are the B_q⁰ mass and decay constant, respectively, η_B = 0.55 is a QCD factor , B̂_Bq is a non-perturbative parameter (with values B̂_Bd = 1.26 and B̂_Bs = 1.33, obtained from recent lattice computations ) and |V_t**q^eff|² = (V_t**b^*V_t**q)². F_t**t^q is given by

F_{tt}^q = S_0(x_t) + \frac{1}{4 r} C_{\text{new}}^{VLL} \\ + \frac{1}{4 r} C_1^{VRR} + \bar{P}_1^{LR} C_1^{LR} + \bar{P}_2^{LR} C_2^{LR} \nonumber \\ + \bar{P}_1^{SLL} \left( C_1^{SLL} + C_1^{SRR} \right) + \bar{P}_2^{SLL} \left( C_2^{SLL} + C_2^{SRR} \right) \nonumber

where r = 0.985 , x_t = \frac{m_t^2}{m_W^2}, with m_t the top quark mass, the P̄ coefficients take the numerical values

\bar{P}_1^{LR} = -0.71 \nonumber \\ \bar{P}_2^{LR} = 0.90 \nonumber \\ \bar{P}_1^{SLL} = -0.37 \nonumber \\ \bar{P}_2^{SLL} = -0.72

and the function

S_0(x_t) = \frac{4 x_t - 11 x_t^2 + x_t^3}{4 (1-x_t)^2} - \frac{3 x_t^3 \log x_t}{2(1-x_t)^3}

was introduced by Inami and Lim in and given, for example, in . Finally, the coefficients in Eq. are related to the D_X**Y^I coefficients in Eq. as

C_{\text{new}}^{VLL} = n^q_\Delta \left( D_{LL}^V - D_{LL}^{V,\text{SM}} \right) \\ C_1^{VRR} = n^q_\Delta D_{RR}^V \\ C_1^{LR} = n^q_\Delta \left( D_{LR}^V + D_{RL}^V \right) \\ C_2^{LR} = n^q_\Delta \left( D_{LR}^S + D_{RL}^S + \delta_2^{LR} \right) \\ C_1^{SLL} = n^q_\Delta \left( D_{LL}^S + \delta_1^{SLL} \right) \\ C_1^{SRR} = n^q_\Delta \left( D_{RR}^S + \delta_1^{SRR} \right) \\ C_2^{SLL} = n^q_\Delta D_{LL}^T \\ C_2^{SRR} = n^q_\Delta D_{RR}^T

where the factor \left(n^q_\Delta\right)^{-1} = \frac{G_F^2 m_W^2}{16 \pi^2} |V_{tq}^{\text{eff}}|^2 normalizes our Wilson coefficients to the ones in . The corrections δ₂^L**R, δ₁^SLL and δ₁^SRR are induced by double penguin diagrams mediated by scalar and pseudoscalar states . These 2-loop contributions may have a sizable impact in some models, and their inclusion is necessary in order to achieve a precise result for Δ**M_Bq. They can be written as

\delta_2^{LR} = - \frac{H_L^{S,P} \left(H_R^{S,P}\right)^\ast}{m_{S,P}^2} \\ \delta_1^{SLL} = - \frac{\left(H_L^{S,P}\right)^2}{2 \, m_{S,P}^2} \\ \delta_1^{SRR} = - \frac{\left(H_L^{S,P}\right)^2}{2 \, m_{S,P}^2}

where H_L^S, P and H_R^S, P are defined in Eq.. The double penguin corrections in Eqs.- are obtained by summing up over all scalar and pseudoscalar states in the model.

NameProcess = "DeltaMBq";
NameObservables = {{DeltaMBs, 1900, "Delta(M_Bs)"},
                   {ratioDeltaMBs, 1901, "Delta(M_Bs)/Delta(M_Bs)_SM"},
                   {DeltaMBq, 1902, "Delta(M_Bd)"},
                   {ratioDeltaMBq, 1903, "Delta(M_Bd)/Delta(M_Bd)_SM"}};

ExternalStates =  {Fd};
NeededOperators = {O4dSLL, O4dSRR, O4dSRL, O4dSLR, O4dVRR, O4dVLL,
                  O4dVLLSM, O4dVRL, O4dVLR, O4dTLL, O4dTLR, O4dTRL, O4dTRR};

IncludeSMprediction["DeltaMBq"] = False;

Body = "DeltaMBq.f90";

Complex(dp) :: MBq, etaB, FBq2, BBq, Ftt, Veff2, r, &
     & P1bLR, P2bLR, P1bSLL, P2bSLL, norm, &
     & CVLLnew, C1VRR, C1LR, C2LR, C1SLL, C1SRR, C2SLL, C2SRR
Real(dp) ::  hbar, xt, GF
Real(dp) :: mS
Complex(dp) :: HL, HR, AL, AR
Integer :: i1, iS

! ----------------------------------------------------------------
! Delta M_{Bd,Bs}
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on A. J. Buras et al, NPB 619 (2001) 434 [hep-ph/0107048]
! and NPB 659 (2003) 3 [hep-ph/0210145]
! ----------------------------------------------------------------

hbar =  6.58211889e-25_dp
xt = mf_u2_160(3)/mw2
r = 0.985_dp
P1bLR = -0.71_dp
P2bLR = 0.90_dp
P1bSLL = -0.37_dp
P2bSLL = -0.72_dp

! QCD factor, see A. J. Buras et al, NPB 47 (1990) 491
! and J. Urban et al, NPB 523 (1998) 40
etaB = 0.55_dp

GF =  (Alpha_160*4._dp*Pi/sinW2_160)/mw**2*sqrt2/8._dp

Do i1 = 1,2

If (i1.eq.1) Then ! Delta M_Bd
 MBq = mass_B0d
 FBq2 = f_B0d_CONST**2
 BBq = 1.26_dp ! see arXiv:0910.2928
 Veff2 = Conjg(Conjg(CKM_160(3,3))*CKM_160(3,1))**2
Else ! Delta M_Bs
 MBq = mass_B0s
 FBq2 = f_B0s_CONST**2
 BBq = 1.33_dp ! see arXiv:0910.2928
 Veff2 = Conjg(Conjg(CKM_160(3,3))*CKM_160(3,2))**2
End if

! normalization factor
norm = GF**2*mw2/(16._dp*Pi**2)*Veff2

! Wilson coefficients
CVLLnew = (O4dVLL(3,i1,3,i1)-O4dVLLSM(3,i1,3,i1))/norm ! we remove the SM contribution
C1VRR = O4dVRR(3,i1,3,i1)/norm
C1LR = (O4dVLR(3,i1,3,i1)+O4dVRL(3,i1,3,i1))/norm
C2LR = (O4dSLR(3,i1,3,i1)+O4dSRL(3,i1,3,i1))/norm
C1SLL = O4dSLL(3,i1,3,i1)/norm
C1SRR = O4dSRR(3,i1,3,i1)/norm
C2SLL  = O4dTLL(3,i1,3,i1)/norm
C2SRR  = O4dTRR(3,i1,3,i1)/norm


! Double Higgs penguins
@ If[getGen[HiggsBoson] > 1, "Do iS = 1, "<>ToString[getGen[HiggsBoson]],""]
@ If[getGen[HiggsBoson] > 1, "HL = OH2qSL(3,i1,iS)",  "HL = OH2qSL(3,i1)"]
@ If[getGen[HiggsBoson] > 1, "HR = OH2qSR(3,i1,iS)",  "HR = OH2qSR(3,i1)"]
@ If[getGen[HiggsBoson] > 1, "mS = "<>SPhenoMassSq[HiggsBoson,iS],  "mS = "<>SPhenoMassSq[HiggsBoson]]
C2LR = C2LR - HL*Conjg(HR)/(mS*norm)
C1SLL = C1SLL - 0.5_dp*HL**2/(mS*norm)
C1SRR = C1SRR - 0.5_dp*HR**2/(mS*norm)
@ If[getGen[HiggsBoson] > 1,"End Do",""]


@ If[getGen[PseudoScalar] > 1, "Do iS = "<>ToString[getGenSPhenoStart[PseudoScalar]]<>", "<>ToString[getGen[PseudoScalar]],""]
@ If[getGen[PseudoScalar] > 1, "AL = OAh2qSL(3,i1,iS)",  "AL = OAh2qSL(3,i1)"]
@ If[getGen[PseudoScalar] > 1, "AR = OAh2qSR(3,i1,iS)",  "AR = OAh2qSR(3,i1)"]
@ If[getGen[PseudoScalar] > 1, "mS = "<>SPhenoMassSq[PseudoScalar,iS],  "mS = "<>SPhenoMassSq[PseudoScalar]]
C2LR = C2LR - AL*Conjg(AR)/(mS*norm)
C1SLL = C1SLL - 0.5_dp*AL**2/(mS*norm)
C1SRR = C1SRR - 0.5_dp*AR**2/(mS*norm)
@ If[getGen[PseudoScalar] > 1,"End Do",""]


Ftt = S0xt(xt) + CVLLnew/(4._dp*r) + &
     & C1VRR/(4._dp*r) + P1bLR*C1LR + P2bLR*C2LR  + &
     & P1bSLL*(C1SLL + C1SRR) + P2bSLL*(C2SLL + C2SRR)

If (i1.eq.1) Then ! Delta M_Bd
   ratioDeltaMBq = Abs(Ftt/S0xt(xt))
   DeltaMBq = G_F**2*mw2/(6._dp*Pi**2)*    &
      & MBq*etaB*BBq*FBq2*Veff2*Abs(Ftt)*1.e-12_dp/hbar
Else ! Delta M_Bs
   ratioDeltaMBs = Abs(Ftt/S0xt(xt))
   DeltaMBs =  G_F**2*mw2/(6._dp*Pi**2)*   &
      & MBq*etaB*BBq*FBq2*Veff2*Abs(Ftt)*1.e-12_dp/hbar
End if

End Do

Contains

  Real(dp) Function S0xt(x) ! See for example hep-ph/9806471
    Implicit None
    Real(dp), Intent(in) :: x
    S0xt = 1._dp - 2.75_dp * x + 0.25_dp * x**2 &
        & - 1.5_dp * x**2 * Log(x) / (1-x)
    S0xt = x*S0xt / (1 -x)**2
  End  Function S0xt

Mass splitting of neutral Koans (Δ****M_K, ε_K)

Δ**M_K and ε_K, the observables associated to K⁰ − K̄⁰ mixing, can be written as 

\Delta M_{K} = 2 \, \text{Re} \, \langle \bar K^0 | H_\text{eff}^{\Delta S = 2} | K^0 \rangle \\ \varepsilon_K = \frac{e^{i \pi/4}}{\sqrt{2} \Delta M_{K}} \, \text{Im} \, \langle \bar K^0 | H_\text{eff}^{\Delta S = 2} | K^0 \rangle \, .

The matrix element in Eqs. and is given by

\langle \bar K^0 | H_\text{eff}^{\Delta S = 2} | K^0 \rangle = f_V \left( D_{LL}^V + D_{RR}^V \right) + f_S \left( D_{LL}^S + D_{RR}^S \right) + f_T \left( D_{LL}^T + D_{RR}^T \right) \nonumber \\ + f_{LR}^1 \left( D_{LR}^S + D_{RL}^S \right) + f_{LR}^2 \left( D_{LR}^V + D_{RL}^V \right) \, .

The f coefficients are

f_V = \frac{1}{3} m_K f_K^2 B_1^{VLL}(\mu) \\ f_S = - \frac{5}{24} \left( \frac{m_K}{m_s(\mu) + m_d(\mu)} \right)^2 m_K f_K^2 B_1^{SLL}(\mu) \\ f_T = - \frac{1}{2} \left( \frac{m_K}{m_s(\mu) + m_d(\mu)} \right)^2 m_K f_K^2 B_2^{SLL}(\mu) \\ f_{LR}^1 = - \frac{1}{6} \left( \frac{m_K}{m_s(\mu) + m_d(\mu)} \right)^2 m_K f_K^2 B_1^{LR}(\mu) \\ f_{LR}^2 = \frac{1}{4} \left( \frac{m_K}{m_s(\mu) + m_d(\mu)} \right)^2 m_K f_K^2 B_2^{LR}(\mu)

where μ = 2 GeV is the energy scale at which the matrix element is computed and f_K the Kaon decay constant. The values of the quark masses at μ = 2 GeV are given by m_d(μ)=7 MeV and m_s(μ)=125 MeV (see table 1 in ), whereas the B_i^X coefficients have the following values at μ = 2 GeV : B₁^VLL(μ)=0.61, B₁^SLL(μ)=0.76, B₂^SLL(μ)=0.51, B₁^L**R(μ)=0.96 and B₂^L**R(μ)=1.3.

As in , we treat the SM contribution separately. We define D_L**L^V = D_L**L^V, S**M + D_L**L^V, BSM. For D_L**L^V, BSM one just subtracts the SM contributions to D_L**L^V, whereas for D_L**L^V, S**M one can use the results in , where the relevant QCD corrections are included,

D_{LL}^{V,SM} = \frac{G_F^2 m_W^2}{4 \pi^2} \left[ \lambda_c^{\ast \, 2} \eta_1 S_0(x_c) + \lambda_t^{\ast \, 2} \eta_2 S_0(x_t) + 2 \lambda_c^\ast \lambda_t^\ast \eta_3 S_0(x_c,x_t) \right] \, .

Here x_i = m_i²/m_w², λ_i = V_i**s^*V_i**d and S₀(x) and S₀(x, y) are the Inami-Lim functions . S₀(x) was already defined in Eq. , whereas S₀(x_c, x_t) is given by 

S_0(x_c,x_t) = x_c \left[ \log \frac{x_t}{x_c} - \frac{3 x_t}{4(1-x_t)} - \frac{3 x_t^2 \log x_t}{4 (1-x_t)^2} \right] \, .

In the last expression we have kept only terms linear in x_c ≪ 1. Finally, the η_i coefficients comprise short distance QCD corrections. Their numerical values are η_1, 2, 3 = (1.44,0.57,0.47)  [2].

NameProcess = "KKmix";
NameObservables = {{DeltaMK, 9100, "Delta(M_K)"},
                   {ratioDeltaMK, 9102, "Delta(M_K)/Delta(M_K)_SM"},
                   {epsK, 9103, "epsilon_K"},
                   {ratioepsK, 9104, "epsilon_K/epsilon_K^SM"}};

NeededOperators = {O4dSLL, O4dSRR, O4dSRL, O4dSLR, O4dVRR, O4dVLL, O4dVRL,
                   O4dVLR, O4dTLL, O4dTLR, O4dTRL, O4dTRR,
                   O4dSLLSM, O4dSRRSM, O4dSRLSM, O4dSLRSM, O4dVRRSM, O4dVLLSM, O4dVRLSM, O4dVLRSM,
                   O4dTLLSM, O4dTLRSM, O4dTRLSM, O4dTRRSM};

Body = "KKmix.f90";

Real(dp) :: b_VLL, b_SLL1, b_SLL2, b_LR1, b_LR2
Real(dp) :: ms_mu, md_mu
Complex(dp) :: CVLL, CVRR, CSLL, CSRR, CTLL, CTRR, CLR1, CLR2
Complex(dp) :: fV, fS, fT, fLR1, fLR2, cVLLSM
Complex(dp) :: f_k, M_K, H2eff, DeltaMK_SM, epsK_SM
Real(dp) :: norm, hbar, xt, xc, GF
Integer :: i1
Real(dp), Parameter :: eta_tt = 0.57_dp, eta_ct = 0.47_dp, &
                         & eta_cc = 1.44_dp
! Parameters from S. Herrlich and U. Nierste NPB 476 (1996) 27

! ----------------------------------------------------------------
! Delta M_K and epsilon_K
! Observables implemented by W. Porod, F. Staub and A. Vicente
! Based on A. Crivellin et al, Comput. Phys. Commun. 184 (2013) 1004 [arXiv:1203.5023]
! ----------------------------------------------------------------

! using globally defined hadronic parameters
M_K = mass_K0
f_K = f_k_CONST

xt = mf_u(3)**2 / mW**2
xc = mf_u(2)**2 / mW**2

GF = (Alpha_160*4._dp*Pi/sinW2_160)/mw**2*sqrt2/8._dp

! Coefficients at mu = 2 GeV
! See A. J. Buras et al, NPB 605 (2001) 600 [hep-ph/0102316]
b_VLL = 0.61_dp
b_SLL1 = 0.76_dp
b_SLL2 = 0.51_dp
b_LR1 = 0.96_dp
b_LR2 = 1.3_dp

! Quark mass values at mu = 2 GeV
! See M. Ciuchini et al, JHEP 9810 (1998) 008 [hep-ph/9808328] - Table 1
md_mu = 0.007_dp
ms_mu = 0.125_dp

fV = 1._dp/3._dp*M_K*f_k**2*b_VLL
fS = -5._dp/24._dp*M_K*f_K**2*(M_K/(ms_mu+md_mu))**2*b_SLL1
fT = -1._dp/2._dp*M_K*f_K**2*(M_K/(ms_mu+md_mu))**2*b_SLL2
fLR1 = -1._dp/6._dp*M_K*f_K**2*(M_K/(ms_mu+md_mu))**2*b_LR1
fLR2 = 1._dp/4._dp*M_K*f_K**2*(M_K/(ms_mu+md_mu))**2*b_LR2

! SM contribution
! Based on the results by S. Herrlich and U. Nierste
! NPB 419 (1994) 292, PRD 52 (1995) 6505 and NPB 476 (1996) 27
cVLLSM = eta_cc * (Conjg(CKM_160(2,2))*CKM_160(2,1))**2 * S0xt(xc)    &
     & + eta_tt * (Conjg(CKM_160(3,2))*CKM_160(3,1))**2 * S0xt(xt)    &
     & + Conjg(CKM_160(2,2)*CKM_160(3,2))*(CKM_160(2,1)*CKM_160(3,1)) &
     &   * 2._dp * eta_ct * S0_2(xc,xt)

cVLLSM = Conjg(cVLLSM)  ! we compute (d\bar{s})(d\bar{s}) and not (\bar{d}s)(\bar{d}s)
cVLLSM = oo4pi2*(GF*mW)**2*cVLLSM ! normalization

! BSM contributions (+SM in CVLL)
CVLL = O4dVLL(2,1,2,1)-O4dVLLSM(2,1,2,1)+cVLLSM
CVRR = O4dVRR(2,1,2,1)
CSLL = O4dSLL(2,1,2,1)
CSRR = O4dSRR(2,1,2,1)
CTLL = O4dTLL(2,1,2,1)
CTRR = O4dTRR(2,1,2,1)
CLR1 = O4dSLR(2,1,2,1)+O4dSRL(2,1,2,1)
CLR2 = O4dVLR(2,1,2,1)+O4dVRL(2,1,2,1)

! BSM
H2eff = fV*(CVLL+CVRR) + fS*(CSLL+CSRR) +fT*(CTLL+CTRR) &
              & + fLR1*CLR1 + fLR2*CLR2

DeltaMK = Abs(2._dp*Real(H2eff,dp))
epsK = 1._dp/(sqrt2*DeltaMK)*Abs(Aimag(H2eff))

! SM
H2eff = fV*cVLLSM

DeltaMK_SM = Abs(2._dp*Real(H2eff,dp))
epsK_SM = 1._dp/(sqrt2*DeltaMK_SM)*Abs(Aimag(H2eff))

ratioDeltaMK = DeltaMK/DeltaMK_SM
ratioepsK = epsK/epsK_SM

Contains

! Inami - Lim functions

 Real(dp) Function S0xt(x)
   Implicit None
   Real(dp), Intent(in) :: x
   S0xt = 1._dp - 2.75_dp * x + 0.25_dp * x**2 - &
              &  1.5_dp * x**2 * Log(x) / (1-x)
   S0xt = x*S0xt / (1 -x)**2
 End  Function S0xt

 Real(dp) Function S0_2(xc, xt)
   Implicit None
   Real(dp), Intent(in) :: xc, xt
   S0_2 = Log(xt/xc) - 0.75_dp *  xt /(1-xt) &
        & - 0.75_dp * xt**2 * Log(xt) / (1-xt)**2
   S0_2 = xc *  S0_2
 End  Function S0_2

Pseudo-scalar decays into charged leptons and neutrinos (P → ℓ****ν)

Although P → ℓν, where P = q**q′ is a pseudoscalar meson, does not violate quark flavor, we have included it in the list of observables for practical reasons, as it can be computed with the same ingredients as the QFV observables. The decay width for the process P → ℓ_αν is given by

\Gamma \left( P \to \ell_\alpha \nu \right) = \frac{|G_F f_P (m_P^2 - m_{\ell_\alpha}^2)|^2}{8 \pi m_P^3} \\ \times \sum_\nu \left| V_{qq'} m_{\ell_\alpha} + \frac{m_{\ell_\alpha}}{2 \sqrt{2}} \left( G_{LL}^V - G_{RL}^V \right) + \frac{m_P^2}{2 \sqrt{2} (m_q + m_{q'})} \left( G_{RR}^S - G_{LR}^S \right) \right|^2 \, . \nonumber

Here f_P is the meson decay constant, m_q and m_q′ are the masses of the quarks in the meson and the Wilson coefficients G_X**Y^I are defined in Eq.. The sum in Eq. is over the three neutrinos (whose masses are neglected).

Each P → ℓ_αν decay width is plagued by hadronic uncertainties. However, by taking the ratios

R_P = \frac{\Gamma \left( P \to e \nu \right)}{\Gamma \left( P \to \mu \nu \right)}

the hadronic uncertainties cancel out to a good approximation, allowing for a precise theoretical determination. In case of R_K, the SM prediction includes small electromagnetic corrections that account for internal bremsstrahlung and structure-dependent effects . This leads to an impressive theoretical uncertainty of δ**R_K/R_K ∼ 0.1%, making R_P the perfect observable to search for lepton flavor universality violation .

NameProcess = "Plnu";
NameObservables = {{BrDmunu, 300, "BR(D->mu nu)"},
                   {ratioDmunu, 301, "BR(D->mu nu)/BR(D->mu nu)_SM"},
                   {BrDsmunu, 400, "BR(Ds->mu nu)"},
                   {ratioDsmunu, 401, "BR(Ds->mu nu)/BR(Ds->mu nu)_SM"},
                   {BrDstaunu, 402, "BR(Ds->tau nu)"},
                   {ratioDstaunu, 403, "BR(Ds->tau nu)/BR(Ds->tau nu)_SM"},
                   {BrBmunu, 500, "BR(B->mu nu)"},
                   {ratioBmunu, 501, "BR(B->mu nu)/BR(B->mu nu)_SM"},
                   {BrBtaunu, 502, "BR(B->tau nu)"},
                   {ratioBtaunu, 503, "BR(B->tau nu)/BR(B->tau nu)_SM"},
                   {BrKmunu, 600, "BR(K->mu nu)"},
                   {ratioKmunu, 601, "BR(K->mu nu)/BR(K->mu nu)_SM"},
                   {RK, 602 ,"R_K = BR(K->e nu)/(K->mu nu)"},
                   {RKSM, 603 ,"R_K^SM = BR(K->e nu)_SM/(K->mu nu)_SM"}};

NeededOperators = {OdulvSLL, OdulvSRR, OdulvSRL, OdulvSLR,
                   OdulvVRR, OdulvVLL, OdulvVRL, OdulvVLR,
                   OdulvSLLSM, OdulvSRRSM, OdulvSRLSM, OdulvSLRSM,
                   OdulvVRRSM, OdulvVLLSM, OdulvVRLSM, OdulvVLRSM
};

Body = "Plnu.f90";

Integer :: gt1, gt2, i1, i2, iP
Complex(dp) :: br, br_SM
Real(dp) :: m_M, f_M, tau_M, mlep, mq1, mq2, hbar, ratio, &
     & BrKenuSM, BRKenu, QED

! ----------------------------------------------------------------
! P -> l nu
! Observable implemented by W. Porod, F. Staub and A. Vicente
! Based on J. Barranco et al, arXiv:1303.3896
! ----------------------------------------------------------------

hbar = 6.58211889e-25_dp

! Electromagnetic correction to R_K
! See V. Cirigliano, I. Rosell, PRL 99 (2007) 231801 [arXiv:0707.3439]
QED = -3.6e-2_dp

! meson parameters

Do iP=1,4
If (iP.eq.1) Then ! Ds-meson
 gt1 = 2
 gt2 = 2
 m_M = mass_Dsp
 f_M = f_DSp_CONST
 tau_M = tau_DSp/hbar
Elseif (iP.eq.2) Then ! B-meson
 gt1 = 3
 gt2 = 1
 m_M = mass_Bp
 f_M = f_Bp_CONST
 tau_M = tau_Bp/hbar
Elseif (iP.eq.3) Then ! Kaon
 gt1 = 2
 gt2 = 1
 m_M = mass_Kp
 f_M =  f_Kp_CONST
 tau_M = tau_Kp/hbar
Elseif (iP.eq.4) Then ! D-meson
 gt1 = 1
 gt2 = 2
 m_M = mass_Dp
 f_M =  f_Dp_CONST
 tau_M = tau_Dp/hbar
End if

 mq1 =  mf_u_160(gt2)
 mq2 =  mf_d_160(gt1)

Do i1=1,3
br = 0._dp
br_SM = 0._dp
mlep = mf_l(i1)

Do i2=1,3
 br = br + ((OdulvVLL(gt1,gt2,i1,i2)-OdulvVLR(gt1,gt2,i1,i2))*mlep/   &
                   &           (2._dp*sqrt2)                          &
   & + m_M**2*(OdulvSRL(gt1,gt2,i1,i2)-OdulvSLL(gt1,gt2,i1,i2))/      &
                   & (2._dp*sqrt2*(mq1+mq2)))
 br_SM = br_SM+ (OdulvVLLSM(gt1,gt2,i1,i2)-OdulvVLRSM(gt1,gt2,i1,i2)) &
                   & *mlep/(2._dp*sqrt2)
End Do

ratio = Abs(br/br_SM)**2
br = oo8pi*tau_M*(f_M)**2*M_M*Abs(br)**2*(1._dp - mlep**2/M_M**2)**2 ! G_F already in coefficients included


If (iP.eq.1) Then  !! Ds-meson
 If (i1.eq.2) Then  ! Ds->mu nu
  BrDsmunu =   br
  ratioDsmunu = ratio
  Elseif (i1.eq.3) Then ! Ds->tau nu
  BrDstaunu =  br
  ratioDstaunu = ratio
 End if
Elseif (iP.eq.2) Then !! B-meson
 If (i1.eq.2) Then  ! B->mu nu
  BrBmunu =   br
  ratioBmunu = ratio
  Else              ! B->tau nu
  BrBtaunu =  br
  ratioBtaunu = ratio
 End if
Else If (iP.eq.3) Then !! Kaon
 If (i1.eq.1) Then  ! K->e nu
  BrKenu =   br
  BrKenuSM = BrKenu*ratio
  Elseif (i1.eq.2) Then  ! K->mu nu
  BrKmunu =  br
  ratioKmunu = ratio
  RK = BrKenu/BrKmunu*(1+QED)
  RKSM = BrKenuSM/BrKmunu*ratio*(1+QED)
 End if
Else If (iP.eq.4) Then  !! D-meson
 If (i1.eq.2) Then  ! D->mu nu
  BrDmunu =   br
  ratioDmunu = ratio
 End if
End if
End Do
End Do

See also

[1] Notice that our effective Lagrangian differs from the one in by a 1/(4π)² factor. This relative factor has been absorbed in the expression for ℳ_ℬℓℓ, see Eq..

[2] Note that we have chosen a value for η₁ which results from our numerical values for α_s(m_Z) and m_c(m_c), see table 5 in .