Hunting for $B^+\rightarrow K^+ \tau ^+\tau ^-$ imprints on the $B^+ \rightarrow K^+ \mu ^+\mu ^-$ dimuon spectrum

Abstract

We investigate the possibility of indirectly constraining the $B^{+}\rightarrow K^{+}\tau ^+\tau ^-$ decay rate using precise data on the $B^{+}\rightarrow K^{+}\mu ^+\mu ^-$ dimuon spectrum. To this end, we estimate the distortion of the spectrum induced by the $B^{+}\rightarrow K^{+}\tau ^+\tau ^-\rightarrow K^{+} \mu ^+\mu ^-$ re-scattering process, and propose a method to simultaneously constrain this (non-standard) contribution and the long-distance effects associated to hadronic intermediate states. The latter are constrained using the analytic properties of the amplitude combined with data and perturbative calculations. Finally, we estimate the sensitivity expected at the LHCb experiment with present and future datasets. We find that constraints on the branching fraction of $O(10^{-3})$, competitive with current direct bounds, can be achieved with the current dataset, while bounds of $O(10^{-4})$ could be obtained with the LHCb upgrade-II luminosity.

1 Introduction

In recent years, discrepancies between the observed values and the Standard Model (SM) predictions of the lepton-flavour universality (LFU) ratios $R_{D^{(*)}}$ [1,2,3,4,5] and $R_{K^{(*)}}$ [6,7,8,9], characterizing the semileptonic transitions $b\rightarrow cl\nu$ and $b\rightarrow sll$, have sparked great interest. The pattern of anomalies seems to point to intriguing new-physics (NP) scenarios, with possible connections to the SM flavour puzzle. A large class of NP models proposed to explain these hints of physics beyond the SM, and in particular those aiming for a combined explanation of the $R_{K^{(*)}}$ and $R_{D^{(*)}}$ anomalies, imply dominant couplings to third-generation fermions, which should also enter other semileptonic b-quark decays.

A general expectation, confirmed by many explicit NP constructions, is that of a large enhancement of $b\rightarrow s\tau ^+\tau ^-$ transitions (see e.g. [10,11,12,13,14,15,16,17]). While flavour-changing neutral-current (FCNC) decays with muon and electron pairs have been observed both at the exclusive and at the inclusive level, probing rare decays with a $\tau ^+\tau ^-$ pair in the final state is experimentally very challenging. The current experimental limits for all processes mediated by the $b\rightarrow s\tau ^+\tau ^-$ amplitude are still very far from the corresponding SM predictions [18, 19], leaving the NP expectation of possible large enhancements unchallenged.

In this work we investigate the possibility of indirectly constraining the $b\rightarrow s\tau ^+\tau ^-$ amplitude via its imprint on the $B^+\rightarrow K^+ \mu ^+\mu ^-$ dimuon spectrum. In presence of a large NP enhancement, the $b\rightarrow s\tau ^+\tau ^-$ amplitude would induce a distinctive distortion of the $B^+\rightarrow K^+ \mu ^+\mu ^-$ spectrum via the (QED-induced) re-scattering process $B^+\rightarrow K^+ \tau ^+\tau ^- \rightarrow K^+ \mu ^+\mu ^-$ [10]. The latter has a discontinuity at $q^2=4m_\tau ^2$ ($q^2\equiv m^2_{\mu \mu }$), namely at the threshold where the tau leptons can be produced on-shell. This gives rise to a “cusp” in the dimuon-invariant mass spectrum, which could in principle be detected with sufficient experimental precision. More generally, the lightness of the $\tau$-leptons implies a well-defined deformation of the $B^+\rightarrow K^+ \mu ^+\mu ^-$ spectrum, which is determined only by the analytic properties of the re-scattering amplitude.

It should be stressed that the phenomenon we are considering here is different from the QED mixing between dimension-six FCNC operators with different lepton species analysed in Ref. [20]. If NP is heavy and the $b\rightarrow s\tau ^+\tau ^-$ amplitude is strongly enhanced, the operator mixing can give rise to sizable modifications of the Wilson coefficients of the dimension-six effective Hamiltonian relevant to $b\rightarrow sl ^+l ^-$ decays ($l =e,\mu$). However, this phenomenon cannot be distinguished in a model-independent way from other NP effects of short-distance origin (at least using low-energy data only). On the contrary, the non-local effect we are interested in can be unambiguously attributed to the re-scattering of light intermediate states characterised by the tau mass, hence it can be translated into a model-independent constraint on the $B^+\rightarrow K^+ \tau ^+\tau ^-$ amplitude.

The main difficulty in extracting such bound is obtaining a reliable description of the $B^+\rightarrow K^+ \mu ^+\mu ^-$ dimuon spectrum within the SM, or better in the limit where the $\tau ^+\tau ^- \rightarrow \mu ^+\mu ^-$ re-scattering is negligible. This is non trivial, given that the $B^+\rightarrow K^+ l ^+ l ^-$ spectrum is plagued by theoretical uncertainties originating from $B\rightarrow K$ form factors and hadronic long-distance contributions. While the former are smooth functions in the $q^2$ region of interest and can be well described using lattice QCD [21, 22] and/or light-cone sum rules [23], long-distance effects induced by intermediate hadronic states, such as the charmonium resonances, are more problematic. They are genuine non-perturbative effects and introduce physical discontinuities below and above the $q^2=4m_\tau ^2$ threshold. Far from the resonance region, these effects can be estimated using perturbative constraints derived at $q^2<0$, with $|q^2| \gg \Lambda _\mathrm{QCD}^2$, combined with a $\Lambda _\mathrm{QCD}^2/q^2$ or $\Lambda _\mathrm{QCD}^2/m_c^2$ expansion to incorporate the leading non-perturbative corrections [24, 25]. However, this approach is not suitable for our purpose, which requires a reliable description of the whole spectrum, and in particular of the resonance region. To achieve this goal, we adopt a data-driven approach which takes full advantage of the known analytic properties of the amplitude: knowing the precise location of all one- and two-particle hadronic thresholds, we use subtracted dispersion relations to describe the $q^2$-dependence of the whole spectrum in terms of a series of ($q^2$-independent) hadronic parameters, which are fitted from data. This method can be considered a generalisation of the approaches proposed in Ref. [26] and, to some extent, in Refs. [27,28,29], with a few key differences, the most notable ones being the use of subtracted dispersion relations and the explicit inclusion of two-particle thresholds. To reduce the number of independent free parameters, perturbative constraints derived from the low-$q^2$ region are also implemented. Proceeding this way we obtain a description of the spectrum that is flexible enough to extract the non-perturbative parameters characterising the various hadronic thresholds from data, but retains a significant predictive power in the smooth region within and below the two narrow charmonium states, allowing us to set useful constraints on the $B^+\rightarrow K^+ \tau ^+\tau ^- \rightarrow K \mu ^+\mu ^-$ re-scattering.

The method we propose is particularly well suited for the LHCb experiment, which has already collected a large sample of $B^+\rightarrow K^+ \mu ^+\mu ^-$ events and has an excellent resolution in the dimuon spectrum [30]. In order to estimate the sensitivity of LHCb in view of the full run II dataset, we generate pseudo-experiments based on the yields and amplitudes obtained in Ref. [30], and calculate the expected limit under the background-only hypothesis using the CLs method [34].

The paper is organised as follows: in Sect. 2 we introduce the theoretical framework necessary to describe the $B^+\rightarrow K^+ \mu ^+\mu ^-$ dimuon spectrum within and beyond the SM, separating short-distance contributions (Sect. 2.1), long-distance contributions due to intermediate hadronic states (Sect. 2.3), and long-distance contributions due to the $\tau ^+\tau ^- \rightarrow \mu ^+\mu ^-$ re-scattering (Sect. 2.4). The analysis of the LHCb sensitivity is presented in Sect. 3. The results are summarised in the Conclusions.

2 Theoretical framework

2.1 Effective Hamiltonian and differential decay rate

The dimension-six effective Langrangian describing $b \rightarrow s ll$ transitions, renormalized at low energies [$\mu =O(m_b)$], can be decomposed as

$$\begin{aligned} {\mathcal {L}}_\mathrm {eff} = \frac{4G_F}{\sqrt{2}}V_{tb}V_{ts}^*\sum _{i}{\mathcal {C}}_i(\mu )\, {\mathcal {O}}_i\,, \end{aligned}$$

(1)

where the leading FCNC effective operators are defined as

$$\begin{aligned} \begin{aligned} {\mathcal {O}}_7&= \frac{e}{16\pi ^2}\left( {\bar{s}}\sigma _{\mu \nu }(m_b P_R+ m_s P_L)b \right) F^{\mu \nu }\,, \\ {\mathcal {O}}^l_9&= \frac{e^2}{16\pi ^2}({\bar{s}} \gamma _\mu P_L b)({\bar{l}} \gamma ^\mu l)\,, \\ {\mathcal {O}}^l_{10}&= \frac{e^2}{16\pi ^2}({\bar{s}} \gamma _\mu P_L b)({\bar{l}} \gamma ^\mu \gamma _5 l)\,, \end{aligned} \end{aligned}$$

(2)

and the most relevant four-quark operators ($q=u,c$) as

$$\begin{aligned} \begin{aligned} {\mathcal {O}}_1^q&= ({\bar{s}} \gamma _\mu P_L q)({\bar{q}}\gamma ^\mu P_L b)\,, \\ {\mathcal {O}}_2^q&= ({\bar{s}}^\alpha \gamma _\mu P_L q^\beta )(\bar{q}^\beta \gamma ^\mu P_L b^\alpha ). \end{aligned} \end{aligned}$$

(3)

Within the class of models we are considering, all relevant NP effects are encoded in the values of the Wilson coefficients ${\mathcal {C}}^l_{7,9,10}$. Given the normalisation in Eq. (1), ${\mathcal {C}}^l_{7,9,10}$ and ${\mathcal {C}}^c_{1,2}$ are real and O(1) within the SM, whereas ${\mathcal {C}}^u_{1,2}= (V_{ub} V_{us}^*/V_{tb} V_{ts}^{*}) \times O(1)$ (see Ref. [26] for the precise values of the Wilson coefficients and the complete basis of operators).

The matrix elements $\langle K^+\mu ^+\mu ^- | {\mathcal {O}}_i | B^+ \rangle$ are non-vanishing at the tree level only in the case of the FCNC operators (with $l=\mu$). Considering only the contribution of the FCNC operators, which can be expressed in terms of the $B\rightarrow K$ form factors, the $B^+\rightarrow K^+\mu ^+\mu ^-$ decay rate can be written as:

$$\begin{aligned}&\left. \frac{d\Gamma }{dq^2} \right| _{{\mathcal {C}}_{7,9,10}} \!\!\! = \frac{\alpha _\mathrm {em}^2G_F^2 |V_{tb}V_{ts}^*|^2}{128\,\pi ^5} \kappa (q^2) \beta (q^2) \nonumber \\&\quad \left\{ \frac{2}{3} \kappa ^2(q^2) \beta ^2(q^2) \left| {\mathcal {C}}^\mu _{10}f_+(q^2)\right| ^2 + \frac{m_\mu ^2(m_B^2-m_K^2)^2}{q^2\, m_B^2}\left| {\mathcal {C}}^\mu _{10}f_0(q^2)\right| ^2\right. \nonumber \\&\qquad \left. + \kappa ^2(q^2)\left[ 1-\frac{1}{3}\beta ^2(q^2) \right] \, \left| {\mathcal {C}}^\mu _9 f_+(q^2)\right. \right. \nonumber \\&\qquad \left. \left. +2{\mathcal {C}}_7 \frac{m_b+m_s}{m_B +m_K} f_T(q^2)\right| ^2 \right\} \,, \end{aligned}$$

(4)

where $\kappa (q^2)=\lambda ^{1/2}(m_B^2,m_K^2,q^2)/2m_B$ is the kaon momentum in the B-meson rest frame, $\beta (q^2)=\sqrt{1-4m_\mu ^2/q^2}$, and $f_i(q^2)$ with $i=+,0,T$ are the vector, scalar and tensor $B\rightarrow K$ form factors.

2.2 Non-local contributions: general considerations

The non-local contributions generated by the non-leptonic operators in ${\mathcal {L}}_\mathrm {eff}$ and by the operator ${\mathcal {O}}^\tau _9$ can be encoded in Eq. (4) by replacing ${\mathcal {C}}^\mu _9$ with a $q^2$-dependent function:

$$\begin{aligned} {\mathcal {C}}^\mu _9\ \rightarrow \ {\mathcal {C}}_9^{\mu ,\mathrm {eff}}(q^2)= & {} {\mathcal {C}}^\mu _9+ Y_{c \bar{c}} (q^2) + Y_\mathrm{light} (q^2) \nonumber \\&+ Y_{\tau \bar{\tau }}(q^2)\,, \end{aligned}$$

(5)

where $Y_{{\mathcal {I}}} (q^2)$ denotes the non-local contributions corresponding to the intermediate state ${\mathcal {I}}$, which can annihilate into a dimuon pair via a single-photon exchange.

The functions $Y_{c {\bar{c}}} (q^2)$ and $Y_\mathrm{light} (q^2)$ encode non-perturbative hadronic contributions, which cannot be estimated reliably in perturbation theory, at least in a large fraction of the accessible $q^2$ spectrum. Adopting a notation similar to that of Ref. [26], we can express $Y_{c {\bar{c}}} (q^2)$ as

$$\begin{aligned} Y_{c {\bar{c}}} (q^2) = \frac{16\pi ^2 }{ f_+(q^2) } {\mathcal {H}}_{c \bar{c}}^{(BK)}(q^2), \end{aligned}$$

(6)

where ${\mathcal {H}}_{c {\bar{c}}}^{(BK)}(q^2)$ is defined by the gauge-invariant decomposition of the following non-local hadronic matrix element

$$\begin{aligned}&i \int d^4 x e^{iq\cdot x} \Bigg \langle K(p) | T\left\{ j_\mu ^\mathrm{em}(x), \sum _{i=1,2} {\mathcal {C}}^c_i {\mathcal {O}}^c_i \right\} | B(p+q) \Bigg \rangle \nonumber \\&\quad = [ (p\cdot q) q_\mu - q^2 p_\mu ] {\mathcal {H}}_{c \bar{c}}^{(BK)}(q^2)~, \end{aligned}$$

(7)

with $j_\mu ^\mathrm{em}= \sum _{q} Q_q {\bar{q}} \gamma ^\mu q$. The function $Y_\mathrm{light} (q^2)$, containing the contribution of the subleading non-leptonic operators in ${\mathcal {L}}_\mathrm {eff}$, is defined in a similar way via the replacement

$$\begin{aligned} \sum _{i=1,2} {\mathcal {C}}^c_i {\mathcal {O}}^c_i ~\rightarrow \sum _{i=3-6,8} {\mathcal {C}}_i {\mathcal {O}}_i + \sum _{i=1,2} {\mathcal {C}}^u_i {\mathcal {O}}^u_i~. \end{aligned}$$

(8)

Our main strategy is to write the non-perturbative functions $Y_{c {\bar{c}}} (q^2)$ and $Y_\mathrm{light} (q^2)$ using hadronic dispersion relations. More precisely, for the leading charm contribution we consider one- (1P) and two-particle (2P) intermediate states (see Fig. 1), using dispersion relations subtracted at $q^2=0$, while for the subleading $Y_\mathrm{light} (q^2)$ function we consider only one-particle intermediate states and use unsubtracted dispersion relations. We stress that these dispersion relations, and the corresponding decomposition of ${\mathcal {H}}_{c {\bar{c}}}^{(BK)}(q^2)$, are not exhaustive of all the discontinuities of the four-point function in (7). However, our goal is not to determine completely ${\mathcal {H}}_{c {\bar{c}}}^{(BK)}(q^2)$ via dispersion relations, but only to describe its functional dependence with respect to the variable $q^2$.

Given these considerations, ${\mathcal {C}}_9^{\mu ,\mathrm {eff}}(q^2)$ in Eq. (5) is finally decomposed according to

$$\begin{aligned} {\mathcal {C}}_9^{\mu ,\mathrm {eff}}(q^2)= & {} {\mathcal {C}}^\mu _9 + Y^{(0)}_{c{\bar{c}}} + \Delta Y^\mathrm{1P}_{c {\bar{c}}} (q^2) + \Delta Y^\mathrm{2P}_{c {\bar{c}}} (q^2) \nonumber \\&+ Y^\mathrm{1P}_\mathrm{light} (q^2) + Y_{\tau \bar{\tau }} (q^2)~, \end{aligned}$$

(9)

with $\Delta Y^\mathrm{1P}_{c {\bar{c}}} (0) = \Delta Y^\mathrm{2P}_{c \bar{c}} (0) = 0$. In the next section we analyse the structure of $\Delta Y^\mathrm{1P}_{c {\bar{c}}} (q^2)$, $\Delta Y^\mathrm{2P}_{c {\bar{c}}} (q^2)$, and $Y^\mathrm{1P}_\mathrm{light} (q^2)$ in detail. The expression of $Y_{\tau \bar{\tau }} (q^2)$, which is the only term in Eq. (9) that can be fully evaluated in perturbation theory, is given in Sect. 2.4.

Fig. 1

Diagrammatic representations of the long-distance contributions to $C_9^{\mu ,\mathrm {eff}}$. The left-hand side depicts the exchange of a single vector resonance. The graph on the right-hand side shows the contribution from two-particle intermediate states

2.3 Long-distance hadronic contributions

The general structure of the subtracted dispersion relation used to determine $\Delta Y_{c {\bar{c}}} (q^2)$ is

$$\begin{aligned} \Delta Y_{c {\bar{c}}} (q^2)= & {} \frac{16 \pi q^2}{ f_+(q^2)} \int _{m^2_{J/\Psi }}^\infty ds \frac{ 1 }{s (s - q^2) } ~ \frac{1}{2i} \mathrm{Disc}\left[ {\mathcal {H}}_{c {\bar{c}}}^{(BK)}(s) \right] \nonumber \\\equiv & {} \frac{q^2}{\pi } \int _{m^2_{J/\Psi }}^\infty ds \frac{ \rho _{c{\bar{c}}}(s) }{s (s - q^2) }. \end{aligned}$$

(10)

The function $\rho _{c{\bar{c}}}(s)$ is the spectral-density function describing the hadronic states ${\mathcal {I}}_{c{\bar{c}}}$, characterized by valence charm-quarks and invariant mass s, contributing as real intermediate states in the re-scattering $B\rightarrow K {\mathcal {I}}_{c{\bar{c}}} \rightarrow K \mu ^+\mu ^-$. As noted before, we decompose $\rho _{c{\bar{c}}}(s)$ into one- and two-particle intermediates states, $\rho _{c{\bar{c}}}(s) = \rho ^\mathrm{1P}_{c{\bar{c}}}(s) + \rho ^\mathrm{2P}_{c{\bar{c}}}(s)$,

$$\begin{aligned} \rho ^\mathrm{1P}_{c{\bar{c}}}(s)\propto & {} \sum _j {\mathcal {A}}( B\rightarrow K V^0_j) {\mathcal {A}}( V^0_j \rightarrow \mu ^+\mu ^-) \delta ( s - m_j)~, \end{aligned}$$

(11)

$$\begin{aligned} \rho ^\mathrm{2P}_{c{\bar{c}}}(s)\propto & {} \sum _j \int d p_j^2 ~\delta ( s- p_j^2 ) \int \frac{d^3 \mathbf {p}_{j_1} d^3 \mathbf {p}_{j_2} }{ 16 \pi ^2 E_{j_1} E_{j_2} } {\mathcal {A}}( B\rightarrow K M^+_{j_1} M^-_{j_2} ) \nonumber \\&\times {\mathcal {A}}( M^+_{j_1} M^-_{j_2} \rightarrow \mu ^+\mu ^-) \delta ^{(4)}( p_j - p_{j_1} -p_{j_2}) ~, \end{aligned}$$

(12)

neglecting the phase-space suppressed contribution with three or more particles.

2.3.1 Charmonium resonances

For the sake of simplicity, in Eq. (11) we have treated the single-particle states as infinitely narrow resonances. The effect of finite widths can be incorporated via Breit–Wigner functions, yielding

$$\begin{aligned} \Delta Y_{c{\bar{c}}}^\mathrm {1P}(q^2)= & {} \sum _{j = \Psi (1S), \ldots ,\Psi (4415) } \eta _j\, e^{i\delta _j} \frac{q^2}{m_j^2} A^\mathrm {res} _j (q^2)~, \nonumber \\ A^\mathrm {res} _j (s)= & {} \frac{m_j\Gamma _j}{(m_j^2-s)-im_j \Gamma _j}\,~, \end{aligned}$$

(13)

where the sum runs over all the charmonium vector resonances in the accessible kinematical range. Here $\eta _j$ and $\delta _j$ are real parameters which must be determined from data, similarly to what has been performed by the LHCb collaboration in [30].

We stress that our scope is not to compute the $\eta _j$ employing a facotrization hypothesis for the hadronic matrix element in (7), combined with experimental data on $c{\bar{c}} \rightarrow e^+e^-$, as originally proposed in [31]. As clearly shown in [27], this approach leads to rather inconsistent results. We refrain from any attempt to compute the $\eta _j$, which are genuine non-perturbative hadronic matrix elements: we treat them as free parameter which need to be determined by data.

In principle, both $\eta _j$ and $\delta _j$ are $q^2$-dependent functions. However, if this dependence is smooth around the resonance poles (or the two-particle thresholds, for the analog parameters that we define in Sect. 2.3.2), we can treat them as constant terms for the purpose of our analysis. On general grounds, we expect the characteristic scale determining the $q^2$-variation of the $\{\eta _j, \delta _j\}$ appearing in $\Delta Y_{c{\bar{c}}}^\mathrm {1P}$ to be twice the charm mass. This $q^2$-dependence is indeed related to the invariant mass of the hadronic intermediate states that can mix into the charmonia via re-scattering processes. Incidentally, we note that this is also what one would infer from a perturbative estimated of the $\eta _j$ using the factorization hypothesis. In view of this argument, we believe that it is a good approximation to treat the $\{\eta _j, \delta _j\}$ as constant terms in our analysis. As stated above, our goal is not to compute these parameters but only to fit them from data in order to have a sufficiently general description of the long-distance part the amplitude, able to reproduce all the known discontinuities related to hadronic intermediate states.

The fitted $\eta _j$’s can be put in one-to-one correspondence with the product of the $B^+ \rightarrow K^+ V^0_j$ and $\ V^0_j \rightarrow \mu ^+\mu ^-$ branching fractions via

$$\begin{aligned}&{\mathcal {B}}( B^+ \rightarrow K^+ V^0_j) \times {\mathcal {B}}( V^0_j \rightarrow \mu ^+\mu ^-)\nonumber \\&\quad = \tau _{B^+} \frac{G_F^2\alpha ^2|V_{tb}V_{ts}^*|^2}{128\pi ^5} \int \limits _{4 m^2_\mu }^{(m_B - m_K)^{2}} dq^2 \kappa (q^2)^3 \nonumber \\&\qquad \times \left[ \beta (q^2)- \frac{1}{3}\beta (q^2)^3 \right] \left| f_+(q^2) \right| ^2 \left| \eta _j \right| ^2 \left| \frac{q^2}{m_j^2} A_j^\mathrm{res} (q^2) \right| ^2 \,.\nonumber \\ \end{aligned}$$

(14)

The expression (13) differs from the decomposition adopted in Ref. [30] by the $q^2/m^2_j$ term, which arises from the subtraction procedure in the dispersion relation. On the one hand, the use of subtracted dispersion relations for the charm contribution is necessary to ensure the convergence of the integral in the two-particle intermediate states (see Sect. 2.3.2). On the other hand, choosing the subtraction point at $q^2=0$ allows us to decouple the determination of the resonance parameters of the spectrum from the overall normalisation of the rate, and hence from the determination of ${\mathcal {C}}_9^\mu$ from data. The price to pay is the appearance of the undetermined constant term $Y^{(0)}_{c{\bar{c}}}= Y_{c {\bar{c}}} (0)$ in Eq. (9). This term plays no role in the description of the dimuon spectrum, but is relevant for the extraction of the value of ${\mathcal {C}}_9^\mu$. To this purpose, we note that the estimate presented in Ref. [26], which is based on a $\Lambda ^2/m_c^2$ expansion and also takes next-to-leading $O(\alpha _s)$ corrections on the pure partonic result into account (see Sect. 2.3.4), yields

$$\begin{aligned} Y^{(0)}_{c{\bar{c}}} \approx -0.10 \pm 0.05~, \end{aligned}$$

(15)

which is about $-(2\pm 1)\%$ of ${\mathcal {C}}_9^{\mu , \mathrm{SM}} \approx 4.23$.

2.3.2 Two-particle intermediate states

Proceeding in a similar way, we can decompose the two-particle contributions as

$$\begin{aligned} \Delta Y^\mathrm {2P}_{c{\bar{c}}}(q^2)&=\sum _{j} \eta _j e^{i\delta _j} A_j^\mathrm {2P}(q^2)\,, \nonumber \\ A_j^\mathrm {2P}(q^2)&= \frac{q^2}{\pi }\int _{s_0^j}^\infty \frac{d s}{s}\frac{ {\hat{\rho }}_j (s)}{(s-q^2)}\,, \end{aligned}$$

(16)

where ${\hat{\rho }}_j(s)$ are normalised spectral densities for the two-body intermediate states characterised by the threshold $s_0^j = (m_{j_1}+m_{j_2} )^2$.

While we do not have a precise estimate of these spectral densities at generic kinematical points, an excellent description of their behaviour around the respective thresholds is obtained by approximating them with powers of the Källén function, with an exponent determined by the lowest partial wave allowed in the $B^+ \rightarrow K^+ M_1 M_2 \rightarrow K^+ \mu ^+\mu ^-$ re-scattering. This is because higher-order partial waves, characterised by higher powers of the Källén function, are both phase-space suppressed and, most importantly, give rise to a less singular behaviour at the threshold. From angular momentum conservation we can then determine the leading partial wave and obtain the following estimates for the normalised spectral densities of the two-particle intermediate states of lowest mass:

$$\begin{aligned} {\hat{\rho }}_{DD}(s)&= \left( 1-\frac{4m_D^2}{s}\right) ^{3/2}\,, \nonumber \\ {\hat{\rho }}_{D^*D^*}(s)&= \left( 1-\frac{4m_{D^*}^2}{s}\right) ^{3/2}\,, \nonumber \\ {\hat{\rho }}_{DD^*}(s)&= \left( 1-\frac{4m_{{\bar{D}}}^2}{s}\right) ^{1/2}. \end{aligned}$$

(17)

In the case of the $DD^*$ intermediate state we have replaced the complete expression depending on both masses with a simplified one depending only on $m_{{\bar{D}}}= (m_D+m_{D^*})/2$, which provides an excellent approximation. With these estimates in place we find:

$$\begin{aligned} \Delta Y^{2\mathrm {P}}_{c{\bar{c}}} (q^2)&= \eta _{{\bar{D}}} e^{i\delta _{{\bar{D}}}} h_{S} \left( m_{{\bar{D}}}, q^2 \right) \nonumber \\&\quad + \sum _{j=D, D^*} \eta _{j} e^{i\delta _j} h_{P} \left( m_{j}, q^2 \right) \,, \end{aligned}$$

(18)

with

$$\begin{aligned} \begin{aligned} h_P \left( m, q^2 \right)&= \frac{2}{3}+ \left( 1 -\frac{4 m^2}{q^2} \right) h_S \left( m, q^2 \right) \,, \\ h_S \left( m, q^2 \right)&= 2- G\left( 1 -\frac{4 m^2}{q^2} \right) \,, \end{aligned} \end{aligned}$$

(19)

and

$$\begin{aligned} G(y)&= \sqrt{|y|} \left\{ \Theta (y) \left[ \ln \left( \frac{1+\sqrt{y}}{ 1-\sqrt{y} } \right) -i \pi \right] \right. \nonumber \\&\quad \left. + 2~\Theta (-y) \arctan \left( \frac{1}{\sqrt{-y}} \right) \right\} \,. \end{aligned}$$

(20)

Fig. 2

Real (solid) and imaginary (dashed) parts of the normalised hadronic two-particle contributions to $Y_{c\bar{c}}(q^2)$, as defined in Eq. (16)

It is worth noting that, while the lowest threshold is at $q^2=4m_D^2$, the contribution from the $DD^*$ intermediate state is the only one which can occur in the S-wave, corresponding to a singular (square-root) behaviour at the threshold (see Fig. 2).

2.3.3 Light resonances

The remaining hadronic contribution we need to estimate is $Y_\mathrm{light} (q^2)$, defined by Eqs. (6) and (7) via the replacement (8). The Wilson coefficients of the effective operators appearing in ${\mathcal {H}}_\mathrm{light}^{(BK)}(q^2)$ are either loop- or CKM-suppressed. As a result, we can limit ourselves to include only one-particle hadronic intermediate states. In principle, such operators describe transitions also to states with valence charm quarks; however, since we fit the hadronic coefficients $\eta _j$ from data, these terms are naturally absorbed in the $\eta _j$ appearing in $\Delta Y_{c{\bar{c}}} (q^2)$. We are thus left only with vector resonances containing light valence quarks. Among them, we can further restrict the attention to the $\rho$, $\omega$, and $\phi$ resonances, since the leptonic decay rates of the heavier states are very small.

There is no clear advantage in using subtracted vs. unsubtracted dispersion relations in describing the contributions of the light vector resonances. The convergence of the dispersive integrals does not pose a problem, and the subtraction at $q^2=0$ is not particularly useful since the light-quark contributions are in a non-perturbative regime at $q^2=0$. However, when fitting data, the subtraction at $q^2=0$ retains the advantage of decoupling the determination of the spectrum from that of the Wilson coefficient. As default option, we adopt unsubtracted dispersion relations. As discussed in Sect. 2.3.5, checking the stability of the result using subtracted vs. unsubtracted dispersion relations for the light vector resonances provides an estimate of the “model error” of the proposed approach.

Given these considerations, we decompose $Y^\mathrm{1P}_\mathrm{light} (q^2)$ as

$$\begin{aligned} Y_\mathrm{light}^\mathrm {1P}(q^2) = \sum _{j = \rho , \omega , \phi } \eta _j\, e^{i\delta _j} A^\mathrm {res} _j (q^2)~, \end{aligned}$$

(21)

in perfect analogy with the decomposition adopted in Ref. [30] for these light states.

2.3.4 Theoretical constraints on the hadronic parameters

The hadronic decompositions in Eqs. (13), (18) and (21) contain 12 free complex parameters: 6 in $\Delta Y_{c{\bar{c}}}^\mathrm {1P}(q^2)$, 3 in $\Delta Y_{c\bar{c}}^\mathrm {2P}(q^2)$, and 3 in $Y_\mathrm{light}^\mathrm {1P}(q^2)$. In principle, since they correspond to different functional forms, they could all be fitted from data. In practice however, an unconstrained fit would leave significant degeneracies in the parameter space. It is therefore useful to restrict the variability of such parameters using theoretical constraints. In the following we discuss three conservative conditions which can be imposed using perturbative arguments.

I. Constraint on the slope of $\Delta Y_{c {\bar{c}}} (q^2)$ at $q^2=0$.

The lowest-order perturbative estimate of $\Delta Y_{c {\bar{c}}} (q^2)$ is obtained by factorising the matrix element $\langle K(p) | {\bar{s}} \gamma ^\mu b | B(p+q) \rangle$ in Eq. (7) and computing the charm-loop at $O(\alpha ^0_s)$:

$$\begin{aligned} \Delta Y^\mathrm{pert}_{c {\bar{c}}} (q^2)= & {} 2 \left( {\mathcal {C}}_2+\frac{1}{3} {\mathcal {C}}_1\right) \nonumber \\&\times Q_c \times q^2 \int _{4 m_c^2}^\infty ds \frac{\sqrt{1-\frac{4 m_c^2}{s} }\left( 1+\frac{2 m_c^2}{s}\right) }{s (s - q^2) } \nonumber \\= & {} 2 \left( {\mathcal {C}}_2+\frac{1}{3} {\mathcal {C}}_1\right) \left[ h_S(m_c,q^2)\right. \nonumber \\&\left. -\frac{1}{3}h_P(m_c,q^2) \right] \,. \end{aligned}$$

(22)

This expression is certainly not a good approximation of $\Delta Y_{c {\bar{c}}} (q^2)$ close to the resonance region; however, it is expected to provide a reasonable approximation at $q^2 \approx 0$, up to $O(\Lambda _\mathrm{QCD}/m_c^2)$ corrections. We can thus use it to set bounds on the slope of $\Delta Y_{c {\bar{c}}} (q^2)$ in the vicinity of $q^2=0$. The perturbative result implies

$$\begin{aligned} \left. \frac{d}{ dq^2} \Delta Y^\mathrm{pert}_{c {\bar{c}}} (q^2) \right| _{q^2=0}= & {} \frac{4}{15}\left( {\mathcal {C}}_2 + \frac{1}{3}{\mathcal {C}}_1 \right) \frac{1}{m_c^2} \nonumber \\\approx & {} (1.7\pm 1.7) \times 10^{-2} ~\mathrm{GeV}^{-2}~, \end{aligned}$$

(23)

where the numerical value has been obtained setting $m_b/2< \mu < 2m_b$ and $m_c=1.3$ GeV. According to the analysis of Ref. [26], the inclusion of $O(\Lambda _\mathrm{QCD}/m_c^2, \alpha _s)$ corrections (which involve new hadronic matrix elements) modifies the above prediction to $-(0.5 \pm 0.2) \times 10^{-2}~\mathrm{GeV}^{-2}$. Given these considerations, in the numerical analysis we employ the following constraints

$$\begin{aligned}&\mathrm{Re}\left[ \sum _{j=\Psi (1S),\ldots } \eta _{j} e^{i \delta _j} \frac{\Gamma _j}{m^3_j} + \eta _{{\bar{D}}} e^{i \delta _j} \frac{1}{6 m_{{\bar{D}}}^2} \right. \nonumber \\&\quad \left. + \sum _{j=D, D^*} \eta _{j} e^{i \delta _j} \frac{1}{10 m_j^2} \right] = (1.7\pm 2.2) \times 10^{-2}~\mathrm{GeV}^{-2}~, \nonumber \\&\quad \left| \sum _{j=\Psi (1S),\ldots } \eta _{j} e^{i \delta _j} \frac{\Gamma _j}{m^3_j} + \eta _{{\bar{D}}} e^{i \delta _j} \frac{1}{6 m_{{\bar{D}}}^2} + \sum _{j=D, D^*} \eta _{j} e^{i \delta _j} \frac{1}{10 m_j^2} \right| \nonumber \\&\qquad \le 5 \times 10^{-2}~\mathrm{GeV}^{-2}~, \end{aligned}$$

(24)

where we slightly enlarged the error from (23), such that the $1\sigma$ range covers the difference between the central value in (23) and the one including $O(\Lambda _\mathrm{QCD}/m_c^2, \alpha _s)$ corrections estimated in Ref. [26].

II. Upper bound on the $|\eta _j|$ in $\Delta Y^\mathrm{2P}_{c {\bar{c}}} (q^2)$.

The comparison of the perturbative result with $\Delta Y^\mathrm{2P}_{c {\bar{c}}} (q^2)$ also allows us to define the natural range for the $\eta _{{\bar{D}}, D, D^*}$ parameters, which are poorly constrained by data. Focusing the attention on the leading S-wave contribution, it turns out that the perturbative quark loop can be saturated by the $DD^*$ meson loop, in the limit $m_c \rightarrow m_{{\bar{D}}}$, setting $\eta _{{\bar{D}}} = 2 ({{\mathcal {C}}}_{2} + {{\mathcal {C}}}_{1/3} ) \approx (0.2 \pm 0.2).$ On general grounds, each of the exclusive meson contributions should be significantly smaller than the inclusive quark contribution. As a result, in the following we set the upper limit

$$\begin{aligned} \left| \eta _{{\bar{D}}, D, D^*} \right| \le 0.2~. \end{aligned}$$

(25)

III. Upper bound on $|Y^\mathrm{1P}_\mathrm{light} (q^2=0)|$.

Using an unsubtracted dispersion relation and taking into account only one-particle intermediate states for the light-quark contributions implies $Y^\mathrm{1P}_\mathrm{light} (q^2) \rightarrow 0$ for large $q^2$, while $Y^\mathrm{1P}_\mathrm{light} (0) \not = 0$. More precisely, one finds a power-like suppression of the type $Y^\mathrm{1P}_\mathrm{light} (q^2) \sim \Lambda _{\mathrm {QCD}}^2/q^2$ at large $q^2$, whereas $Y^\mathrm{1P}_\mathrm{light} (0)$ is not parametrically suppressed by any scale ratio. However, since the Wilson coefficients entering $Y^\mathrm{1P}_\mathrm{light}$ are strongly suppressed, either by loop factors or by subleading CKM factors, $|Y^\mathrm{1P}_\mathrm{light} (0)|$ cannot be too large. Parametrically we expect

$$\begin{aligned} | Y^\mathrm{1P}_\mathrm{light} (0) | < O(1) \times \mathrm{max} \{ |{\mathcal {C}}_{3\ldots 6}|, | {\mathcal {C}}_{1,2}^{u} | \}~. \end{aligned}$$

(26)

Taking the size of the ${\mathcal {C}}_i$ into account, we set the conservative bound¹

$$\begin{aligned} \left| Y^\mathrm{1P}_\mathrm{light} (0) \right| \approx \left| \sum _{j = \rho , \omega , \phi } \eta _j \frac{\Gamma _j}{m_j} \right| \le 0.1~, \end{aligned}$$

(27)

which should be interpreted as a constraint on the relative phases of the light resonances.

2.3.5 Estimate of the “model error”

Despite not being entirely dictated by first principles, the parameterisation of long-distance hadronic contributions discussed so far contains all the relevant one- and two-particle discontinuities of the amplitude, with free coefficients to be fixed by data. It should therefore provide a sufficiently general (and unbiased) description of the impact of hadronic contributions on the $B^+\rightarrow K^+ \mu ^+\mu ^-$ spectrum. Still, it may be worthwhile to assess whether the proposed parameterisation influences the extraction of information on $Y_{\tau \bar{\tau }} (q^2)$ and, correspondingly, the extraction of a bound on ${\mathcal {B}}(B^+\rightarrow K^+ \tau ^+\tau ^-)$. An estimate of this “model error” can be obtained by examining the stability of the obtained bound on ${\mathcal {B}}(B^+\rightarrow K^+ \tau ^+\tau ^-)$ under small variations of the model assumptions. The latter include: (i) the use of subtracted vs. unsubtracted dispersion relations for the light resonances; (ii) the use of $q^2$-dependent widths for both charmonium and/or light resonances; (iii) strengthening or relaxing the theoretical constraints in Eqs. (24), (25), and (27).

2.4 Tau-lepton contribution

The contribution from the intermediate $\tau$-leptons can be computed in perturbation theory, yielding

$$\begin{aligned} Y_{\tau {\bar{\tau }}}(q^2)=-\frac{\alpha _\mathrm {em}}{2\pi }\,{\mathcal {C}}_9^\tau \, \left[ h_S(m_\tau ,q^2)-\frac{1}{3}h_P(m_\tau ,q^2) \right] \,, \end{aligned}$$

(28)

with the functions $h_L(m,s)$ defined in Eq. (19). The functional form is identical to the one of the perturbative charm contribution and, to a large extent, to the one of the $DD^*$ contribution, illustrated in Fig. 2. However, the cusp is located at $q^2=4m_\tau ^2$, sufficiently well separated from the various hadronic thresholds.

In principle, the short-distance $b\rightarrow s\tau ^+\tau ^-$ amplitude does not need to be controlled by the CKM matrix in a generic NP model. However, in most realistic scenarios the weak phases of all $b\rightarrow s l ^+ l ^-$ amplitudes are aligned to the SM one, implying $\mathrm{Im}({\mathcal {C}}_9^\tau )= \mathrm{Im}({\mathcal {C}}_9^\mu )=0$. In the following, we adopt this (motivated) simplifying assumption.

An estimate of the maximal allowed size of $|{\mathcal {C}}_9^\tau |$ can be derived from the experimental upper bound on ${\mathcal {B}}(B^+ \rightarrow K^+ \tau ^+\tau ^-)< 2.25\times 10^{-3}$ at 90% CL by Babar [19], which is more than four orders of magnitude larger than ${\mathcal {B}}(B^+ \rightarrow K^+ \tau ^+\tau ^-)_\mathrm {SM}\approx 1.5 \times 10^{-7}$ [32]. Neglecting the contributions from operators other than ${\mathcal {O}}^\tau _9$ and ${\mathcal {O}}^\tau _{10}$, we find

$$\begin{aligned} {\mathcal {B}}(B^+ \rightarrow K^+ \tau ^{+}\tau ^-) {\approx } \left\{ \begin{array}{ll} 8.7 {\times } 10^{-9} {\times } | {\mathcal {C}}^\tau _{9} |^2; \qquad &{} C_{10}^\tau = - C_{9}^\tau ,\\ 2.7 {\times } 10^{-9} {\times } {|} {{\mathcal {C}}}^\tau _{9} |^2; \qquad &{} {\mathcal {C}}^\tau _{10}=0~. \\ \end{array} \right. \end{aligned}$$

(29)

In the case ${\mathcal {C}}^\tau _9={\mathcal {C}}^\tau _{10}$ (${\mathcal {C}}^\tau _{10}=0$) the Babar result then implies $|{\mathcal {C}}^\tau _{9} | \le 5.1 \times 10^2~ (9.1 \times 10^2)$, to be compared to ${\mathcal {C}}_9^{\tau ,\mathrm {SM}}\approx 4.2$. As we discuss below, saturating this bound leads to a pronounced ditau cusp in the spectrum (see Fig. 3), opening the possibility of extracting a more stringent bound on ${\mathcal {B}}(B^+ \rightarrow K^+ \tau ^+\tau ^-)$ from a precise measurement of the $B^+\rightarrow K^+ \mu ^+\mu ^-$ dimuon spectrum.

3 Analysis of the expected sensitivity at LHCb

In order to assess the sensitivity to the branching ratio ${\mathcal {B}}(B^+ \rightarrow K^+ \tau ^+\tau ^-)$ at the LHCb experiment, we generate pseudo-experiments corresponding to the signal yields obtained in Ref. [30] and scaled to the full run II dataset, taking into account the collected luminosity and b-hadron cross-section increase at 13 TeV [33]. This leads to around 40,000 $B^{+}\rightarrow K^{+} \mu ^{+}\mu ^{-}$ candidates (cutting the two narrow resonances). As the efficiency is reasonably flat as a function of dimuon mass and the background level is very low, we neglect these effects. Figure 3 shows the fit model with a dataset generated at the expected yield. This illustrates the visible sensitivity to a hypothetical signal component generated according to the current experimental limit [19].

Fig. 3

Example pseudodata expected from the full run II dataset collected by the LHCb experiment assuming the SM. The distribution expected if the $B^+\rightarrow K^+ \tau ^+\tau ^-$ branching fraction were present at the current experimental limit of $2.25\times 10^{-3}$ is overlaid

Table 1 Sensitivity to $C_9^{\tau }$ according to various LHCb scenarios [${\mathcal {C}}_9^{\tau ,\mathrm {SM}}\approx 4.2$]

The size and phase of the one-particle resonant contributions are determined from the branching fractions reported in Ref. [33], which are used to determine the initial values of $\eta _j$ and $\delta _j$ for the data to be generated. Due to the complicated experimental resolution effects near the $J/\psi$ and $\psi (2S)$ resonances, the regions $9.2< q^{2} < 10.0\,\mathrm{GeV}^{2}/c^{4}$ and $13.2< q^{2} < 13.95 \,\mathrm{GeV}^{2}/c^{4}$ are excluded from the fit and the phase differences associated with these resonances are constrained to the uncertainties in Ref. [30]. Outside of this region, finite-resolution effects in $q^2$ are ignored as all the components are broad. In order to mimic the sensitivity one would have when fitting the data, Gaussian constraints are applied to the $J/\psi$ and $\psi (2S)$ resonant parameters according to the uncertainties reported in Ref. [33].

For the two-particle hadronic contribution, we conservatively allow the magnitude and phase of them to vary in the fit. As the shape of the $\hat{\rho }_{DD}$ and $\hat{\rho }_{D*D*}$ spectral densities are very similar, we combine them with an equal contribution to avoid large correlations in the fit.

The form factor uncertainties are taken from Ref. [22] and are implemented in the fit as a multivariate Gaussian constraint. The data slightly helps constrain the form factor parameters, but this affects the sensitivity on ${\mathcal {C}}_9^\tau$ only in a mild way.

The expected sensitivity on the ${\mathcal {C}}_9^\tau$ contribution is determined using the CLs method [34]. The sensitivity with the current dataset is reported in Table 1, along with two other potential future scenarios corresponding to the LHCb upgrade-II luminosity and a hypothetical improvement of the form factor uncertainties by a factor of three. The estimated sensitivity utilising the run I–II datset corresponds to a limit on the $B^+\rightarrow K^+\tau ^+\tau ^-$ branching ratio which is slightly more stringent than the current constraints placed by the BaBar collaboration and is expected to compete with the projected sensitivity of the Belle-II experiment when more data is collected.

3.1 Interplay with the hadronic contributions

The $B^+\rightarrow K^+ \tau ^+\tau ^- \rightarrow K \mu ^+\mu ^-$ re-scattering leads to two main features in the $B^+\rightarrow K^+\mu ^+\mu ^-$ dimuon spectrum: (i) the cusp in-between the $J/\psi$ and $\psi (2S)$ resonances, and (ii) a distortion in the shape of the spectrum before the two resonances. The effect after the $\psi (2S)$ peak is less relevant since in that region the spectrum is rather discontinuous due to the various one- and two-particle thresholds. In order to investigate the sensitivity to the cusp feature, we have performed a fit limited to the region between the $J/\psi$ and $\psi (2S)$ resonances: this leads to a sensitivity to ${\mathcal {C}}_9^\tau$ diluted by a factor of four. We thus conclude that is the deformation of the spectrum, in particular before the $J/\psi$, that generates the largest sensitivity to ${\mathcal {C}}_9^\tau$. This implies that neglecting the resolution is justified.

Since the deformation of the spectrum at low $q^2$ plays a relevant role, we deduce that the assumed shape of the charmonium contribution is an important ingredient in constraining the $B^+\rightarrow K^+ \tau ^+\tau ^-$ signal. The component which most closely resembles the signal is the contribution from two-particle hadronic intermediate states. This is reflected in a correlation coefficient between this amplitude and the signal of about 0.6. However, the two are clearly distinct given the different location of the thresholds. We also explicitly checked that the theoretical constraints described in Eqs. (24) and (25) do not affect our sensitivity estimate: the best fit value of the two-particle hadronic contribution lies far from these bounds. This leaves open the possibility of a further increase of sensitivity with more stringent constraints on the two-particle hadronic contribution, which could be derived using $B \rightarrow D D^* K$ data. We finally note that the correlation between the hadronic contribution and the signal is fully taken into account in the sensitivity estimates using the CLs method.

4 Conclusions

If the branching ratio $\mathcal {B}(B^+ \rightarrow K^+\tau ^+\tau ^-)$ were significantly enhanced over its SM value, it would induce a peculiar distortion of the $B^+\rightarrow K^+ \mu ^+\mu ^-$ spectrum, characterised by a cusp at $q^2=4m_\tau ^2$ and by a distortion of the dimuon distribution. In this work we have proposed a method that uses this effect as a tool to extract a bound on ${\mathcal {B}}(B^+ \rightarrow K^+\tau ^+\tau ^-)$ from future precise measurements of $\mathrm{d} \Gamma (B^+\rightarrow K^+ \mu ^+\mu ^-)/ \mathrm{d} q^2$.

A necessary ingredient to achieve this goal is a reliable description of the $B^+\rightarrow K^+ \mu ^+\mu ^-$ dimuon spectrum, within the SM, in the full kinematical range, especially in the region before and within the narrow charmonium states. As we have shown, this can be obtained by means of a data-driven approach which takes full advantage of the known analytic properties of the decay amplitude, supplemented by robust theoretical constraints. Our approach differs from previous attempts of including non-local hadronic contributions to the $B^+\rightarrow K^+ \mu ^+\mu ^-$ decay amplitude by three main points: (i) the use of dispersion relations subtracted at $q^2=0$ for the charmonium states; (ii) the inclusion of two-particle thresholds; (iii) the use of short-distance constraints at low $q^2$ to reduce the number of free parameters. In this way one separates the problem of the normalisation of the $B^+\rightarrow K^+ \mu ^+\mu ^-$ rate, and the corresponding extraction of short-distance Wilson coefficients, from the problem of obtaining a reliable description of the dimuon spectrum. While within our approach there is no significant progress on the first problem, there is a tangible advantage on the second one. The parameterisation of the amplitude we propose is flexible enough to allow the extraction of all the relevant parameters characterising hadronic thresholds in the dimuon spectrum from data, while retaining significant predictive power in the smooth region within and below the two narrow charmonium resonances. This fact is the key property which allows us to set useful constraints on the $B^+\rightarrow K^+ \tau ^+\tau ^- \rightarrow K^+ \mu ^+\mu ^-$ re-scattering from future precise measurements of $\mathrm{d} \Gamma (B^+\rightarrow K^+ \mu ^+\mu ^-)/ \mathrm{d} q^2$.

The method we have proposed is particularly well suited for the LHCb experiment, which has already collected a large sample of $B^+\rightarrow K^+ \mu ^+\mu ^-$ events and has an excellent resolution in the dimuon spectrum [30]. As we have shown, the data already collected in run II should allow to set a bound on ${\mathcal {B}}(B \rightarrow K\tau ^+\tau ^-)$ of $O(10^{-3})$, competitive with current direct bounds (see Table 1). Bounds of $O(10^{-4})$ could be obtained with the LHCb upgrade-II luminosity.