1 Introduction

The incredible joint theoretical and experimental effort carried out in the last years allows us to probe the Standard Model (SM) of particle physics with an unprecedented precision. This brought to light some deviations between theoretical predictions and experimental measurements in semileptonic B meson decays [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. These discrepancies, the so-called B anomalies, hint at Lepton Flavour Universality (LFU) violation. This is quite surprising, as LFU is one of the foundation of the SM.

The B anomalies can be split into two classes: (i) deviations in \(\mu /e\) universality in \(b\rightarrow s\ell ^+\ell ^-\) and (ii) deviations in \(\tau \) vs. light leptons universality in the \(b\rightarrow c\ell {\bar{\nu }}\) transitions. These exciting findings may indicate the presence of New Physics (NP) particles, and have inspired a plethora of theoretical and experimental work. The NP explanations for B anomalies span a broad class of new heavy particles, from vectors to scalars states [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60]. However, a common feature in all of these models is the prediction of sizeable effects for Lepton Flavour Violating (LFV) B, \(\tau \) and \(\mu \) decays. Current upper bounds on these modes largely constrain the allowed parameter space for NP models, and upcoming experimental analyses will be fundamental to corroborate or falsify these NP hypotheses.

When searching for LFV decays mediated by \(b\rightarrow s\ell _1\ell _2\) transitions, it is crucial to consider both mesonic and baryonic decays. Although mediated by the same underlying partonic transition, these two types of decays provide orthogonal information on possible NP models. A striking example of this is the NP analysis of \(\Lambda _b\rightarrow \Lambda \mu ^+\mu ^-\) decays [61], which shows that even though \(B\rightarrow K\mu ^+\mu ^-\) angular distribution seems to be affected by some short-distance NP [62,63,64,65], the latter is not visible in the \(\Lambda _b\rightarrow \Lambda \mu ^+\mu ^-\) angular observables. It is therefore natural to assume that \(\Lambda _b\rightarrow \Lambda \ell _1^-\ell _2^+\) decays provide complementary information compared to their mesonic counterparts \(B^+\rightarrow K^+\ell _1\ell _2\) or \({{\bar{B}}}_s\rightarrow \ell _1\ell _2\) decays. More precisely, the spin structure of the \(\Lambda _b\rightarrow \Lambda \) decays induces a richer set of hadronic matrix elements. Therefore, \(\Lambda _b\rightarrow \Lambda \) decays probe different parameter space than their mesonic counter parts. In addition, the \(\Lambda _b\) baryon is copiously produced at LHCb and the dataset collected with Run 1 and Run 2 allow to make precision measurement of observables constructed from \(\Lambda _b\) decays (see e.g. [66]).

In this paper, we calculate for the first time the angular distribution of \(\Lambda _b\rightarrow \Lambda \ell _1^-\ell _2^+\) decays using a full base of NP operators (partial results are available in [67, 68]). To achieve this, we use the decomposition of the \(\Lambda _b\rightarrow \Lambda \) hadronic matrix elements in [69] and the lattice QCD determination of the corresponding form factors [70]. We then use model independent constraints to derive upper bounds for the branching ratio of \(\Lambda _b\rightarrow \Lambda \ell _1^-\ell _2^+\) decays and specific models to provide predictions in a few scenarios [44, 45].

This paper is organised as follows: in Sect. 2 we highlight the main steps of our calculation and provide numerical results in a generic scenario. In Sect. 3 we use constraints on various LFV mesonic decays to put bounds on the branching ratio of \(\Lambda _b\rightarrow \Lambda \ell _1^-\ell _2^+\), for different choices of leptons in the final state and make predictions for specific models. We conclude in Sect. 4.

2 The angular distribution of \( \Lambda _b\rightarrow \Lambda \, \ell _1^{-} \ell _2^{+}\)

In this section, we introduce the concepts that we need for the study of phenomenological aspects in Sect. 3.

We consider the following effective Hamiltonian for LFV \(b\rightarrow s\ell _1^-\ell _2^+\) transitions:

$$\begin{aligned} {\mathcal {H}}_{\mathrm {eff}}= & {} -\frac{4 G_F}{\sqrt{2}}V_{tb}V_{ts}^* \frac{\alpha _{\mathrm {em}}}{4\pi } \sum _{i=9,10,S,P} \Bigg (C^{\ell _1\ell _2}_i(\mu ){\mathcal {O}}^{\ell _1\ell _2}_i(\mu )\nonumber \\&\quad +C_i^{\prime \ell _1\ell _2}(\mu ){\mathcal {O}}_i^{\prime \ell _1\ell _2}(\mu )\Bigg ), \end{aligned}$$
(1)

where the relevant operators are defined by

$$\begin{aligned} {\mathcal {O}}_{9}^{\ell _1\ell _2}&=({\bar{s}}\gamma _\mu P_L b)(\bar{\ell _1}\gamma ^\mu \ell _{2}),&{\mathcal {O}}_{10}^{\ell _1\ell _2}&= ({\bar{s}}\gamma _\mu P_L b)(\bar{\ell _1}\gamma ^\mu \gamma ^5\ell _{2}),\nonumber \\ {\mathcal {O}}_{S}^{\ell _1\ell _2}&= ({\bar{s}}P_R b)(\bar{\ell _1}\ell _{2}),&{\mathcal {O}}_{P}^{\ell _1\ell _2}&= ({\bar{s}}P_R b)(\bar{\ell _1}\gamma _5\ell _{2}), \nonumber \\ {\mathcal {O}}_{T}^{\ell _1\ell _2}&= ({\bar{s}} \sigma ^{\mu \nu } b) (\bar{\ell }_{1} \sigma _{\mu \nu } \ell _2),&{\mathcal {O}}_{T5}^{\ell _1\ell _2}&= ({\bar{s}} \sigma ^{\mu \nu } b) (\bar{\ell }_{1} \sigma _{\mu \nu }\gamma _5 \ell _2), \end{aligned}$$
(2)

and the operators with flipped chirality \({\mathcal {O}}^{\prime \ell _1\ell _2}_{i}\) are obtained from \({\mathcal {O}}^{\ell _1\ell _2}_{i}\) by replacing \(P_L \leftrightarrow P_R\), where \(P_{L/R}=\frac{1}{2}(1\mp \gamma _5)\). Notice that the operator \({\mathcal {O}}_7\):

$$\begin{aligned} {\mathcal {O}}_{7} = \frac{m_b}{e}({\bar{s}}\sigma _{\mu \nu }P_R b)F^{\mu \nu } \end{aligned}$$
(3)

cannot generate LFV contributions due to the universality of electromagnetic interactions. We parametrise the hadronic matrix elements for \(\Lambda _b(p,s_{{\Lambda _b}})\rightarrow \Lambda (k,s_{\Lambda })\) decays using an helicity decomposition [69,70,71,72]:

$$\begin{aligned}&\langle \Lambda (k,s_{\Lambda }) | {\bar{s}} \,\gamma ^\mu \, b | \Lambda _b(p,s_{{\Lambda _b}}) \rangle \nonumber \\&\quad =\,+ {\bar{u}}_\Lambda (k,s_{\Lambda }) \bigg [ f_0(q^2)\, (m_{\Lambda _b}-m_\Lambda )\frac{q^\mu }{q^2} \nonumber \\&\qquad + f_+(q^2) \frac{m_{\Lambda _b}+m_\Lambda }{s_+}\left( p^\mu + k^{\mu } - (m_{\Lambda _b}^2-m_\Lambda ^2)\frac{q^\mu }{q^2} \right) \nonumber \\&\qquad + f_\perp (q^2) \left( \gamma ^\mu - \frac{2m_\Lambda }{s_+} p^\mu - \frac{2 m_{\Lambda _b}}{s_+} k^{ \mu } \right) \bigg ] u_{\Lambda _b}(p,s_{{\Lambda _b}}), \end{aligned}$$
(4)
$$\begin{aligned}&\langle \Lambda (k,s_{\Lambda }) | {\bar{s}} \, \gamma ^\mu \gamma _5\, b | \Lambda _b(p,s_{{\Lambda _b}}) \rangle \, \nonumber \\&\quad =-{\bar{u}}_\Lambda (k,s_{\Lambda }) \,\gamma _5 \bigg [ g_0(q^2)\, (m_{\Lambda _b}+m_\Lambda )\frac{q^\mu }{q^2} \nonumber \\&\qquad + g_+(q^2)\frac{m_{\Lambda _b}-m_\Lambda }{s_-}\left( p^\mu + k^{\mu } - (m_{\Lambda _b}^2-m_\Lambda ^2)\frac{q^\mu }{q^2} \right) \nonumber \\&\qquad + g_\perp (q^2) \left( \gamma ^\mu + \frac{2m_\Lambda }{s_-} p^\mu - \frac{2 m_{\Lambda _b}}{s_-} k^{\mu } \right) \bigg ] u_{\Lambda _b}(p,s_{{\Lambda _b}}), \end{aligned}$$
(5)
$$\begin{aligned}&\big \langle {\Lambda (k,s_\Lambda )}\big \vert {\bar{s}} i \sigma ^{\mu \nu } b\big \vert {\Lambda _b(p,s_{\Lambda _b})}\big \rangle \nonumber \\&\quad = \, + {\bar{u}}_{\Lambda }(k, s_\Lambda ) \Big \{2h_+(q^2)\frac{p^\mu k^{ \nu }-p^\nu k^{\mu }}{s_{+}} \nonumber \\&\qquad +h_\perp (q^2)\Big [\frac{m_{\Lambda _{b}}+m_{\Lambda }}{q^2}(q^\mu \gamma ^\nu -q^\nu \gamma ^\mu ) \nonumber \\&\qquad -2\left( \frac{1}{q^2}+\frac{1}{s_{+}}\right) (p^\mu k^{\nu }-p^\nu k^{\mu }) \Big ] \nonumber \\&\qquad +\widetilde{h}_+ (q^2)\Big [i\sigma ^{\mu \nu }-\frac{2}{s_{-}}(m_{\Lambda _{b}}(k^{\mu }\gamma ^\nu -k^{\nu }\gamma ^\mu )\nonumber \\&\qquad -m_{\Lambda }(p^\mu \gamma ^\nu -p^\nu \gamma ^\mu )+p^\mu k^{\nu }-p^\nu k^{\mu }) \Big ] \nonumber \\&\qquad +\widetilde{h}_\perp (q^2) \frac{m_{\Lambda _{b}}-m_{\Lambda }}{q^2 s_{-}}\Big [(m_{\Lambda _{b}}^2-m_{\Lambda }^2-q^2)(\gamma ^\mu p^\nu - \gamma ^\nu p^\mu )\nonumber \\&\qquad -(m_{\Lambda _{b}}^2-m_{\Lambda }^2+q^2)(\gamma ^\mu k^{\nu }-\gamma ^\nu k^{\mu })\nonumber \\&\qquad +2(m_{\Lambda _{b}}-m_{\Lambda })(p^\mu k^{\nu }-p^\nu k^{\mu }) \Big ] \Big \} u_{\Lambda _b}(p, s_{\Lambda _b})\, \end{aligned}$$
(6)

with \(q=p-k\), \(s_\pm =(m_{\Lambda _b} \pm m_\Lambda )^2-q^2\), and \(s_{{\Lambda _b}}\) and \(s_{\Lambda }\) are the spin of the \(\Lambda _b\) and \(\Lambda \) baryons, respectively. Applying equations of motion to Eqs. (4)–(6), we obtain the following matrix elements for scalar and pseudoscalar operators:

$$\begin{aligned}&\langle \Lambda (k,s_{\Lambda })| {\bar{s}} \, b|\Lambda _b(p,s_{{\Lambda _b}})\rangle \nonumber \\&\quad = \frac{m_{{\Lambda _b}}-m_{\Lambda }}{m_b(\mu ) -m_s(\mu )}f_0{\bar{u}}_\Lambda (k,s_{\Lambda }) u_{{\Lambda _{b}}}(p,s_{{\Lambda _b}}), \end{aligned}$$
(7)
$$\begin{aligned}&\langle \Lambda (k,s_{\Lambda })| {\bar{s}} \gamma _5 b|\Lambda _b(p,s_{{\Lambda _b}})\rangle \nonumber \\&\quad =\frac{m_{{\Lambda _b}}+m_{\Lambda }}{m_b(\mu ) +m_s(\mu )}g_0{\bar{u}}_\Lambda (k,s_{\Lambda })\gamma _5 u_{{\Lambda _{b}}}(p,s_{{\Lambda _b}}), \end{aligned}$$
(8)

which agree with the expressions in Ref. [69]. In the following, we take the masses in \(\overline{\mathrm{MS}}\) using \({\overline{m}}_b({\overline{m}}_b)=4180\) MeV [73] and \({\overline{m}}_s({\overline{m}}_b)=78\,\text {MeV}\)  [74].

2.1 Differential decay width and numerical analysis

We decompose the spin-independent double-differential decay width as

$$\begin{aligned}&\frac{1}{\Gamma ^{(0)}}\frac{d\Gamma (\Lambda _b(p,s_{{\Lambda _b}})\rightarrow \Lambda (k,s_{\Lambda }) \ell _1^-(p_1) \ell _2^+(p_2))}{d\cos \theta d q^2}\nonumber \\&\quad = a + b \cos \theta _\ell + c \cos ^2\theta _\ell , \end{aligned}$$
(9)

with \(\Gamma ^{(0)}=\frac{ \alpha _{\mathrm {em}}^2 G_{\mathrm {F}}^2 |V_{tb}V_{ts}^*|^2}{2048\pi ^5m_{{\Lambda _b}}^3 q^2}\sqrt{\lambda _H}\sqrt{\lambda _L}\). We define \(\lambda _H\equiv \lambda (m_{{\Lambda _b}}^2, m_{\Lambda }^2, q^2)\) and \(\lambda _L\equiv \lambda (q^2, m_{\ell _1}^2, m_{\ell _2}^2)\), where \(\lambda \) is the usual Källén function defined as \(\lambda (a,b,c) = a^2+b^2+c^2-2 a (b+c)-2bc\). Here \(\cos \theta _\ell \) is the helicity angle in the dilepton frame as defined in Appendix A. The coefficients a, b and c are one of the main result of this work and have been calculated using the operator base in Eq. (1) and the decomposition for the hadronic matrix elements in Eqs. (4)–(6). We find:

$$\begin{aligned} a&=-\frac{1}{q^2} \bigg \{|f_0|^2\frac{(m_{{\Lambda _b}}-m_{\Lambda })^2}{q^2}s_+[|C^{\ell _1\ell _2}_{10+} |^2(m_{\ell _{1}}+m_{\ell _{2}})^2q_{- } \nonumber \\&\quad + |C^{\ell _1\ell _2}_{9+}|^2 (m_{\ell _{1}}- m_{\ell _{2}})^2q_{+ }]\nonumber \\&\quad +|f_\perp |^2 s_-[|C^{\ell _1\ell _2}_{10+}|^2 (\lambda _L+2q^2 q_{+ })\nonumber \\&\quad +|C^{\ell _1\ell _2}_{9+}|^2(\lambda _L+2q^2 q_{-})]\nonumber \\&\quad + |f_+|^2(m_{{\Lambda _b}}+m_{\Lambda })^2 s_-[|C^{\ell _1\ell _2}_{10+}|^2 \,q_{+} +|C^{\ell _1\ell _2}_{9+}|^2 \, q_{-}]\nonumber \\&\quad +|g_0|^2\frac{(m_{{\Lambda _b}}+m_{\Lambda })^2}{q^2}s_-[|C^{\ell _1\ell _2}_{10-}|^2(m_{\ell _{1}} + m_{\ell _{2}})^2q_{- }\nonumber \\&\quad +|C^{\ell _1\ell _2}_{9-}|^2(m_{\ell _{1}}- m_{\ell _{2}})^2q_{+ }] \nonumber \\&\quad +|g_\perp |^2 s_+[|C^{\ell _1\ell _2}_{10-}|^2(\lambda _L+2q^2 q_{+ })\nonumber \\&\quad +|C^{\ell _1\ell _2}_{9-}|^2(\lambda _L+2q^2 q_{- })]\nonumber \\&\quad +|g_+|^2(m_{{\Lambda _b}}-m_{\Lambda })^2 s_+[|C^{\ell _1\ell _2}_{10-}|^2 q_{+ } +|C^{\ell _1\ell _2}_{9-}|^2 q_{- }]\nonumber \\&\quad +16|h_+|^2s_-[(m_{\ell _1}+m_{\ell _2})^2 q_- |C^{\ell _1\ell _2}_{T}|^2\nonumber \\&\quad +(m_{\ell _1}-m_{\ell _2})^2 q_+|C^{\ell _1\ell _2}_{T5}|^2] \nonumber \\&\quad +16|{\tilde{h}}_+|^2s_+[(m_{\ell _1}-m_{\ell _2})^2 q_- |C^{\ell _1\ell _2}_{T}|^2\nonumber \\&\quad +(m_{\ell _1}+m_{\ell _2})^2 q_+ |C^{\ell _1\ell _2}_{T5}|^2]\nonumber \\&\quad +16|h_\perp |^2\frac{s_-}{q^2}(m_{{\Lambda _b}}+m_{\Lambda })^2\nonumber \\&\quad \times [q_-((m_{\ell _1}+m_{\ell _2})^2+q^2) |C^{\ell _1\ell _2}_{T}|^2\nonumber \\&\quad +q_+((m_{\ell _1}+m_{\ell _2})^2+q^2)|C^{\ell _1\ell _2}_{T5}|^2] \nonumber \\&\quad +16|{\tilde{h}}_\perp |^2\frac{s_+}{q^2}(m_{{\Lambda _b}}-m_{\Lambda })^2[q_+((m_{\ell _1} -m_{\ell _2})^2+q^2)\nonumber \\&\quad \times |C^{\ell _1\ell _2}_{T}|^2+q_-((m_{\ell _1}+m_{\ell _2})^2+q^2) |C^{\ell _1\ell _2}_{T5}|^2]\bigg \}\nonumber \\&\quad -(q_- |C_{P-}^{\ell _1\ell _2}|^2+q_+|C_{S-}^{\ell _1\ell _2}|^2)\frac{s_- (m_{{\Lambda _b}}+m_{\Lambda })^2}{(m_b+m_s)^2}|g_0|^2\nonumber \\&\quad -(q_- |C_{P+}^{\ell _1\ell _2}|^2+q_+|C_{S+}^{\ell _1\ell _2}|^2)\frac{s_+ (m_{{\Lambda _b}}-m_{\Lambda })^2}{(m_b-m_s)^2}|f_0|^2\nonumber \\&\quad -\frac{2s_+}{q^2} \frac{(m_{{\Lambda _b}}-m_{\Lambda })^2}{m_b-m_s}[{\mathrm {Re}}(C^{\ell _1\ell _2}_{10+} C^{*\ell _1\ell _2}_{P+})(m_{\ell _{1}}+m_{\ell _{2}})q_-\nonumber \\&\quad +{\mathrm {Re}}(C^{\ell _1\ell _2}_{9+} C^{*\ell _1\ell _2}_{S+})(m_{\ell _{1}}-m_{\ell _{2}})q_+]|f_0|^2 \nonumber \\&\quad -\frac{2s_-}{q^2} \frac{(m_{{\Lambda _b}}+m_{\Lambda })^2}{m_b+m_s}[{\mathrm {Re}}(C^{\ell _1\ell _2}_{10-} C^{*\ell _1\ell _2}_{P-})(m_{\ell _{1}}+m_{\ell _{2}})q_-\nonumber \\&\quad +{\mathrm {Re}}(C^{\ell _1\ell _2}_{9-} C^{*\ell _1\ell _2}_{S-})(m_{\ell _{1}}-m_{\ell _{2}})q_+]|g_0|^2 \nonumber \\&\quad - \frac{8}{q^2} \bigg \{ (m_{\ell _1} +m_{\ell _2})(m_{{\Lambda _b}}+m_{\Lambda }) s_- q_- \text {Re}(C_{9+}^{\ell _{1} \ell _2} \,C_T^{* \ell _{1} \ell _2})\nonumber \\&\quad \times \big [ \text {Re}(f_+ \,h_+^{*}) +2\text {Re}(f_\perp \,h_\perp ^{*}) \big ] \nonumber \\&\quad + (m_{\ell _1} - m_{\ell _2}) (m_{{\Lambda _b}}-m_{\Lambda }) s_+ q_+ \text {Re}(C_{10-}^{\ell _{1} \ell _2} \,C_T^{* \ell _{1} \ell _2}) \nonumber \\&\quad \times \big [ \text {Re}(g_+ \, \tilde{h}_+^{*}) +2\text {Re}(g_\perp \,\tilde{h}_\perp ^{*}) \big ] \nonumber \\&\quad + (m_{\ell _1} +m_{\ell _2}) (m_{{\Lambda _b}}- m_{\Lambda }) s_+ q_- \text {Re}(C_{9-}^{\ell _{1} \ell _2} \,C_{T5}^{* \ell _{1} \ell _2})\nonumber \\&\quad \times \big [ \text {Re}(g_+ \,\tilde{h}_+^{*} +2\text {Re}(g_\perp \, \tilde{h}_\perp ^{*})) \big ] \nonumber \\&\quad +(m_{\ell _1} - m_{\ell _2}) (m_{{\Lambda _b}}+ m_{\Lambda }) s_- q_+ \text {Re}(C_{10+}^{\ell _{1} \ell _2} \,C_{T5}^{* \ell _{1} \ell _2})\nonumber \\&\quad \times \big [ \text {Re}(f_+ \,\tilde{h}_+^{*} +2\text {Re}(f_\perp \, \tilde{h}_\perp ^{*})) \big ] \bigg \}, \end{aligned}$$
(10)
$$\begin{aligned} b&= \frac{2}{q^4} \big [{\mathrm {Re}} (f_0 f_+^*) (|C^{\ell _1\ell _2}_{9+}|^2+|C^{\ell _1\ell _2}_{10+}|^2) \nonumber \\&\quad +{\mathrm {Re}} (g_0 g_+^*) (|C^{\ell _1\ell _2}_{9-}|^2+|C^{\ell _1\ell _2}_{10-}|^2) \big ]\sqrt{\lambda _H\lambda _L}\nonumber \\&\quad \times (m_{\ell _{2}}^2-m_{\ell _{1}}^2)(m_{{\Lambda _b}}^2-m_{\Lambda }^2)-4[ {\mathrm {Re}}(C^{\ell _1\ell _2}_{9-} C^{*\ell _1\ell _2}_{10+})\nonumber \\&\quad +{\mathrm {Re}}(C^{\ell _1\ell _2}_{9+}C^{*\ell _1\ell _2}_{10-})] {\mathrm {Re}}(f_\perp g_\perp ^*) \sqrt{\lambda _H\lambda _L}\nonumber \\&\quad -\frac{2(m_{{\Lambda _b}}^2-m_{\Lambda }^2)}{q^2(m_b-m_s)}\sqrt{\lambda _H\lambda _L} [{\mathrm {Re}}(C^{\ell _1\ell _2}_{10+} C^{*\ell _1\ell _2}_{P+})(m_{\ell _{1}}-m_{\ell _{2}})\nonumber \\&\quad +{\mathrm {Re}}(C^{\ell _1\ell _2}_{9+} C^{*\ell _1\ell _2}_{S+})(m_{\ell _{1}}+m_{\ell _{2}})]{\mathrm {Re}}(f_0 f_+^*) \nonumber \\&\quad -\frac{2(m_{{\Lambda _b}}^2-m_{\Lambda }^2)}{q^2(m_b+m_s)}\sqrt{\lambda _H\lambda _L}\nonumber \\&\quad \times [{\mathrm {Re}}(C^{\ell _1\ell _2}_{10-} C^{*\ell _1\ell _2}_{P-})(m_{\ell _{1}}-m_{\ell _{2}})+{\mathrm {Re}}(C^{\ell _1\ell _2}_{9-} C^{*\ell _1\ell _2}_{S-})\nonumber \\&\quad \times (m_{\ell _{1}}+m_{\ell _{2}})]{\mathrm {Re}}(g_0 g_+^*) \nonumber \\&\quad +\frac{64}{q^4}(|C^{\ell _1\ell _2}_{T}|^2+|C^{\ell _1\ell _2}_{T5}|^2)\nonumber \\&\quad \times (m_{\ell _2}^2-m_{\ell _1}^2)(m_{{\Lambda _b}}^2-m_{\Lambda }^2)\sqrt{\lambda _H \lambda _L}\text {Re}(h_\perp {\tilde{h}}_\perp ^*) \nonumber \\&\quad - \frac{8}{q^2} \bigg \{ \text {Re}(C_{9+}^{\ell _1 \ell _2}\, C_T^{*\ell _1 \ell _2}) (m_{\ell _1} -m_{\ell _2}) (m_{{\Lambda _b}}- m_{\Lambda }) \sqrt{\lambda _H \lambda _L}\nonumber \\&\quad \times \big [ \text {Re}(f_0 h_+^{*}) + 2\text {Re}( f_\perp \tilde{h}_\perp ^{*}) \big ] + \text {Re}(C_{10-}^{\ell _1 \ell _2}\, C_T^{*\ell _1 \ell _2})\nonumber \\&\quad \times (m_{\ell _1} +m_{\ell _2}) (m_{{\Lambda _b}}+ m_{\Lambda }) \sqrt{\lambda _H \lambda _L} \big [ \text {Re}(g_0 \tilde{h}_+^{*}) \nonumber \\&\quad + 2\text {Re}( g_\perp h_\perp ^{*}) \big ] + q^2 \frac{(m_{{\Lambda _b}}- m_{\Lambda })}{(m_b- m_s)} \text {Re}(C_{S+}^{\ell _1 \ell _2}\, C_T^{*\ell _1 \ell _2}) \nonumber \\&\quad \times \sqrt{\lambda _H \lambda _L} \text {Re}(f_0 \,h_+^{*}) +q^2 \frac{(m_{{\Lambda _b}}+ m_{\Lambda })}{(m_b + m_s)} \text {Re}(C_{P-}^{\ell _1 \ell _2}\, C_T^{*\ell _1 \ell _2})\nonumber \\&\times \sqrt{\lambda _H \lambda _L} \text {Re}(g_0 \,\tilde{h}_+^{*} ) + \text {Re}(C_{9-}^{\ell _1 \ell _2}\, C_{T5}^{*\ell _1 \ell _2}) \nonumber \\&\quad \times (m_{\ell _1} -m_{\ell _2}) (m_{{\Lambda _b}}+ m_{\Lambda }) \sqrt{\lambda _H \lambda _L} \big [ \text {Re}(g_0 \tilde{h}_+^{*}) \nonumber \\&\quad + 2\text {Re}( g_\perp h_\perp ^{*}) \big ]+ \text {Re}(C_{10+}^{\ell _1 \ell _2}\, C_{T5}^{*\ell _1 \ell _2})\nonumber \\&\quad \times (m_{\ell _1} +m_{\ell _2}) (m_{{\Lambda _b}}- m_{\Lambda }) \sqrt{\lambda _H \lambda _L} \nonumber \\&\quad \times \big [ \text {Re}(f_0 h_+^{*}) + 2\text {Re}( f_\perp \tilde{h}_\perp ^{*}) \big ] \nonumber \\&\quad +q^2 \frac{(m_{{\Lambda _b}}+ m_{\Lambda })}{(m_b + m_s)} \text {Re}(C_{S-}^{\ell _1 \ell _2} C_{T5}^{*\ell _1 \ell _2}) \sqrt{\lambda _H \lambda _L} \text {Re}(g_0 \tilde{h}_+^{*})\nonumber \\&\quad + q^2 \frac{(m_{{\Lambda _b}}- m_{\Lambda })}{(m_b -m_s)}\text {Re}(C_{P+}^{\ell _1 \ell _2} C_{T5}^{*\ell _1 \ell _2}) \sqrt{\lambda _H \lambda _L} \text {Re}(f_0 h_+^{*}) \bigg \} ,\end{aligned}$$
(11)
$$\begin{aligned} c&=+ (|C^{\ell _1\ell _2}_{9+}|^2+|C^{\ell _1\ell _2}_{10+}|^2)\nonumber \\&\quad \times \frac{\lambda _L\lambda _H}{q^2 s_+}\bigg [-|f_+|^2 \frac{(m_{{\Lambda _b}}+m_{\Lambda })^2}{q^2} +|f_\perp |^2\bigg ] \nonumber \\&\quad + (|C^{\ell _1\ell _2}_{9-}|^2+|C^{\ell _1\ell _2}_{10-}|^2)\nonumber \\&\quad \times \frac{\lambda _L\lambda _H}{q^2 s_-}\bigg [-|g_+|^2 \frac{(m_{{\Lambda _b}}-m_{\Lambda })^2}{q^2 }+|g_\perp |^2\bigg ]\nonumber \\&\quad +(|C^{\ell _1\ell _2}_{T}|^2+|C^{\ell _1\ell _2}_{T5}|^2) \frac{4\lambda _L}{q^2}\bigg [s_- |h_+|^2+s_+|{\tilde{h}}_+|^2 \nonumber \\&\quad -\frac{(m_{{\Lambda _b}}+m_{\Lambda })^2s_-}{q^2}|h_\perp |^2 -\frac{(m_{{\Lambda _b}}-m_{\Lambda })^2s_+}{q^2}|{\tilde{h}}_\perp |^2\bigg ], \end{aligned}$$
(12)

with \(C^{\ell _1 \ell _2}_{X \pm } = (C^{\ell _1 \ell _2}_{X} \pm C^{\prime \ell _1 \ell _2}_{X})\), \(q_{\pm }= (m_{\ell _{1}}\pm m_{\ell _{2}})^2-q^2\) and \(\sigma ^{\mu \nu } = i/2 [\gamma ^\mu ,\gamma ^\nu ]\). Our results agree with [72] when setting \(m_{\ell _{1}}=m_{\ell _{2}}=0\) and with [75] for \(m_{\ell _{1}}=m_{\ell _{2}}\). However, we note that our convention for the helicity angle \(\cos \theta _\ell \) has the opposite sign that the one in [72] . We also note that our formulae above disagree with Ref. [68], in the specific in terms proportional to \((m_{\ell _1}-m_{\ell _2})^2|C^{\ell _1\ell _2}_{9+}|^2 |f_0|^2\) and \((m_{\ell _1}-m_{\ell _2})^2|C^{\ell _1\ell _2}_{9+}|^2 |g_0|^2\). We ascribe these differences to an incorrect treatment of mass effects in [68]. Our formulae for the angular coefficients a, b, c also hold for \(\Lambda _b\rightarrow \Lambda ^*\ell _1^-\ell _2^+\) decays, when setting the additional perpendicular \(\Lambda _b\rightarrow \Lambda ^*\) form factor to zero. This is a reasonable approximation as in the Heavy-Quark-Expansion, this form factor is suppressed by \(\Lambda _{\mathrm {QCD}}/m_b\) [76,77,78].

In the following, we focus on the branching ratio and forward–backward asymmetry:

$$\begin{aligned} \frac{d{\mathcal {B}}^{\ell _1\ell _2}}{d q^2}&= 2\Gamma ^{(0)} \tau _{\Lambda _b}\left( a + \frac{c}{3}\right) , \end{aligned}$$
(13)
$$\begin{aligned} \frac{d A_{\mathrm {FB}}^{\ell _1\ell _2}}{dq^2}&= \frac{\int _0^1 d\cos \theta \frac{d\Gamma }{d\cos \theta d q^2}-\int _{-1}^0 d\cos \theta \frac{d\Gamma }{d\cos \theta d q^2 }}{\int _0^1 d\cos \theta \frac{d\Gamma }{d\cos \theta d q^2}+\int _{-1}^0 d\cos \theta \frac{d\Gamma }{d\cos \theta d q^2 }}\nonumber \\&= \frac{b}{2\left( a+\frac{c}{3}\right) }, \end{aligned}$$
(14)

where the branching ratio \({\mathcal {B}} = \tau _{\Lambda _b} \Gamma \), where \(\tau _{\Lambda _b}\) is the mean life of the \(\Lambda _b\) baryon [79] and \(\Gamma \) the total width. To evaluate the size of LFV \(\Lambda _b\rightarrow \Lambda \ell _1^-\ell _2^+\) decay, we provide the \(q^2\)-integrated quantities of Eqs. (13)–(14). For simplicity, we set \(C_i^{\prime \ell _1\ell _2}=0\). We further set \(C^{\ell _1\ell _2}_{T(5)}=0\). This choice is discussed at the beginning of Sect. 3. Using the values for the masses from PDG [73], CKM factors from the UT-fit collaboration [80] and lattice QCD inputs for the form factors [70], we obtain

$$\begin{aligned} 10^{8} \cdot {\mathcal {B}}^{\ell _1\ell _2}&=\, \xi _9^{\ell _1\ell _2} |C_{9}^{\ell _1\ell _2}|^2 +\xi _{10}^{\ell _1\ell _2} |C_{10}^{\ell _1\ell _2}|^2 +\xi _S^{\ell _1\ell _2} |C_{S}^{\ell _1\ell _2}|^2\nonumber \\&\quad + \xi _P^{\ell _1\ell _2}|C_{P}^{\ell _1\ell _2}|^2 +\xi _{9S}^{\ell _1\ell _2}{\mathrm {Re}}( C_{9}^{\ell _1\ell _2} C_{S}^{*\ell _1\ell _2} ) \nonumber \\&\quad +\xi _{10P}^{\ell _1\ell _2} \, {\mathrm {Re}}( C_{10}^{\ell _1\ell _2} C_{P}^{*\ell _1\ell _2} ), \end{aligned}$$
(15)
$$\begin{aligned}&A_{\mathrm {FB}}^{\ell _1\ell _2} = \bigg [\rho ^{\ell _1\ell _2} (|C_{10}^{\ell _1\ell _2}|^2 +|C_{9}^{\ell _1\ell _2}|^2 )+\rho _{910}^{\ell _1\ell _2} {\mathrm {Re}}( C_{9}^{\ell _1\ell _2} C_{10}^{*\ell _1\ell _2}) \nonumber \\&\quad +\rho _{9S}^{\ell _1\ell _2}{\mathrm {Re}}( C_{9}^{\ell _1\ell _2} C_{S}^{*\ell _1\ell _2} ) +\rho _{10P}^{\ell _1\ell _2}{\mathrm {Re}}( C_{10}^{\ell _1\ell _2} C_{P}^{*\ell _1\ell _2} )\bigg ]\bigg / \Gamma ^{\ell _1 \ell _2}, \end{aligned}$$
(16)

with the numerical values for the coefficients \(\xi ^{\ell _1\ell _2}_i\) and \(\rho ^{\ell _1\ell _2}_i\) listed in Tables 1 and 2 and \(\Gamma ^{\ell _1\ell _2}\) is the integrated width. We present only explicit results for the final states \(\tau ^{\pm }\mu ^{\mp }\) and \(\mu ^{\pm }e^{\mp }\). The results for \(\tau ^\pm e^{\mp }\) can easily be obtained from the above results and agree within 1\(\sigma \) with those of \(\tau ^\pm \mu ^\mp \). The quoted uncertainties only include those from the form factors, which are the dominant ones. The correlation matrices between the two sets of \(\lbrace \xi ^{\ell _1\ell _2}_i,\rho ^{\ell _1\ell _2}_i\rbrace \) coefficients are given in Appendix B. The \(\xi _i^{\ell _1\ell _2}\) coefficients in Table 1 do not depend on the charges of the final state leptons, except for \( \xi ^{\mu \tau }_{9S}\) which depends on \((m_{\ell _{1}}-m_{\ell _{2}})\) and thus switches sign when switching the charges of the final state leptons, i.e. \(\xi ^{\tau \mu }_{9S} = - \xi ^{\mu \tau }_{9S}\). Besides, we note that for \(\mu e\) final states, \(\xi _9^{\ell _1\ell _2} = \xi ^{\ell _1\ell _2}_{10}\) and \(\xi _S^{\ell _1\ell _2}=\xi _P^{\ell _1\ell _2}\), such that only the combination \(|C_9^{\ell _1\ell _2}|^2 + |C_{10}^{\ell _1\ell _2}|^2\) and \(|C_P^{\ell _1\ell _2}|^2 + |C_{S}^{\ell _1\ell _2}|^2\) can be constrained. The coefficients \(\rho _i^{\ell _1\ell _2}\) in Table 2 are reported in units of \(10^{-21}\) \(\hbox {GeV}^{-1}\), which is then compensated in \(A_{\mathrm {FB}}^{\ell _1\ell _2}\) by the size of the decay width.

Table 1 Numerical values for the parameters of Eq. (15). The coefficients do not depend on the charges of the final state leptons, except for \(\xi ^{\ell _1\ell _2}_{9S}\) which changes sign when switching the charges of the leptons, i.e. \(\xi ^{\tau \mu }_{9S} = - \xi ^{\mu \tau }_{9S}\) and \(\xi ^{e \mu }_{9S} = - \xi ^{\mu e}_{9S}\). The uncertainties only include those from the form factor which are the dominant ones
Table 2 Coefficients for the numerator of \(A_{\mathrm {FB}}^{\ell _1 \ell _2}\). We give the values in units of \(10^{-21}~\text {GeV}^{-1}\). This factor is compensated by the size of the decay width in Eq. (16)

3 Phenomenological implications

In the following, we discuss the implications of the available constraints on LFV B-meson decays and which bounds they imply on the observables in the baryonic modes. In order to do so, we need to choose which NP operators are present. Since no NP particles have been observed so far above the electroweak scale, we choose to work with the SMEFT:

$$\begin{aligned} {\mathcal {L}}_\text {eff}= & {} {\mathcal {L}}_\text {SM} - \frac{1}{M^2}\bigg \{[{\mathcal {C}}_{lq}^{(3)}]^{ij\alpha \beta }({\bar{Q}}^i \gamma ^\mu \sigma ^a Q^j)({\bar{L}}^\alpha \gamma _\mu \sigma ^a L^\beta )\nonumber \\&+[{\mathcal {C}}_{lq}^{(1)}]^{ij\alpha \beta }({\bar{Q}}^i \gamma ^\mu Q^j)({\bar{L}}^\alpha \gamma _\mu L^\beta )\nonumber \\&+[{\mathcal {C}}_{ledq}]^{ij\alpha \beta }({\bar{Q}}^i d_R^j)({\bar{e}}_R^\alpha L^\beta ) \bigg \} , \end{aligned}$$
(17)

where we adopt the so-called Warsaw basis [81]. Here we denote with Q and L the left-handed quark and lepton doublets, respectively, and with \(e_R\) and \(d_R\) the right-handed charged leptons and down-type quarks, respectively. We further denote \(\epsilon = i\sigma _2\) and M is the effective scale which can be associated with the mass of the heavy NP degrees of freedom.

The operators in Eq. (17) are the complete set of dimension-6 semileptonic operators that can contribute to \(b\rightarrow s\ell _1\ell _2\) transitions. We note that none of these operators contain a tensor current; nonetheless, at low energy, the operator \({\mathcal {O}}_{T(5)}^{\ell _1\ell _2}\) defined in Eq. (2) could be generated through effective operators containing a covariant derivative [82]. However, as tensor operators provide a poor explanation for B anomalies (see e.g. [83]), we do not consider them in our analysis.

The Wilson coefficients of the operators in Eq. (17) can be constrained from low-energy processes as well as high-\(p_T\) data, and in general a flavour structure has to be assumed to reduce the number of independent NP parameters. In the following, we choose to consider only constraints from low-energy data and first do not to assume any hierarchy for the NP couplings. In Sect. 3.2 we then study particular scenarios, where a more complex structure for NP couplings in flavour space is assumed. For the \(b\rightarrow s \ell _1 \ell _2\) transition we are interested in, we set \(i=2\) and \(j=3\) and generic \(\alpha =\ell _1\) and \(\beta =\ell _2\) in Eq. (17). Performing now the tree-level matching onto Eq. (1), we have

$$\begin{aligned} \begin{aligned} C_{9}^{\ell _1 \ell _2} = - \,{C}_{10}^{\ell _1 \ell _2} =&+\frac{v^2}{\Lambda ^2}\frac{\pi }{\alpha _\text {em} |V_{tb}V^*_{ts}|}\\&\quad \times \left( [{\mathcal {C}}_{lq}^{(3)}]^{23\ell _1 \ell _2}+[{\mathcal {C}}_{lq}^{(1)}]^{23\ell _1 \ell _2}\right) , \\ C_{9}^{\prime \,\ell _1 \ell _2} = +\, {C}_{10}^{\prime \, \ell _1\ell _2} =&+ \frac{v^2}{\Lambda ^2}\frac{\pi }{\alpha _\text {em} |V_{tb}V^*_{ts}|} [{\mathcal {C}}_{l d}]^{23\ell _1 \ell _2}, \\ {C}_S^{\ell _1 \ell _2} = -\,{C}_{P}^{\ell _1 \ell _2} =&+\frac{v^2}{\Lambda ^2}\frac{\pi }{\alpha _\text {em} |V_{tb}V^*_{ts}|} \,[{\mathcal {C}}_{leqd}]^{23\ell _1 \ell _2}, \\ {C}_S^{\prime \,\ell _1 \ell _2} = +\,{C}_{P}^{\prime \, \ell _1\ell _2} =&+ \frac{v^2}{\Lambda ^2}\frac{\pi }{\alpha _\text {em} |V_{tb}V^*_{ts}|} \,[{\mathcal {C}}_{leqd}^*]^{32\ell _1 \ell _2}. \end{aligned} \end{aligned}$$
(18)

3.1 Model-independent approach

Table 3 Experimental upper limits for LFV B decays at \(90\%\) CL
Table 4 Predictions for the coefficients describing \(B^+\rightarrow K^+ \ell _1^-\ell _2^+\) decays using the hardonic form factors from Refs. [89, 90]. We note that these coefficients are independent of the charges of the leptons, except for \(c_{\ell _1\ell _2}^{S9}\) which changes sign depending on the charge of the heavier lepton
Fig. 1
figure 1

Model independent constrains on different combinations of Wilson coefficients obtained from the \(90\%\) CL upper limits on meson \(b\rightarrow s \mu ^\pm \tau ^\mp \) transitions

constraints on several combinations of Wilson coefficients from measurements of mesonic LFV decays. We consider the branching ratios of the decay modes \({\bar{B}}_s \rightarrow \ell _1^-\ell _2^+\) and \(B\rightarrow K \ell _1^-\ell _2^+\), for which the experimental upper limits at \(90\%\) CL are reported in Table 4. Using Eq. (1), we have:

$$\begin{aligned}&{\mathcal {B}}({\bar{B}}_s\rightarrow \ell ^-_1 \ell ^+_2)=\frac{\tau _{B_s}}{64\pi ^3}\frac{\alpha _\text {em}^2G_F^2 |V_{tb}V_{ts}^*|^2}{m_{B_s}^3} f_{B_s}^2 \, \lambda ^{1/2}(m_{B_s}^2,m_{\ell _1}^2,m_{\ell _2}^2)\nonumber \\&\quad \times \left\{ [m_{B_s}^2-(m_{\ell _1}-m_{\ell _2})^2]\left| (m_{\ell _1}+m_{\ell _2}){C}^{\ell _1\ell _2}_{10-}+\frac{m_{B_s}^2}{m_b+m_s}{C}^{\ell _1\ell _2}_{P-}\right| ^2\right. \nonumber \\&\quad +\left. [m_{B_s}^2-(m_{\ell _1}+m_{\ell _2})^2]\left| (m_{\ell _1}-m_{\ell _2})({C}^{\ell _1\ell _2}_{9-}) +\frac{m_{B_s}^2}{m_b+m_s}({C}^{\ell _1\ell _2}_{S-})\right| ^2\right\} , \nonumber \\ \end{aligned}$$
(19)

and

$$\begin{aligned}&{\mathcal {B}}(B^+\rightarrow K^+\ell _1^-\ell _2^+) \nonumber \\&\quad = 10^{-8} \bigg \{c^{9+}_{\ell _1\ell _2}\left| C_{9+}^{\ell _1\ell _2} \right| ^2 + c^{10+}_{\ell _1\ell _2}\left| C_{10+}^{\ell _1\ell _2} \right| ^2 + c^S_{\ell _1\ell _2}\left| C_{S+}^{\ell _1\ell _2} \right| ^2\nonumber \\&\qquad +c^P_{\ell _1\ell _2}\left| C_{P+}^{\ell _1\ell _2}\right| ^2 +c^{S9}_{\ell _1\ell _2}\,{\mathrm {Re}}[C_{S+}^{*\ell _1\ell _2} C_{9+}^{\ell _1\ell _2}] \nonumber \\&\qquad + c^{P10}_{\ell _1\ell _2}\,{\mathrm {Re}}[C_{P+}^{*\ell _1\ell _2} C_{10+}^{\ell _1\ell _2}] \bigg \}, \end{aligned}$$
(20)

Both Eqs. (19) and (20) agree with previous results in the literature [91, 92]. Using again the values for the masses from PDG [73], CKM factors from the UT-fit collaboration [80], \(f_{B_s}=215\,\text {MeV}\) [93] and Lattice QCD/Light Cone Sum Rule results in Refs. [89, 90], we find the coefficients \(c_{\ell _1\ell _2}^i\) in Eq. (20) as listed in Table 4. Similar as for \(\xi _{9S}\), the coefficient \(c_{\ell _1\ell _2}^{S9}\) is proportional to \(m_{\ell _1}-m_{\ell _2}\) and thus changes sign depending on charge of the heavier lepton. We stress that the numbers in Table 4 are strongly dependent on the choice for \(\alpha _{\mathrm {em}}\). Here we take \(\alpha _{\mathrm {em}}=1/133\). A different choice can be implemented by rescaling the \(c_{\ell _1\ell _2}^i\) coefficients.

Finally, we use the experimental upper bounds listed in Table 3 and Eqs. (19) and (20) to constrain different combinations of couplings \(C_{i}^{\ell _1, \ell _2}\). As stated before, we do not consider \(\tau e\) decays as the constraints coming from these decays are similar to those from the \(\tau \mu \) channel. Furthermore, for simplicity we also do not consider the \({\mathcal {O}}^{\prime \ell _1\ell _2}_i\) operators. This choice is motivated by the fact that these operators are unappealing when trying to fit \(b\rightarrow s \ell \ell \) data [62,63,64,65, 94]. Nevertheless, we stress that the baryonic channels have a different dependence on the primed operators with respect to the mesonic ones, which may be interesting to consider once scenarios involving these operators become more interesting to explain the B anomalies.

Fig. 2
figure 2

Model independent constrains on different combinations of Wilson coefficients obtained from the \(90\%\) CL upper limits on meson \(b\rightarrow s \mu ^\pm e^\mp \) transitions

The obtained bounds for \(\tau \mu \) and \(\mu e\) finals states are given in Fig. 1 and in Fig. 2, respectively. We consider three 2-dimensional scenarios, in which we allow only some combinations of NP Wilson coefficients to be non-zero: \(C_9^{\ell _1\ell _2}\) and \(C_{10}^{\ell _1\ell _2}\), \(C_S^{\ell _1\ell _2}\) and \(C_P^{\ell _1\ell _2}\) and the SMEFT inspired one, where \(C_9^{\ell _1\ell _2}=-C_{10}^{\ell _1\ell _2}\) and \(C_S^{\ell _1\ell _2}=-C_P^{\ell _1\ell _2}\). For the \(C_9^{\ell _1\ell _2}-C_{10}^{\ell _1\ell _2}\) and \(C_S^{\ell _1\ell _2}-C_P^{\ell _1\ell _2}\) scenarios, which are independent of the charge configuration in the final state, we only consider the strongest bound in Table 3. As the interference between \(C_9^{\ell _1\ell _2}\) and \(C_S^{\ell _1\ell _2}\) depends on the charge configuration of the leptons in the final state, we present plots for both the \(\tau ^+\mu ^-\) and \(\tau ^-\mu ^+\) final states. We note that the \({\bar{B}}_s\rightarrow \tau ^- \tau ^+\) decay only gives a very weak constraint in the \(C_9^{\ell _1\ell _2}=-C_{10}^{\ell _1\ell _2}\) plane ranging from \(-200\) to 200. From comparison of the plots in Fig. 1, we find large differences between the \(\tau ^+\mu ^-\) and the \(\tau ^-\mu ^+\). Hence, we stress that it is important to analyse these final states separately. For the electron, the differences between \(\mu ^-e^+\) and \(\mu ^+e^-\) are negligible and we only present one figure.

As the mesonic \(B^+\rightarrow K^+\ell _1\ell _2\) and \({{\bar{B}}}_s\rightarrow \ell _1\ell _2\) are mediated by the same quark level transition, we can use the obtained upper limits on combinations of Wilson coefficients and convert those into upper limit on the branching ratio and forward–backward asymmetry for \(\Lambda _b\rightarrow \Lambda \ell _1\ell _2\) decays using Eq. (16). When allowing for only one NP Wilson coefficient to be nonzero at a time, for example allowing only \(C_9^{\ell _1\ell _2}\ne 0\), the corresponding bounds can be easily obtained by calculating the scale factor between \(c^i_{\ell _1\ell _2}\) of the meson \(B\rightarrow K\) LFV decay and \(\xi _i^{\ell _1\ell _2}\) of the baryon \(\Lambda _b\rightarrow \Lambda \) decay using Tables 1 and 4 and re-scaling the upper limit of the mesonic decay accordingly. In addition, comparing the coefficients in these tables, we observe that the ratios \(c^{i}_{\ell _1\ell _2}/c^{j}_{\ell _1\ell _2}\) and \(\xi _{i}^{\ell _1\ell _2}/\xi _{j}^{\ell _1\ell _2}\) are very similar for \(i,j=9,10\) and \(i,j=S,P\). Therefore, the sensitivities for LFV \(B\rightarrow K\) and \(\Lambda _b\rightarrow \Lambda \) decays are rather similar when considering the \(C_9^{\ell _1\ell _2}-C_{10}^{\ell _1\ell _2}\) only and \(C_S^{\ell _1\ell _2}-C_P^{\ell _1\ell _2}\) only scenarios. Upper limits (at \(90\%\) CL) for the branching ratio of \(\Lambda _b\rightarrow \Lambda \ell _1\ell _2\) derived from their mesonic counter parts for the three scenarios are presented in Table 5. These values should be interpreted as follows: any future experimental upper limit on the baryonic mode below the quoted value gives stronger constraints on the Wilson coefficients than those obtained from the current mesonic upper limits.

Table 5 Upper limits for the branching ratio of \(\Lambda _b \rightarrow \Lambda \) LFV decays obtained in a model independent way by considering their mesonic counter parts. Bounds are at \(90\%\) CL. For the first two scenarios, the branching ratios are independent of the charge configuration. However for the SMEFT scenario this is not the case anymore, hence we present both branching ratios for \(\mu ^-\tau ^+\) and in brackets \(\tau ^-\mu ^+\)
Fig. 3
figure 3

Illustration of the orthogonality between current mesonic and possible future baryonic constraints

The complementarity of the mesonic and the baryonic LFV channels specifically arises when considering both (axial)vector and (pseudo)scalar operators. This complementarity is caused the difference between the ratios \(c^{i}_{\ell _1\ell _2}/c^{j}_{\ell _1\ell _2}\) and \(\xi _{i}^{\ell _1\ell _2}/\xi _{j}^{\ell _1\ell _2}\) for \(i=S9, P10\) and \(j=S,P,9,10\). We expect a similar complementarity also when both tensor operators and (axial)vector operators are present. We illustrate quantitatively this in Fig. 3 for the SMEFT scenario where \(C_9^{\ell _1\ell _2}=-C_{10}^{\ell _1\ell _2}, C_S^{\ell _1\ell _2}=-C_P^{\ell _1\ell _2}\). We present the current meson constraints combined with two possible constraints on the \(\Lambda _b \rightarrow \Lambda \mu ^-\tau ^+\) branching ratio and observe that the mesonic modes place strong constraints on scalar/pseudoscalar interactions while the baryonic channel is more sensitive to \(C_9^{\mu \tau }\) and \(C_{10}^{\mu \tau }\).

Finally, we consider the integrated forward–backward asymmetry \( A_{\mathrm {FB}}^{\ell _1\ell _2}\) which provides orthogonal information compared to the branching ratio. From Eq. (16) we note the following properties: \(A_{\mathrm {FB}}^{\ell _1\ell _2}\) is identically zero if \(C_9^{\ell _1\ell _2} = C_{10}^{\ell _1\ell _2}=0\), and in the case in which only \(C_9^{\ell _1\ell _2} = -C_{10}^{\ell _1\ell _2}\ne 0\) \(A_{\mathrm {FB}}^{\ell _1\ell _2}\) is independent on the values of \(C_9^{\ell _1\ell _2}\) and \(C_{10}^{\ell _1\ell _2}\). In the latter scenario, we find for \( C_9^{\ell _1\ell _2}= - C_{10}^{\ell _1\ell _2}\) and \( C_S^{\ell _1\ell _2} = C_P^{\ell _1\ell _2} =0\)

$$\begin{aligned}&A_{\mathrm {FB}}^{\tau \mu } = 0.14\pm 0.01,\quad A_{\mathrm {FB}}^{\mu \tau } = 0.40\pm 0.03,\nonumber \\&A_{\mathrm {FB}}^{e\mu } = A_{\mathrm {FB}}^{\mu e} = 0.33\pm 0.04. \end{aligned}$$
(21)

A measurement or an upper limit different from these values provides interesting complementary information. This is illustrated in Fig. 4, where we consider for the \(\mu ^- \tau ^+\) final state a future scenario in which an upper limit of \(A_{\mathrm{FB}}^{\mu \tau }<0.3\) and \({\mathcal {B}}^{\mu \tau } < 7.7\times 10^{-5}\) are considered. As we can see from Fig. 4, the information on \(A_{\mathrm {FB}}^{\tau \mu }\) helps to rule out a large part of the allowed space in the \(C_9^{\ell _1\ell _2}-C_{10}^{\ell _1\ell _2}\) plane.

Fig. 4
figure 4

Illustration of how the forward–backward asymmetry provides orthogonal constraints to the branching ratio of \(\Lambda _b\rightarrow \Lambda \mu ^- \tau ^+\). The shaded area present the allowed region for an upper limit \(A_{\mathrm{FB}}^{\mu \tau }<0.3\), compared to an upper limit for the branching ratio of \(7.7\times 10^{-5}\)

3.2 Explicit models

As mentioned, in many models that explain LFU violation, also LFV naturally occurs. Since our aim is not to perform a detailed analysis of all the observables in low-energy phenomenology, we choose to focus here on two specific models that explain the B anomalies. We choose two interesting solutions, which are the most favourite in the literature: the combination of the scalar leptoquarks \(S_1\) and \(S_3\) and the vector leptoquark \(U_1\). For these models we provide predictions for observables in \(\Lambda _b\rightarrow \Lambda \ell _1^-\ell _2^+\) decays.

The \({{\varvec{S}}}_{{\varvec{1}}}+{{\varvec{S}}}_{{\varvec{3}}}\) scalar leptoquarks scenario [44]

Here we focus on the \(S_1+S_3\) scenario, following the analysis in Ref. [44].Footnote 1 The main idea there is to apply the Froggatt–Nielsen mechanism [95], that explains the hierarchies of quark masses, as a power counting for NP operators, and thereby providing simultaneously an explanation for the B-anomalies and the flavour puzzle. Converting the formalism of Ref. [44] to the Wilson coefficients defined in Eq. (1), we find:

$$\begin{aligned} C_9^{\ell _1\ell _2} = - C_{10}^{\ell _1\ell _2}=\frac{v^2}{M^2} \frac{\pi }{\alpha _{\mathrm {em}}|V_{tb}V_{ts}^*|} |g_3|^2 {\tilde{S}}_{QL}^{3\ell _2}{\tilde{S}}_{QL}^{*2\ell _1} \end{aligned}$$
(22)

where M is the mass of the heavy scalar leptoquarks, \({\tilde{S}}_{QL}^{i\ell _i}\) is the spurion associated with the \(S_3\) scalar leptoquark and encodes the Froggatt–Nielsen power counting, and \(g_3\) is an overall coupling which is expected to be of \({\mathcal {O}}(1)\). Notice that the scalar leptoquark \(S_1\) does not contribute to \(b\rightarrow s\ell _1^-\ell _2^+\) transitions. With this, we find

$$\begin{aligned} C_9^{\mu \tau } = - C_{10}^{\mu \tau }=&\, -(0.41\pm 0.07) , \nonumber \\ C_9^{\tau \mu } = - C_{10}^{\tau \mu }=&\,(10\pm 2) . \end{aligned}$$
(23)

For the modes with electrons and muons in the final states, we find \(C_9^{e\mu }\propto 10^{-3}\) and a even lower value for \(C_9^{\mu e}\). Therefore, we conclude that the corresponding branching ratios are too small to be measured by any experiment in the near future.

Focusing then on the final states with muons and taus, using the Wilson coefficients in (23) and our results in Sect. 2 gives

$$\begin{aligned} \begin{aligned} {\mathcal {B}}^{\mu \tau } =\,&(7.1\pm 2.5 )\times 10^{-9}, \\ {\mathcal {B}}^{\tau \mu } =\,&(4.2\pm 1.7)\times 10^{-6}, \end{aligned} \end{aligned}$$
(24)

where the errors are dominated by the ones on the NP Wilson coefficients. Note that since this model predicts \(C_9^{\ell _1\ell _2}=-C_{10}^{\ell _1\ell _2}\), \(A_{\mathrm {FB}}^{\ell _1\ell _2}\) is independent from any Wilson coefficients and assumes the value in Eq. (21). From Fig. 4 we can conclude that the \(C_9^{\ell _1\ell _2}=-C_{10}^{\ell _1\ell _2}\) scenario would be excluded by the measurement of \(A_{\mathrm {FB}}^{\ell _1\ell _2}\). Hence, this stresses the importance of obtaining experimental constraints on this observable.

The \({{\varvec{U}}}_{{\varvec{1}}}\) vector leptoquark scenario [45]

Other interesting NP models are those with a vector leptoquark, usually denoted \(U_1\). In fact, this NP particle is the only one able to accommodate both classes of B anomalies on its own. Among the various possibilities available in the literature, we focus on [45], where the vector leptoquark is a massive state originating from the Spontaneous Symmetry Breaking of a gauge group larger than the SM one. As a consequence of the gauge representation of the \(U_1\) vector leptoquark, not only vector and axial vector couplings, but also scalar and pseudoscalar couplings are generated. In particular, the latter are very useful in explaining the large discrepancies in \(b\rightarrow c\tau {\bar{\nu }}\) data and as a consequence, generate sizeable \(b\rightarrow s\ell _1\ell _2\) interactions. Therefore, we expect very different signatures for the \(U_1\) model than the ones in the scalar leptoquark case. In Ref. [45] several cases are taken into account, where the flavour structure of the NP couplings has a \(U(2)^5\) flavour symmetry [96] or not, and where (pseudo-)scalar couplings are present or not. In the following we report results for the case in which no \(U(2)^5\) flavour symmetry is assumed. We note that using the scenario based on the \(U(2)^5\) flavour symmetry yields very similar results. We also note that given the flavour structure assumed in Ref. [45], the couplings of the vector leptoquark to electrons is zero, hence no effect is predicted for \(\Lambda _b\rightarrow \Lambda e^{\pm }\mu ^{\mp }\). In the notation of Ref. [45] we have

$$\begin{aligned} C_{9}^{\ell _1\ell _2} = - C_{10}^{\ell _1\ell _2}= & {} +\frac{2\pi }{\alpha _{\mathrm {em}}|V_{tb}V^*_{ts}|}C_U\beta _L^{2\ell _1}(\beta _L^{3\ell _2})^*, \nonumber \\ C_{S}^{\ell _1\ell _2} = - C_{P}^{\ell _1\ell _2}= & {} +\frac{4\pi }{\alpha _{\mathrm {em}}|V_{tb}V^*_{ts}|}C_U\beta _L^{2\ell _1}(\beta _R^{3\ell _2})^*, \end{aligned}$$
(25)

where \(C_U\) is a normalisation constant which contains the mass of the vector leptoquark normalised to the electroweak vacuum-expectation value and the gauge coupling of the leptoquark. The factor \(\beta _{L(R)}^{j\beta }\) represents the coupling in flavour space to left(right)-handed fermions. In the following we neglect the uncertainties on the fitted parameters obtained from [45] due to their large and asymmetric distributions. Either way, this scenario provides a useful benchmark that allows us to predict the size of LFV \(\Lambda _b\rightarrow \Lambda \) decays. We first look at the case \(\beta _R^{3\beta }=0\). We find

$$\begin{aligned} C_9^{\tau \mu } = -C_{10}^{\tau \mu }= & {} -5.93,\nonumber \\ C_9^{\mu \tau } = -C_{10}^{\mu \tau }= & {} +2.90. \end{aligned}$$
(26)

The predictions for \(A_{\mathrm {FB}}^{\ell _1\ell _2}\) in this case are the same as in Eq. (21). The corresponding integrated branching ratios are:

$$\begin{aligned} {\mathcal {B}}^{\tau \mu }= & {} 1.5 \times 10^{-6}\nonumber \\ {\mathcal {B}}^{\mu \tau }= & {} 3.6 \times 10^{-7} . \end{aligned}$$
(27)

In the case where \(\beta _R^{3\beta }\ne 0\), we find

$$\begin{aligned} C_9^{\tau \mu }= & {} -C_{10}^{\tau \mu } =-4.47,\quad C_S^{\tau \mu } = -C_{P}^{\tau \mu } =0,\nonumber \\ C_9^{\mu \tau }= & {} -C_{10}^{\mu \tau } =2.03,\quad C_S^{\mu \tau } = -C_{P}^{\mu \tau }= 4.06, \end{aligned}$$
(28)

which yields

$$\begin{aligned} {\mathcal {B}}^{\tau \mu }= & {} 8.5 \times 10^{-7}\quad \text {and}\quad A_{\mathrm {FB}}^{\tau \mu } = 0.14,\nonumber \\ {\mathcal {B}}^{\mu \tau }= & {} 5.3 \times 10^{-7}\quad \text {and} \quad A_{\mathrm {FB}}^{\mu \tau } = 0.12. \end{aligned}$$
(29)

Some comments are in order. In the scenario where \(C_{S(P)}^{\ell _1\ell _2}=0\), \(A_{\mathrm {FB}}^{\ell _1\ell _2}\) is independent from the Wilson coefficients and its value is given in Eq. (21). We find that \({\mathcal {B}}^{\tau \mu }>{\mathcal {B}}^{\mu \tau }\) due to a factor of two between the respective NP couplings. In the case where we have also \(C_{S}^{\mu \tau }=-C_{P}^{\mu \tau }\ne 0\), we find that \({\mathcal {B}}^{\mu \tau }\) is surprisingly small due to the negative interference between \(C_{S}^{\mu \tau }\) and \(C_{P}^{\mu \tau }\). On the other hand, \(A_{\mathrm {FB}}^{\mu \tau }\) is found to be smaller than \(A_{\mathrm {FB}}^{\tau \mu }\), hence providing a possible way to distinguish the different scenarios.

3.3 LHCb prospects

The results found in the above sections indicate that \(\Lambda _b\rightarrow \Lambda \ell _1\ell _2\) decays are very good probes of physics beyond the SM and provide in certain scenarios complementary bounds with respect to the ones from \({\bar{B}}\rightarrow \ell _1\ell _2\) and \(B^+\rightarrow K^+\ell _1\ell _2\) decays. Here we want to comment on the prospective for measurement of \(\Lambda _b\rightarrow \Lambda \ell _1\ell _2\) decays at the LHCb experiment.

If we consider measurement carried out with the same dataset, we expect for the measured yields:

$$\begin{aligned}&\frac{{\mathcal {N}}(\Lambda _b\rightarrow \Lambda (\rightarrow p\pi )\ell _1\ell _2)}{{\mathcal {N}}(B^+\rightarrow K^+\ell _1\ell _2)}\nonumber \\&\quad = \frac{{\mathcal {B}}(\Lambda _b\rightarrow \Lambda (\rightarrow p\pi )\ell _1\ell _2)| _{\mathrm {theory}}}{{\mathcal {B}}(B^+\rightarrow K^+\ell _1\ell _2)|_{\mathrm {theory}}} \frac{f_{\Lambda _b}}{f_{B^+}}r_{\Lambda _b/B^+}, \end{aligned}$$
(30)

where \(f_{\Lambda _b}/f_{B^+}\) is the ratios of the fragmentation functions for the \(\Lambda _b\) and the \(B^+\) modes, respectively, and \(r_{\Lambda _b/B^+}\) is a correction factor due to different reconstruction efficiencies. In Ref. [97], the ratio \(f_{\Lambda _b}/{(f_u+f_d)}\) is measured. Using isospin relations, we can write \(f_{\Lambda _b}/{f_{B^+}} = 2 f_{\Lambda _b}/{(f_u+f_d)} = (0.518 \pm 0.036)\). The ratio of the predicted values of the theoretical branching ratios depends on the NP model and final state leptons. However, as we noted already in Sect. 3.1, the branching ratios for the baryon and the meson case are very similar in size: therefore, for an order of magnitude estimate, we consider them to be equal. The last piece of information needed is the ratio \(r_{\Lambda _b/B^+}\), that is difficult to estimate without a thorough simulation of the LHCb detector. However, in order to give an estimate, we use the information in Refs. [6, 98], that are based on the same integrated luminosity. From these papers we extract

$$\begin{aligned} \frac{{\mathcal {N}}(\Lambda _b\rightarrow \Lambda (\rightarrow p \pi )\mu ^-\mu ^+)}{{\mathcal {N}}(B^+\rightarrow K^+\mu ^-\mu ^+) } \approx 0.31. \end{aligned}$$
(31)

This means that we expect the efficiency for the reconstruction of the \(\Lambda _b\) to be roughly 1.67 times less than that of the \(B^+\), when also taking into account the fragmentation fractions effect. Hence we set \(r_{\Lambda _b/B^+}=1.67\). We expect that all the other correction factors due to the reconstruction of the leptons in the final state cancel out since we are comparing the same leptonic final states in both decays. This yields

$$\begin{aligned} {\mathcal {B}}(\Lambda _b\rightarrow \Lambda (\rightarrow p\pi )\ell _1\ell _2)\approx 1.67 \frac{f_{\Lambda _b}}{f_{B^+}}{\mathcal {B}}(B^+\rightarrow K^+\ell _1\ell _2). \end{aligned}$$
(32)

Using the current upper limit on \({\mathcal {B}}(B^+\rightarrow K^+\tau ^+\mu ^-)\), we thus expect that LHCb can reach the following sensitivity:

$$\begin{aligned} {\mathcal {B}}(\Lambda _b\rightarrow \Lambda (\rightarrow p\pi )\mu ^-\tau ^+)\lesssim 6.5\times 10^{-5}, \end{aligned}$$
(33)

for Run 1 and 2 datasets. In the above estimate, we have not included any correction for the trigger efficiency, which can be different for the baryonic and mesonic mode. The estimate in Eq. (33) can be compared to the model dependent and model independent bounds on \({\mathcal {B}}(\Lambda _b\rightarrow \Lambda \tau ^+\mu ^-)\) found in the previous sections. In particular, the expected upper bound from LHCb would already give better constraints than the corresponding ones from the mesonic decays, as illustrated in Fig. 3. We also stress that future runs will improve the upper limit in Eq. (33) of at least a factor of roughly two with Run 3 and a factor of three with further runs [66].

4 Conclusions

In this paper, we present the first full analysis of \(\Lambda _b\rightarrow \Lambda \ell _1^{-} \ell _2^{+}\) lepton flavour violating (LFV) decays in terms of possible new physics operators. The main results of this paper are Eqs. (10)–(12), where the coefficients of the angular distributions for \(\Lambda _b\rightarrow \Lambda \ell _1^{-} \ell _2^+\) decays are given. We study the interplay between the baryonic and mesonic searches for LFV, where for the latter upper limits are already available. We convert these upper limits into constraints on the branching ratio and forward–backward asymmetry for \(\Lambda _b\rightarrow \Lambda \ell _1^{-} \ell _2^{+}\) decays. We find that the \(\Lambda _b\rightarrow \Lambda \ell _1^{-} \ell _2^{+}\) decays provide different constraints on the new physics Wilson coefficients than \({\bar{B}}_s\rightarrow \ell _1^{-} \ell _2^{+}\) and \(B^+\rightarrow K^+ \ell _1^{-} \ell _2^{+}\) decays, and have the potential to reduce the allowed parameter space for new physics models. We then analyse quantitatively the size of \(\Lambda _b\rightarrow \Lambda \ell _1^{-} \ell _2^{+}\) decays in specific scenarios that can address B anomalies, using as a Refs. [44, 45]. Our findings indicate that the predicted branching ratio for \(\Lambda _b\rightarrow \Lambda \ell _1^{-} \ell _2^{+}\) for these scenarios are such that they can further constrain the new physics couplings. As a final prospective, we estimate the reach of LHCb for \(\Lambda _b\rightarrow \Lambda \ell _1^{-} \ell _2^{+}\) decays, finding that an upper limit of \({\mathcal {B}}(\Lambda _b\rightarrow \Lambda \mu ^-\tau ^+)\lesssim 6.5\cdot 10^{-5}\) can be reached with Run 1 and Run 2 data.

Table 6 Correlation matrix for the \(\Lambda _b\rightarrow \Lambda \mu ^-e ^+\) parameters