Searching new physics in rare $B$-meson decays into multiple muons

New heavy vector bosons and light scalars are predicted in a plethora of models of new physics. In particular, in new strongly interacting sectors they play the role of the $\rho$ and $\pi$ mesons in QCD. We show that some of their interactions, for example those required for the explanation of the $B$ anomalies and the $g-2$ of the muon, can be only probed in $B$ meson decays. We highlight new golden channels not yet studied experimentally, including $B^+ \rightarrow K^+(D^+)\mu^+\mu^-\mu^+\mu^-$ and $B^0\rightarrow K^{*0}\mu^+\mu^-\mu^+\mu^-$. Relying on generator level simulations for data taking with the LHCb detector, we determine the reach of this facility to the aforementioned processes. We show that branching ratios as small as $9\times 10^{-12}$ ($3.2\times 10^{-10}$) and $2.7\times 10^{-11}$ can be tested at the $95\, \%$ CL respectively.

anomalies. However, our results are of much broader applicability. The paper is organized as follows. In Section 2 we introduce the generic Lagrangian we are interested in; we comment on constraints on the different parameters and compute the amplitudes for the different B decays. In Section 2.1, we match a particular composite Higgs model to the Lagrangian above. We study the new rare B meson decays in Section 3. We subsequently interpret these results in the model introduced before. Finally, we conclude in Section 4.

Generic Lagrangian
Let us extend the SM with a new vector boson V and a new scalar a with masses of the order of TeV and GeV, respectively. They are both singlets of the SM gauge group. At energies below the electroweak scale v ∼ 246 GeV, the Lagrangian we are interested in is where u i L , d i L stand for the i-th family of left-handed up and down quarks, and l i L , e i R stand for the i-th family of left-and right-handed leptons, respectively. Inspired by recently observed flavour anomalies, we will focus mostly on the case in which the only non-vanishing g coupling is g d 23 .
To a lesser extent, we will also consider g d 13 = 0. These couplings are constrained by measurements of ∆M s and ∆M d [27], respectively. Thus, for m V = 1 TeV, g d 23 0.002 and g d 13 is bounded to be about one order of magnitude smaller.
Together with a non-vanishing g , these couplings trigger rare B decays as shown in Fig. 1. The decay width for B 0 s → aa is with f B ∼ 0.23 GeV [28]. Similar expressions hold for other decay modes, e.g. B 0 → aa. The amplitude for B + → K + aa reads: The transferred momentum is q 2 = (p − p 3 ) 2 , and varies between q 2 min = 4m 2 a and q 2 max = (m B − m K ) 2 . The contraction of this matrix element with (p 1 + p 2 ) = q in the amplitude annihilates the f P + (q 2 ) part. Altogether, we obtain In the approximation m K , m a → 0, f 0 (q 2 ) → 1, one easily obtains Following Ref. [29], we parametrize the form factor as f 0 (q 2 ) = r 2 /(1 − q 2 /m 2 fit ), with r 2 = 0.330 and m 2 fit = 37.46 GeV 2 ; see Fig. 2. Similar expressions hold for other processes, e.g. B 0 → K * 0 aa or B + c → D + aa. The latter is however hard to test at the LHCb and we will not consider it. The reason is that the B + c production cross section is much smaller and the B + c width is larger (which reduces both the impact of the new interactions and the experimental efficiency).
We note also that final states containing one meson and aa can probe effective operators containing four quarks and two light scalars (12 of these operators are present in the SM effective field theory extended with a [30].) One can easily estimate Γ(B + → K + aa) ∼ (few GeV) 9 /Λ 8 , which is of the order of 10 −11 provided Λ 1 TeV. For this reason, we will also consider the channel B + → D + aa. It tests operators such as We will assume that a decays into muons with a width smaller than ∼ 10 MeV and a lifetime shorter than ∼ 10 fs. In this case, it will appear not to have any experimentally measurable flight distance and will appear to have zero width. This is easily achieved if a is muonphilic with 10 −5 g e 22 1. The processes discussed so far leads therefore to four-muon final states with and without additional mesons and with the muons forming pairs of two identical masses.

Explicit model
Light scalars are natural within CHMs, for they are approximate Nambu-Goldstone Bosons (NGBs) arising from the spontaneous symmetry breaking G/H in the confinement of a new strong sector at a scale f ∼ TeV. The simplest coset delivering the four Higgs degrees of freedom as well as a new scalar singlet a is SO(6)/SO(5) [31]. Interestingly, it can be UV completed in four dimensions [32].
In this model, the SM fermions do not couple directly to the Higgs boson. They rather mix with other composite resonances that do interact with the Higgs boson. Thus, the Yukawa Lagrangian depends on the quantum numbers of the aforementioned resonances. As a simple yet realistic example, we assume that the second generation leptons mix with two fundamental representations 6 of SO (6). An equivalent description is the embedding of the elementary leptons into incomplete fundamental representations of SO (6). The most general such embedding depends on a single positive parameter γ to give (2.8) Using the corresponding Goldstone matrix , (2.9) one obtains the leading-order Yukawa Lagrangian The coupling y µ ∼ 6 × 10 −4 stands for the muon Yukawa. The subindex in (U T L R ) 1 indicates the projection of the fundamental representation of SO(6) into the singlet of SO(5) according to the decomposition 6 = 1 + 5. The ellipsis stands for terms containing higher powers of a. Likewise, the one-loop induced potential for a reads: where we have neglected terms not involving a. Λ R stands for L R /µ R ; analogously for Λ I L with I the flavour index. c R is a free parameter encoding the details of the strong sectors. Its size can be estimated using naive power counting [33], c R ∼ g 2 * /(32π 2 ), with 1 g * √ 4π the typical coupling between resonances. All in all, we obtain The scalar a decays 100 % into muons. The other fermions respect this phenomenology provided they do not break the shift symmetry a → a + constant, nor a → −a. These two conditions can hold simultaneously if the left (right) chiralities mix with e.g. 6 (1), 6 (15) or 6 (20 ). The scalar defined above can explain the longstanding anomaly on the magnetic moment of the muon [25]. In this concrete model, we can fit the experimental measurement (∆a µ ) obs = (2.74±0.73)×10 −9 [34] within two standard deviations for g * = 2 and f = 800 GeV and γ 10. Fitting the experimental observation within one standard deviation is in principle possible, but it requires even larger values of γ, too small values of g * (which would contradict the strongly coupled nature of the composite sector) and f 800 GeV (which is in tension with Higgs and electroweak precision data [35]). Therefore, the value m a ∼ 1 GeV is a very likely value in this setup.
New composite vector bosons V explaining the observed anomalies in R K ( * ) appear also naturally in this framework [36][37][38][39][40][41][42][43][44][45]. These particles decay preferably into composite states [44], being too broad and too heavy to be directly detected unless very dedicated LHC analyses are performed for masses m V 3 TeV [44]. 1 The interaction with the SM fermions takes the form of Eq. 2.1, with g d 23 ∼ 0.002 m 2 V /TeV 2 [44,46]. Relying on these results, we consider for reference the benchmark point BP : m a = 1 GeV , m V = 4 TeV , g d 23 ∼ 0.03 . (2.14) For g 0.1, this point satisfies all current constraints from LHC searches and measurements of ∆M s . Constraints set by the latter observable could be competitive if the more recent predictions of the SM value are confirmed [46,47]. However, we will show that a signal should be observable with the future upgrades of the LHCb experiment. In particular, we obtain B(B 0 Note that the final state with K + meson is of the same order of magnitude as the four body final state. Together with the fact that the B + meson has a higher production cross section than B 0 s in pp collisions, this suggests that B + → K + µ + µ − µ + µ − is a key signature to explore for this kind of models. This channel has not been experimentally explored though. Similar conclusions were pointed out in Ref. [48] in the context of dark sectors 2 .

Reach of the LHCb
LHCb has searched for the decays B 0 (s) → µ + µ − µ + µ − [26] and has set the limits B(B 0 → µ + µ − µ + µ − ) < 7 × 10 −10 and B(B 0 s → µ + µ − µ + µ − ) < 2.5 × 10 −9 with a 3 fb −1 dataset at the collision energies of √ s = 7 and √ s = 8 TeV. However, there is a limitation in this analysis as it places a veto on the mass of the muon pairs to be close to the ϕ or J/ψ mass, to avoid background from the decay B 0 s → J/ψϕ followed by both vector mesons decaying to a pair of muons. As the a mass is likely very close to the ϕ mass, the current analysis is not sufficiently general. If a new analysis removes the veto around the ϕ mass, and instead imposes a requirement that two opposite muon combinations recombine to the same invariant mass, the limit should stay the same and the background from B 0 s → J/ψϕ would still be avoided. 1 Irrespectively of g , the interaction between V and the heavy resonances L triggers the decay B 0 s → V * → µ + L * , L * → aµ − , where * denotes off-shellness. For couplings equal to the unit, the corresponding width at tree level reads exactly , with m L the mass of L. Even for m L ∼ 500 GeV and m V ∼ 1 TeV, the corresponding branching ratio is ∼ 10 −13 and therefore beyond the reach of our analysis. 2 In this case, rare B meson decays are triggered by flavour-violating scalars.  . Upper right) Same as before but for the transverse momentum of the softest particle. Bottom left) Same as before but for the pseudorapidity of the most central particle. Bottom right) Same as before but for the pseudorapidity of the most forward particle.
Due to the four muons in the final state for the LHCb analysis, the combinatorial background to a possible signal is extremely low. As an essentially background free analysis, even in the far future, the branching fraction limit can be expected to scale inversely with the number of B mesons produced.
In the run periods of Upgrade-I and Upgrade-II of LHCb, the collision energy will be √ s = 14 TeV. As the b cross-section is scaling more or less in direct proportion to the collision energy, the amount of B mesons produced per inverse fb, can be expected to be around a factor 14/7.7 = 1.8 higher compared to the average Run-1 conditions of the LHC. Expectations below will be quoted for the end of LHCb (9 fb −1 ), end of Upgrade-I (50 fb −1 ) and end of upgrade-II (300 fb −1 ). The naive scaling factors, compared to the current Run-1 for these, and taking the different b cross-sections into account, are for LHCb a factor 4.4, for Upgrade-I a factor 29, and for Upgrade-II a factor 180. Thus for the limits on B(B 0 s → µ + µ − µ + µ − ) we should expect 6 × 10 −10 , 9 × 10 −11 and 1.4 × 10 −11 , respectively. This assumes no changes to the trigger or tracking performance in the upgrades of LHCb.
Irreducible backgrounds to the decay have to be considered. The decay B 0 s → ϕϕ with ϕ → µ + µ − is one of these. Using the measured branching fractions [49,50] we get B(B 0 s → (ϕ → µ + µ − )(ϕ → µ + µ − )) = 1.84 × 10 −5 × (2.89 × 10 −4 ) 2 = 1.5 × 10 −12 . As can be seen from the expected limits above, even at the end of LHCb Upgrade II, this is not relevant. For the equivalent decay mode of the B 0 , the measured branching fraction limit for the B 0 → ϕϕ decay is three orders of magnitude below the B 0 s mode and thus even less of a concern. The decay B 0 s → ϕµ + µ − has a measured differential rate of 2.6 × 10 −8 GeV −2 in the region of the squared dimuon mass close to the ϕ mass [17]. Letting the ϕ decay to a muon pair and considering a mass region with width of around 20 MeV, corresponding to a realistic mass resolution, this will give a background at the 10 −13 level and is thus not relevant.
A simplified model for which limits the LHCb experiment can be made for similar decay modes. When comparing different B hadrons, the relative production fractions as measured at √ s = 7 TeV collisions in the LHCb acceptance [51] are taken into account. The relative production fractions are not expected to change significantly with collision energy as they are determined from the fragmentation process. From this, we conclude that the production of B + and B 0 are the same and that the production of B 0 s mesons is a factor 3.7 less common. For the reconstruction in LHCb, it is assumed that the efficiency is 95 % per track inside the fiducial volume defined by the pseudorapidity 2.5 < η < 5.0 and that tracks have a transverse momentum above 0.5 GeV with respect to the beam axis to be reconstructed. For the trigger it is assumed that at least one reconstructed muon should have a transverse momentum above 1.7 GeV. The effect of these criteria is that final states with a larger number of particles have a lower efficiency, both due to the requirement that all tracks have to be reconstructed but also due to that the muons turn softer and the trigger efficiency thus is getting lower; see Fig. 3.
For the positive identification of muons, it is assumed that the efficiency is 100 % for muons with a total momentum above 2.5 GeV. It is assumed that no or only very loose particle identification is required on the charged hadrons. For the D + reconstruction, it is assumed that only the D + → K − π + π + final state is used. This is the easiest decay mode to reconstruct and has a branching fraction of 9.0 % [52]. The final state with a semileptonic decay of the D + could also be considered, in an analysis similar to the B + → µ + µ − µ + ν [53] analysis carried out by LHCb. However, to estimate the reconstruction efficiency of this five charged lepton final state with a neutrino would require a full detector level simulation which is beyond this paper. For the K * 0 case, Table 1. Expected branching fraction limits for the different decays under consideration in units of 10 −11 . It is assumed that the D + meson is only reconstructed in the K − π + π + final state.
we consider the decay into K + π − , which has a branching ratio of ∼ 67 %. All overall efficiencies for a given final state are evaluated relative to the published analysis on the decays B 0 (s) → µ + µ − µ + µ − [26]. As the trigger and main selection are the same for all these decays, this provides a robust normalisation method.
Simulations are carried out using Pythia 8 [54] for the production of B mesons in pp collisions and EvtGen [55] for the subsequent decays. The decays are assumed to be of the type B → (M )aa with a → µ + µ − and with M a possible meson in the final state. The B meson decay is simulated with a flat phase space distribution. If the hadron is unstable, it is decayed to stable particles using the default model in EvtGen and with branching fractions taken from the PDG [52]. A summary of the expected limits that can be set are given in Table 1.
Translated to the plane (g d 23 , m V ) for given values of g , the limits from B 0 s → µ + µ − µ + µ − and B + → K + µ + µ − µ + µ − are compared with those from ∆M s in Fig. 4. Interestingly, we see that scales of several tens of TeV not yet probed by current experiments could be tested in the Upgrade II of the LHCb with our analysis.

Conclusions
We have considered scenarios involving new heavy and flavour-violating vectors V as well as light scalars a. We have shown that these particles give rise to rare B meson decays that are not yet probed. As the preferred mass of the scalar a lies inside the window vetoed by current LHCb searches, namely [950, 1090] MeV, even the simplest decay mode B 0 (s) → µ + µ − µ + µ − is not fully probed. Other decay modes of interest are B + → K + (D + )µ + µ − µ + µ − and B 0 → K * 0 µ + µ − µ + µ − . We have shown that the five-body final state can be as significant as the four-body. Relying on simulations, we have estimated the reach of the LHCb experiment for these processes in the current  Upper right) Same as before but for g = 0.5. Bottom left) Same as before but for the LHCb Upgrade II and g = 0.1 Bottom right) Same as before but for g = 0.5. In all cases, the benchmark point defined in section 2.1 is shown with a star for reference. run and in Upgrades I and II. In the Upgrade II scenario we expect that branching fraction limits in the 10 −11 region can be reached.
Finally, we emphasize again that the decays into a meson and aa are the only sensible probe of different effective operators in the SM effective-field theory extended with a scalar singlet. We therefore encourage the experimental collaborations to consider these processes in future analyses.