Enhanced CP asymmetries in $B \rightarrow K \mu ^+ \mu ^-$

Abstract

We show that the current values of $R_K^\mathrm {exp}$ and $R_{K^*}^\mathrm {exp}$ can be accommodated by allowing a nonzero New Physics coupling $\delta C_9^{\mu \mu }$ to be complex, both in the scenario in which only $\delta C_9^{\mu \mu }$ is affected, and in the scenario with complex $\delta C_{9,10}^{\mu \mu }$ satisfying $\delta C_{9}^{\mu \mu }=-\delta C_{10}^{\mu \mu }$. A presence of the weak CP-violating phase can then be tested by measuring the CP-asymmetry, $\mathcal {A}_\mathrm {CP}$. We show that this asymmetry is enhanced around the peak of each $c\bar{c}$-resonance, and in fact more pronounced in the close vicinity of $J/\psi$ and $\psi (2S)$. Therefore, measuring $\mathcal {A}_\mathrm {CP}$ before and after the resonances’ peak could be revelatory of the CP-violation that originates from beyond the Standard Model, or to be a significant constrain when building a realistic scenario of New Physics.

1 Introduction

One of the most intriguing observations made at LHCb so far is the indication of the lepton flavor universality violation (LFUV). Regarding the $b\rightarrow s\ell \ell$ processes it was found that the measured $R_{K^{(*)}} = {{\mathcal B}' (B\rightarrow K^{(*)} \mu \mu )}/{{\mathcal B}' ( B\rightarrow K^{(*)} e e)}$, where ${\mathcal B}'$ stands for the partial branching fraction, is smaller than predicted in the Standard Model (SM). For example, in the bin of $q^2 \in [1.1,6]\mathrm {\,GeV}^2$ the measured values, $R_K = 0.846^{+0.062}_{-0.056}$ [1, 2] and $R_{K^*} = 0.69^{+0.12}_{-0.09}$ [3], are both about $2.5\,\sigma$ smaller than their SM estimate, $R_{K^{(*)}}=1.00(1)$ [4].

To understand the origin of such a discrepancy between theory and experiment one readily extends the effective field theory beyond the SM and from a fit to the experimental data extracted from the full angular distribution of $B\rightarrow K^{(*)} \mu \mu$ [5], one can deduce which scenario of New Physics is preferred [6,7,8]. On the basis of such analyses and the available experimental data it is reasonable to assume that the LFUV in $R_{K^{(*)}}^\mathrm {exp}<R_{K^{(*)}}^\mathrm {SM}$ is due to a more pronounced coupling of New Physics to the muon pair in the final state. Being much smaller, such a coupling to the electron pair can be assumed to be zero. Of all the possibilities tested through the global analyses of the exclusive $b\rightarrow s\ell \ell$ modes the most favored scenarios are those in which the signs of physics beyond the Standard Model (BSM) arise from coupling to the vector muonic current, i.e. $\delta C_9\ne 0$, and the one in which $\delta C_9=-\delta C_{10}\ne 0$, where $\delta C_i\equiv \delta C_i^{\mu \mu } = C_i^{\mu \mu } - C_i^\mathrm {SM}$ stands for the contribution arising solely from BSM. In order to keep the number of free parameters minimalistic a common assumption is that $\delta C_i \in \mathbb {R}$. Starting from results based on the effective field theory approach one can then build a specific model verifying either $\delta C_9 \ne 0$ or $\delta C_9 =-\delta C_{10} \ne 0$.

Of course, the assumption $\delta C_i \in \mathbb {R}$ helps keeping the number of BSM (real) parameters minimalistic but that assumption should be scrutinized experimentally. A physical consequence of allowing $\delta C_i \in \mathbb {C}$ is that the BSM effects of CP violation can be visible from the di-muon invariant mass spectra of the decays we consider here . A presence of the imaginary part in the BSM contribution to the Wilson coefficient(s) could be tested through various CP-asymmetries which are known to be tiny in the SM [9,10,11,12,13,14]. LHCb already attempted measuring direct CP asymmetry

$$\begin{aligned} \mathcal {A}_\mathrm {CP}^{K^{(*)}} = {{\mathcal B}(\overline{B} \rightarrow \overline{K}^{(*)} \mu \mu )- {\mathcal B}( B \rightarrow K^{(*)} \mu \mu ) \over {\mathcal B}(\overline{B} \rightarrow \overline{K}^{(*)} \mu \mu )+ {\mathcal B}( B \rightarrow K^{(*)} \mu \mu )} \,, \end{aligned}$$

(1)

and reported [15]¹:

$$\begin{aligned} \mathcal {A}_\mathrm {CP}^{K^+} =0.012(17)(1),\quad \mathcal {A}_\mathrm {CP}^{K^{0 *}} =-0.035(24)\,, \end{aligned}$$

(2)

which suggest that the imaginary part of the Wilson coefficient(s) is likely to be small and is currently consistent with zero. Measuring this quantity is complicated either because it is small, or because it varies within a bin so that its integrated value over the size of one bin is very small. Furthermore, the complications arising from the overwhelming presence of the $c\bar{c}$-resonances in the $\mu \mu$-spectrum makes this measurement even more challenging. A common practice in experiment (and also in theory) is not to measure around the narrow resonances, $q^2\approx m_{J/\psi ,\psi (2S)}^2$, while in the region of large $q^2$’s – where the $c\bar{c}$-resonances are broader – the spectrum is measured but for the comparison with theory it is important to integrate over a sufficiently large bin in order to rely on the quark-hadron duality. Contrary to that practice, as we argue in this work, it turns out that measuring $\mathcal {A}_\mathrm {CP}(q^2)$ around the resonance region can be more beneficial than measuring it away from resonances.

In the LHCb analysis [16] the CP-averaged $B \rightarrow K \mu ^+ \mu ^-$ spectrum was carefully measured including the regions close to the $c\bar{c}$-resonances. The spectrum was modeled by amplitude which contained Breit–Wigner parameterization of each of the resonances and the fit results revealed the strength of the resonances and the corresponding strong phases, which are a crucial ingredient to predicting the direct CP asymmetry. We use the LHCb results to interpret existing measurements of $\mathcal {A}_\mathrm {CP}$ in the scope of BSM model with V-A interaction ($\delta C_9 = -\delta C_{10}$), and in the vector current BSM model ($\delta C_9$). In both cases the effective BSM couplings are allowed to be complex in order to entail a non-zero $\mathcal {A}_\mathrm {CP}$. The goal of this paper is to present the constraints on the complex BSM couplings from existing measurements and to show that, based on the current knowledge of resonant strong phases, $\mathcal {A}_\mathrm {CP}(q^2)$ is enhanced in the regions close to resonances.

2 Constraints on complex $\delta C_9$

In this section we present how the measurements of theoretically and experimentally well understood quantities, $R_{K^{(*)}}$, $\mathcal {B}(B_s \rightarrow \mu ^+ \mu ^-)$, and CP-averaged $d\Gamma /dq^2(B \rightarrow K \mu ^+ \mu ^-)$, translate into allowed regions in the complex $\delta C_9$ plane. The mentioned constraints are CP-even and are only sensitive to $(\mathrm {Im}\,\delta C_9)^2$. Furthermore we will include the constraint stemming from the measured $\mathcal {A}_\mathrm {CP}$ of $B \rightarrow K \mu ^+ \mu ^-$, which depends linearly on $\mathrm {Im}(\delta C_9)$. Constraints on $\delta C_9$ will be presented for a BSM scenario that admits $\delta C_9 = -\delta C_{10}$, as well as for the scenario where BSM is present only in $\delta C_9$.

2.1 $B\rightarrow K \ell ^+ \ell ^-$ effective Hamiltonian

Here we remind the reader of the basic ingredients needed to compute the differential decay rate of $B\rightarrow K \ell ^+ \ell ^-$, where $\ell$ stands for one of the lepton flavors. As usual, the starting point is the effective Hamiltonian describing the $b \rightarrow s \ell ^+ \ell ^-$ transitions, namely

$$\begin{aligned} \mathcal {H}^{b\rightarrow s \ell \ell }_{\mathrm {eff}} = -\dfrac{4 G_F V_{tb}V_{ts}^*}{\sqrt{2}} \sum _{i=7,9,10} C_i(\mu )\mathcal {O}_i(\mu )\,, \end{aligned}$$

(3)

where the short distance physics is encoded in the Wilson coefficients $C_{7,9,10}$, while the long distance part is described by the hadronic matrix elements of the effective operators

$$\begin{aligned} \mathcal {O}_7&= \dfrac{e m_b}{4\pi } \, (\bar{s}_R \sigma _{\mu \nu } b_R ) F^{\mu \nu }\,, \end{aligned}$$

(4)

$$\begin{aligned} \mathcal {O}_{9(10)}&= \dfrac{e^2}{(4\pi )^2} \,(\bar{s}_L\gamma _\mu b_L ) (\bar{\ell }\gamma ^\mu (\gamma ^5)\ell )\,. \end{aligned}$$

(5)

Notice that we focus here only to the operators the Wilson coefficients of which are non-zero in the Standard Model. While discussing the effects of New Physics relevant to $B\rightarrow K\mu ^+\mu ^-$ we will assume them to either modify only $C_9$ as $C_9=C_9^\mathrm {SM}+\delta C_9$, or to modify both $C_9$ and $C_{10}$ in such a way that $\delta C_9 = - \delta C_{10}$. The differential decay rate of $\bar{B}(p)\rightarrow \bar{K}(k)\ell ^+ \ell ^-$ can then be compactly written as [9]

$$\begin{aligned} {d\bar{\Gamma }\over dq^2}&= 2 \mathcal {N}(q^2)\,\biggl [ \frac{1}{6}\left( 1+ \frac{2m_\ell ^2}{q^2}\right) \lambda (q^2) \left( \vert F_V\vert ^2 + \vert F_A\vert ^2 \right) \nonumber \\&\quad + 4m_\ell ^2 m_B^2 \vert F_A\vert ^2 - q^2 \vert F_P\vert ^2 \nonumber \\&\quad + 2 m_\ell (m_B^2 -m_K^2-q^2) \mathrm {Re}\left( F_P\,F_A^*\right) \biggr ] \,, \end{aligned}$$

(6)

where $q^2=(p-k)^2$, $\lambda (q^2) = [q^2 - (m_B-m_K)^2] [q^2 - (m_B+m_K)^2]$, while the explicit expressions of the $q^2$-dependent functions $F_{V,A,P}$ read:

$$\begin{aligned} F_V&= \ C_9\, f_+(q^2) + \frac{2 m_b}{m_B+m_K} C_7\, f_T(q^2)\,, \nonumber \\ F_A&= \ C_{10}\, f_+(q^2) \,,\nonumber \\ F_P&= \ C_{10}\, m_\ell \left[ f_+(q^2) - \frac{m_B^2 - m_K^2}{q^2} \left( f_0(q^2) - f_+(q^2) \right) \right] \,. \end{aligned}$$

(7)

The normalization is also $q^2$-dependent:

$$\begin{aligned} \mathcal {N}(q^2)&= \frac{G_F^2 \alpha ^2 |V_{tb} V_{ts}^{*} |^2}{512 \pi ^5 m_B^3} \, \sqrt{\lambda (q^2) }\, \sqrt{1-\frac{4 m_\ell ^2}{q^2}} \,. \end{aligned}$$

(8)

In the above expressions, besides the Wilson coefficients $C_{7,9,10}$, we used the hadronic form factors $f_{+,0,T}(q^2)$ which parametrize the hadronic matrix elements as follows:

$$\begin{aligned} \langle K(k)| \bar{s} \gamma _{\mu } b |B(p)\rangle&= \left[ (p + k)_\mu - {m_B^2 - m_K^2 \over q^2} q_\mu \right] f_{+}(q^2) \, \nonumber \\&\quad + {m_B^2-m_K^2 \over q^2} q_{\mu } f_{0}(q^2) \, , \nonumber \\ \langle K(k)| \bar{s}\sigma _{\mu \nu }b |B(p) \rangle&= -i \left( p_\mu k_\nu - p_\nu k_\mu \right) \frac{2 f_T(q^2)}{m_B + m_K} \, . \end{aligned}$$

(9)

The above form factors have been computed by means of numerical simulations of QCD on the lattice [17,18,19] and we will use those results in our phenomenological discussion.

2.2 Resonant $B \rightarrow K \ell ^+ \ell ^-$ spectrum

A non-resonant contribution of the $c\bar{c}$-pairs, as well as those arising from the light quarks, is usually included by promoting the Wilson coefficient $C_9$ to $C_9(q^2)= C_9+ Y(q^2)$ where, owing to the quark-hadron duality, the function $Y(q^2)$ is computed perturbatively. That obviously cannot account for the $c\bar{c}$-resonances, present in the $q^2$-spectrum of the decay. For that reason, when comparing theory with experiment, one vetoes regions around the prominent resonances, such as $J/\psi$, $\psi (2S)$ and $\psi (3770)$, or by working in the low $q^2$ region ($q^2 \lesssim m_{J/\psi }^2$) in order to avoid the $c\bar{c}$-resonances altogether. Since the purpose of this work is to discuss the potential effects of CP-violation, we will not follow the usual description but, instead, we will adopt the model of Ref. [16] in which the authors actually reconstructed contributions from the $c \bar{c}$-resonances through the fit to the experimental data. ² More specifically they trade $C_9$ for

$$\begin{aligned} C_9^\mathrm {eff} (q^2)&= C_9 + C_9^\mathrm {res}(q^2) \nonumber \\&= C_9 + \sum _{j} \frac{m_j \,\Gamma _j \,\eta _j \, e^{i \delta _j}}{m_j^2-q^2 - i\, m_j \Gamma _j(q^2)}, \end{aligned}$$

(10)

with $j \in \{J/\psi , \psi (2S),\psi (3770),\psi (4040), \psi (4160), \psi (4415)\}$, for which the masses ($m_j$) and widths ($\Gamma _j$) are well known [20]. The $q^2$-dependent width function reads $\Gamma _j(q^2) = \Gamma _j \times \sqrt{1-4m_\ell ^2/q^2}/ \sqrt{1-4m_\ell ^2/m_i^2}$. By using $C_9^\mathrm {eff} (q^2)$ of Eq. (10) to fit the measured CP-averaged spectrum the authors of Ref. [16] were able to determine the value of $\eta _j$ and the strong phase $\delta _j$ for each of the six $c\bar{c}$-resonances. The Wilson coefficients $C_9$ and $C_{10}$, assumed to be real, have been also fitted in the analysis, with their values in line with results of the global fits. It appears though that there is a fourfold ambiguity related to the choice of the signs of the first two resonances. By using the results for all $\eta _j$ and $\delta _j$ from Ref. [16] we were able to reconstruct their model. Notice also that in Eq. (10), on the right hand side, we have $C_9=C_9^\mathrm {SM}+\delta C_9$. While the authors of Ref. [16] assumed $\delta C_9$ to be real, we will allow it to be complex. In other words, we allow for a possible BSM weak phase.

2.3 $R_{K^{(*)}}$ and $\mathcal {B}(B_s \rightarrow \mu ^+ \mu ^-)$

Fig. 1

Constraints on $\bigl ( \mathrm {Re}(\delta C_9), \mathrm {Im}(\delta C_9)\bigr )$ compatible with $R_K^\mathrm {exp}$ and $R_{K^*}^\mathrm {exp}$, as well as with the error bars on $\mathcal {A}_\mathrm {CP}$ measured in the bins between $2 \mathrm {\,GeV}^2< q^2 < 8\mathrm {\,GeV}^2$, as reported in Ref. [15]. The horizontal lines correspond to the $\mathcal {A}_\mathrm {CP}$ constraint for branches 1 and 2, see text for details. Vertical dashed lines enclose the $2\sigma$ region of $\mathrm {Re}(\delta C_9)$ obtained from CP-averaged spectra in Ref. [16]. The scenario considered in this case is the one in which all of the BSM effects are described by $\delta C_9$ only. The fitted $1\sigma$ values of $\delta C_9$ are ellipses in gray (red) for strong phases in branch 1 (2). By a star we denote our benchmark point which we chose to be $\delta C_9 = -0.85-0.73 i$

As we already mentioned above, we consider two scenarios of New Physics. In the first one we will allow only $C_9$ to receive an extra contribution, $\delta C_9 \in \mathbb {C}$, while leaving other Wilson coefficients at their Standard Model values. In such a situation it suffices to use $R_{K^{(*)}}$ to determine $\delta C_9$. Knowing that $R_K^\mathrm {exp}$ and $R_{K^*}^\mathrm {exp}$ are obtained from the partial branching fractions integrated in the interval $q^2\in [1.1,6]~\mathrm{GeV}^2$, we obtain the simple formulas:

$$\begin{aligned} R_K&= 1.003 + 0.244\, \mathrm {Re}(\delta C_9)\nonumber \\&\quad + 4.01\times 10^{-3}\, \mathrm {Im}(\delta C_9) + 0.028\, |\delta C_9|^2,\end{aligned}$$

(11a)

$$\begin{aligned} R_{K^*}&= 0.997 + 0.202\,\mathrm {Re}(\delta C_9)\nonumber \\&\quad + 1.65\times 10^{-3}\, \mathrm {Im}(\delta C_9) + 0.033\, |\delta C_9|^2, \end{aligned}$$

(11b)

in both of which most of the $1\%$ overall error affects the first term. Using these two expressions, together with $R_K^\mathrm {exp}$ and $R_{K^*}^\mathrm {exp}$, we get a region of allowed values for $\bigl ( \mathrm {Re}(\delta C_9), \mathrm {Im}(\delta C_9)\bigr )$ shown in Fig. 1.

Fig. 2

In addition to the constraints mentioned in the caption of Fig. 1, here we also account for $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp}$. The scenario we consider in this plot is the one where the BSM effects are described by $\delta C_9$ and $\delta C_{10}$, satisfying $\delta C_9=-\delta C_{10}$. Our benchmark point in this case, $\delta C_9 = -0.48-0.7 i$, is depicted by a star

In the second scenario considered in this work, we allow both $C_{9,10}$ to receive contributions from BSM, $C_{9,10}=C_{9,10}^\mathrm {SM}+\delta C_{9,10}$, but by respecting the left-handedness, i.e. $\delta C_9 = -\delta C_{10}$. In this case, the above formulas for $R_{K^{(*)}}$ become:

$$\begin{aligned} R_K&= 1.003 + 0.477\, \mathrm {Re}(\delta C_9)\nonumber \\&\quad + 4.01\times 10^{-3}\, \mathrm {Im}(\delta C_9) + 0.057 |\delta C_9|^2, \end{aligned}$$

(12a)

$$\begin{aligned} R_{K^*}&= 0.997 + 0.472 \,\mathrm {Re}(\delta C_9)\nonumber \\&\quad + 1.65\times 10^{-3}\, \mathrm {Im}(\delta C_9) + 0.066 |\delta C_9|^2. \end{aligned}$$

(12b)

Since we allow in this scenario $\delta C_{10}\ne 0$, the recently updated $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp}=(2.69^{+0.37}_{-0.35}) \times 10^{-9}$ [21] becomes an important constraint too. Notice that $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp}$ is $2.5\,\sigma$ smaller than predicted in the SM, $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {SM}=(3.66\pm 0.11) \times 10^{-9}$ [22]. For the purpose of determining $\delta C_{10} \in \mathbb {C}$, we use

$$\begin{aligned} \mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {th}&= \tau _{B_s}\dfrac{\alpha ^2 G_F^2 m_{B_s}}{16 \pi ^3} \left| V_{tb}V_{ts}^*\right| ^2 m_\mu ^2 \nonumber \\&\quad \times \sqrt{1-\frac{4 m_\mu ^2}{m_{B_s}^2}} \left| C_{10}\right| ^2 f_{B_s}^2, \end{aligned}$$

(13)

where, thanks to the lattice QCD efforts, the uncertainty in the decay constant is not anymore an obstacle to constraining the BSM contribution, $f_{B_s}= (230.3\pm 1.3)\,\mathrm{MeV}$ [17]. In order to compare Eq. (13) to $\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp}$, one also needs to account for the effect of $B_s- \overline{B}_s$ oscillations which, to a good approximation, amounts to [23]

$$\begin{aligned} \mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {exp} \approx \dfrac{1}{1-y_s}\mathcal {B}(B_s\rightarrow \mu ^+\mu ^-)^\mathrm {th}, \end{aligned}$$

(14)

where $y_s=\Delta \Gamma _{B_s}/(2 \Gamma _{B_s})=0.061(7)$ [24]. Finally, a comparison between theory and experiment in the scenario with $\delta C_9 = -\delta C_{10}$ results in a region of allowed values in the $\bigl ( \mathrm {Re}(\delta C_9), \mathrm {Im}(\delta C_9)\bigr )$ plane, which we plot in Fig. 2. We should also note that a possibility of the complex Wilson coefficients was recently discussed in Refs. [25, 26].

2.4 $\mathcal {A}_\mathrm {CP}$ constraint

The measurement of $\mathcal {A}_\mathrm {CP}$ in non-resonant regions of $B^\pm \rightarrow K^\pm \mu ^+ \mu ^-$ has been presented by the LHCb collaboration in [15]. We will employ their bin-by-bin results for $\mathcal {A}_\mathrm {CP}$ in the region of $2 \mathrm {\,GeV}^2< q^2 < 8\mathrm {\,GeV}^2$, with each bin-width being $1\mathrm {\,GeV}^2$. With binned data the CP-asymmetry is defined as

$$\begin{aligned} \mathcal {A}_\mathrm {CP}(\mathrm {bin})&= \frac{\bar{\Gamma }(\mathrm {bin}) - \Gamma (\mathrm {bin})}{\bar{\Gamma }(\mathrm {bin}) + \Gamma (\mathrm {bin})}\,,\nonumber \\ \mathrm {bin}&\equiv [q_1^2,q_2^2]\,, \end{aligned}$$

(15)

where $\bar{\Gamma }\, (\Gamma )$ refers to decay $B^- \rightarrow K^- \mu ^+ \mu ^-\,(B^+ \rightarrow K^+ \mu ^+ \mu ^-)$. Importantly, the theoretical prediction of $\mathcal {A}_\mathrm {CP}$ in each of the $[2,3]\mathrm {\,GeV}^2,\ldots ,[7,8]\mathrm {\,GeV}^2$ bins depends on the resonant spectrum too and therefore the results of Ref. [16] should be accounted for carefully. More specifically, as we shall see below, the theoretical prediction for $\mathcal {A}_\mathrm {CP}$ is proportional to $\mathrm {Im}(\delta C_9)\times \mathrm {Im}(C_9^\mathrm {res})$ in a given bin. With the fourfold sign ambiguity regarding the first two resonances, as reported in Ref. [16], it turns out that it suffices to consider two distinct solutions with negative strong phase $\delta _{J/\psi }$ and either negative (“Branch 1”) or positive value (“Branch 2”) of the strong phase $\delta _{\psi (2S)}$:

$$\begin{aligned} \text {Branch 1: }\quad \delta _{J/\psi }&= -1.66\,,\quad \delta _{\psi (2S)} = -1.93\,, \end{aligned}$$

(16a)

$$\begin{aligned} \text {Branch 2: }\quad \delta _{J/\psi }&= -1.50\,,\quad \delta _{\psi (2S)} = 2.08\,. \end{aligned}$$

(16b)

The two fit branches are shown in Figs. 1,2. The other two solutions would result in flipping the sign of $\mathrm {Im}(\delta C_9)$. The allowed region of $\delta C_9$, obtained from the fit to binned $\mathcal {A}_\mathrm {CP}$ in the region $2\,\mathrm{GeV}^2 \le q^2\le 8 \mathrm{GeV}^2$, results in horizontal lines in Figs. 1 and 2.³

3 Behavior of direct CP asymmetry at the resonance

The differential CP asymmetry is defined as

$$\begin{aligned} \mathcal {A}_\mathrm {CP}(q^2) = \frac{d\bar{\Gamma }/dq^2-d\Gamma /dq^2}{d\bar{\Gamma }/dq^2+d\Gamma /dq^2}, \end{aligned}$$

(17)

where, as mentioned above, $d\bar{\Gamma }/dq^2$ and $d\Gamma /dq^2$ refer to the differential decay rates of $\bar{B} \rightarrow \bar{K} \mu ^+ \mu ^-$ and $B \rightarrow K \mu ^+ \mu ^-$. In the following, for simplicity, we will neglect the muon mass in Eq. (6) and set $C_{10} = C_{10}^\mathrm {SM} + \delta C_{10}$, $C_9 = C_9^\mathrm {SM} + \delta C_9$, $C_7 = C_7^\mathrm {SM}$, to get:

$$\begin{aligned} \begin{aligned} \frac{d\bar{\Gamma }}{dq^2}&= \frac{\mathcal {N} \lambda }{3} \left[ f_+(q^2) \right] ^2 \Bigg \{ |C_{10}^\mathrm {SM} + \delta C_{10}|^2 \\&\quad + \left| C_9 + C_9^\mathrm {res}(q^2)+ \frac{2m_b}{m_B + m_K} \frac{f_T(q^2)}{f_+(q^2)} C_7^\mathrm {SM} \right| ^2 \Bigg \}, \end{aligned} \end{aligned}$$

(18)

where, again for notational simplicity, we omit the argument in $\mathcal {N}$ and $\lambda$. Note that $C_9^\mathrm {res}(q^2)$ is complex and contains a CP-even (strong) phase, cf. Eq. (10). $C_{10}^\mathrm {SM}$ and $C_7^\mathrm {SM}$ are assumed to be real so that the only potential sources of CP violation are in the imaginary parts of $\delta C_{9,10}$. The rate for CP-conjugated decay is obtained from (18) by replacing $\delta C_{9,10} \rightarrow \delta C_{9,10}^*$. The necessary ingredient for a non-zero $\mathcal {A}_\mathrm {CP}$ is interference between two terms of the amplitude which have different strong and weak phases. Such effect is possible only for $\mathrm {Im}(\delta C_9)$ that interferes with $\mathrm {Im}(C_9^\mathrm {res}(q^2))$ and drives the numerator of $\mathcal {A}_\mathrm {CP}$, namely:

$$\begin{aligned} \frac{d\bar{\Gamma }}{dq^2} - \frac{d\Gamma }{dq^2} = \frac{4\mathcal {N} \lambda }{3}\left[ f_+(q^2)\right] ^2\, \mathrm {Im}(C_9^\mathrm {res}(q^2)) \,\mathrm {Im}(\delta C_9).\nonumber \\ \end{aligned}$$

(19)

Fig. 3

Behavior of $\mathcal {A}_\mathrm {CP}\equiv \mathcal {A}_\mathrm {CP}(q^2)$ around the $J/\psi$-resonance for which $\eta _{J/\psi }=8500$, $\delta _{J/\psi } = -1.66$. Due to the large strong phase the $\mathcal {A}_\mathrm {CP}$ is antisymmetric with respect to the position of the $J/\psi$ peak, $q^2=9.58~\mathrm{GeV}^2$. Its maximal values are attained far away from the peak (full line). We also show the behavior of $\mathcal {A}_\mathrm {CP}$ if the strong phase was 2- (dashed) or 10-times (dotted) smaller. Dashed gridlines denote the positions and heights of $\mathcal {A}_\mathrm {CP}$ extrema according to Eq. (26)

The contribution of $\delta C_{10}$ is not important, as it only modifies the denominator of $\mathcal {A}_\mathrm {CP}$ and we neglect it in this discussion. To make the argument clearer, let us now assume that $\delta C_9$ is strictly imaginary and that $|\delta C_9|^2$ is negligible with respect to $|C_9^\mathrm {SM}|^2$ in the denominator of $\mathcal {A}_\mathrm {CP}(q^2)$. We are interested in the behavior of $\mathcal {A}_\mathrm {CP}\equiv \mathcal {A}_\mathrm {CP}(q^2)$ close to one of the aforementioned $\bar{c} c$-resonances, where $C_9^\mathrm {res}(q^2)$ is approximately

$$\begin{aligned} C_9^\mathrm {res}(q^2)\approx \frac{m_j \Gamma _j \eta _j e^{i\delta _j}}{m_j^2 - q^2 -i m_j \Gamma _j}\,. \end{aligned}$$

(20)

The expression for $\mathcal {A}_\mathrm {CP}$ then boils down to

$$\begin{aligned} \mathcal {A}_\mathrm {CP}&= \mathrm {Im}(\delta C_9)\, \frac{2\eta _j \left( \cos \delta _j- x \sin \delta _j\right) }{\eta _j^2 - 2 \eta _j B \left[ \sin \delta _j+ x \cos \delta _j\right] + A \left[ 1+ x^2 \right] }\,, \end{aligned}$$

(21)

where $x \equiv (q^2-m_j^2)/(m_j\Gamma _j)$ measures the distance from the resonance peak, while B and A are defined as

$$\begin{aligned} B&= C_{9}^\mathrm {SM} + \frac{2m_b}{m_B+m_K} \frac{ f_T(q^2)}{ f_+(q^2) } C_7^\mathrm {SM} \approx 3.8\,,\nonumber \\ A&= (C_{10}^\mathrm {SM})^2+B^2 \approx 31\,. \end{aligned}$$

(22)

The above two quantities are almost constant throughout the whole range of physical $q^2$’s. If the strong phase is large, $|\delta _j| \approx \pi /2$, then the imaginary part of $C_9^\mathrm {res}(q^2)$ vanishes on the resonant peak, leading to zero $\mathcal {A}_\mathrm {CP}$. In the limit of small strong phase, $\delta _j= 0$, the value of the asymmetry on the peak is

$$\begin{aligned} \bigl . \mathcal {A}_\mathrm {CP}(x \rightarrow 0) \biggr |_{\delta _j= 0} = \mathrm {Im}(\delta C_9)\, \frac{2 \eta _j}{\eta _j^2 + A}. \end{aligned}$$

(23)

The lowest lying and the most prominent $c\bar{c}$-resonances are $J/\psi$ and $\psi (2S)$. Their narrow widths and relatively large branching fractions to the $\mu ^+ \mu ^-$ channel imply large $\eta _j$ parameters [$\eta _{J/\psi } \approx 8.5\times 10^{3}$, $\eta _{\psi (2S)} \approx 1.4\times 10^{3}$]. Therefore, $\mathcal {A}_\mathrm {CP}$ is suppressed by $\eta _j$ at the resonant peak. Farther away from the peak, however, $\mathcal {A}_\mathrm {CP}$ is enhanced by the large values of $\eta _j$, as can be seen in the following limit

$$\begin{aligned} \mathcal {A}_\mathrm {CP}(|x| \rightarrow \infty ) = \mathrm {Im}(\delta C_9)\, \frac{2 \eta _j (\cos \delta _j-x \sin \delta _j)}{ A \, x^2}\,. \end{aligned}$$

(24)

In the large x regime the leading asymptotic term is proportional to $\sin \delta _j$,

$$\begin{aligned} \mathcal {A}_\mathrm {CP}&\approx \mathrm {Im}(\delta C_9)\,\frac{-2\sin \delta _j\,\eta _j}{Ax} \nonumber \\&= \mathrm {Im}(\delta C_9) \, \frac{-2\sin \delta _j\,\eta _j \Gamma _j m_j}{A(q^2-m_j^2)}\,. \end{aligned}$$

(25)

Finally, $\mathcal {A}_\mathrm {CP}$ also develops two extrema around the resonant peak. Explicit expression for their positions ($q^2_{1,2}$) can be derived from Eq. (21). We obtain:

$$\begin{aligned} q^2_{1,2}&= m_j^2 \pm m_j \Gamma _j \left( \frac{\eta _j}{\sqrt{A}} +\frac{B}{\sqrt{A} \sin \delta _j}\right) \nonumber \\&\quad \qquad + m_j \Gamma _j \cot \delta _j+ \mathcal {O}(1/\eta _j)\,, \end{aligned}$$

(26a)

$$\begin{aligned} \mathcal {A}_\mathrm {CP}\left( q^2_{1,2}\right)&= \mathrm {Im}( \delta C_9 )\,\frac{\sin \delta _j}{\sqrt{A} \pm B \cos \delta _j}+\mathcal {O}(1/\eta _j)\,. \end{aligned}$$

(26b)

To leading order in $1/\eta _j$ the positions of the extrema of $\mathcal {A}_\mathrm {CP}$ depend on the product $\eta _j \Gamma _j$. The maximum value of $\mathcal {A}_\mathrm {CP}$, instead, depends only on the strong phase, as can be seen from Eq. (26b). The conclusion that we can draw from this analysis is that the narrow resonances with large $\eta _j$ enhance $\mathcal {A}_\mathrm {CP}$ much farther away from the resonant peak than one would naïvely expect from the small resonance width. This feature is explicitly shown in Fig. 3 where $\mathcal {A}_\mathrm {CP}$ is plotted around the $J/\psi$-resonance. As we see the extrema of $\mathcal {A}_\mathrm {CP}$ are positioned at $x \approx \pm 1500$, i.e. at $q^2 - m_{J/\psi }^2 \approx \pm 1500\ m_{J/\psi } \Gamma _{J/\psi }$. A second message is that a large strong phase $\delta _j$ makes the distribution of $\mathcal {A}_\mathrm {CP}$ antisymmetric around the peak, whereas with smaller strong phase values the shape of $\mathcal {A}_\mathrm {CP}$ around the peak would become asymmetric and its size much smaller. Qualitative features of $\mathcal {A}_\mathrm {CP}$ around the $\psi (2S)$ are similar to those presented in Fig. 3 for $J/\psi$, as both these resonances have large $\delta _j$, barring a sign ambiguity, as determined in Ref. [16]. Thus, in the presence of the imaginary part of $\delta C_9$ the $\mathcal {A}_\mathrm {CP}$ would be enhanced around these two resonances. Furthermore, the shape of $\mathcal {A}_\mathrm {CP}(q^2)$ would be approximately antisymmetric with respect to each of the resonant peaks.

Fig. 4

In the left panel we plot $\mathcal {A}_\mathrm {CP}\equiv \mathcal {A}_\mathrm {CP}(q^2)$ in the full physical region, and for the benchmark point in the BSM scenario satisfying $\delta C_9 = - \delta C_{10}$. Two branches correspond to two sets of solutions for strong phases obtained from the fit to the data of Ref. [16]. Other two branches would correspond to flipping the sign of the above-depicted $\mathcal {A}_\mathrm {CP}(q^2)$ in the region $q^2\lesssim m_{\psi (2S)}^2$. In the right panel we plot the CP-averaged differential decay rate in this scenario to show that the effect of $\delta C_9\in \mathbb {C}$ is extremely difficult to disentangle from this quantity. Indeed, the (dashed) curve corresponding to the CP-average of decay widths and $\mathrm {Im}\left( \delta C_9\right) = - 0.7$ is hardly distinguishable from the full black curve in which $\mathrm {Im}\left( \delta C_9\right)$ is set to zero. Here the “branch 1” set of resonant parameters has been used

3.1 Shape and size of $\mathcal {A}_\mathrm {CP}(q^2)$ for benchmark values of $\delta C_9$

So far in this Section we introduced several simplifications to make our discussion clearer. We now use the full formulas to draw $\mathcal {A}_\mathrm {CP}(q^2)$ for the two benchmark points we have chosen in both of our scenarios, denoted by stars in Figs. 1,2. These benchmark points correspond to a large value of $|\mathrm {Im}( \delta C_9)|$. Following the line of the discussion above, besides a non-zero value of $\mathrm {Im}( \delta C_9)$, it is essential to have a good handle over $\mathrm {Im}(C_9^\mathrm {res}(q^2))$, which we were able to get thanks to the results of Ref. [16] in which the resonant parameters are given for four possible solutions (branches). In the region $q^2 \lesssim m_{\psi (2S)}^2$ only two of those would give different shapes of $\mathcal {A}_\mathrm {CP}(q^2)$, while the other two solutions would simply flip the overall sign of those $\mathcal {A}_\mathrm {CP}(q^2)$ that we already obtained from the first two branches. Regarding the first two branches the results are plotted in Fig. 4. We note immediately that due to the fact that $\mathcal {A}_\mathrm {CP}(q^2)$ essentially depends on $\mathrm {Im}( \delta C_9)$, and since in both benchmark points $\mathrm {Im}( \delta C_9)\approx -0.7$, the two scenarios for each branch practically coincide. This is why we decided to plot $\mathcal {A}_\mathrm {CP}(q^2)$ only in the scenario with $\delta C_9 = - \delta C_{10}$. We need to emphasize once again that it is essential to measure $\mathcal {A}_\mathrm {CP}$ in a bin before and in another bin after the peak. Another remark is that locally, near the first two resonances, the value of $\mathcal {A}_\mathrm {CP}(q^2)$ can be appreciable, about $\pm 15\%$ for our benchmark $\mathrm {Im}( \delta C_9)$. It should be reiterated, however, that the effects of CP-violation would be extremely difficult to disentangle from the differential CP-averaged decay width, as we show also in Fig. 4.

4 Conclusion

In most of the phenomenological studies of the exclusive $B\rightarrow K^{(*)}\ell ^+\ell ^-$ decays, the potential effects of New Physics are described by the shift of some of the Wilson coefficients. From the experimental data on $B\rightarrow K^{(*)}\mu ^+\mu ^-$ the most favored scenarios seem to be those modifying $C_9$, or those in which both $C_9$ and $C_{10}$ are modified but in such a way that $\delta C_9 = - \delta C_{10}$, where $\delta C_{i}$ refers to the BSM contribution. So far this shift was considered to be real-valued in order to make the number of parameters minimal. However, the BSM contributions can be complex, $\delta C_{i} \in \mathbb {C}$, and the measurement of the CP-asymmetry could reveal the presence of that BSM (weak) phase. From the current data regarding the B-anomalies, namely $R_{K^{(*)}}^\mathrm {exp} < R_{K^{(*)}}^\mathrm {SM}$, as well as from the measured $\mathcal {B}( B_s\rightarrow \mu ^+\mu ^-)$, we were able to constrain the regions in the $\bigl (\mathrm {Re}(\delta C_9), \mathrm {Im}(\delta C_9)\bigr )$ plane, clearly showing that having $\mathrm {Im}(\delta C_9)\ne 0$ is, in both scenarios, perfectly plausible and consistent with data. That effect, however, is small and its measurement along the lines presented in Ref. [15] would require extremely high experimental precision. In this work we showed that, if measured closely to the peaks of the $c\bar{c}$-resonances, the effects of $\mathrm {Im}(\delta C_9)\ne 0$ can be further amplified by the strong phase of each resonance and the resulting CP-asymmetry, $\mathcal {A}_\mathrm {CP}(q^2)$, measured before and/or after the resonance’s peak can be easier to distinguish. We showed, through a simplified example, the details of how this enhancement actually occurs. We focused on the $B\rightarrow K \mu ^+\mu ^-$ decay, but the discussion can straightforwardly be extended to $B\rightarrow K^*\mu ^+\mu ^-$. ⁴