$B_s\to\mu^+\mu^-$ as current and future probe of new physics

The rare flavour-changing neutral current decay $B_s\to\mu^+\mu^-$ is among the most important indirect probes of new physics at the LHC, as it is strongly suppressed in the Standard Model, very sensitive to new physics effects, and theoretically exceptionally clean. We present a thorough state-of-the-art analysis of the constraints on new physics from present and future measurements of this decay, focusing on scalar operators. We show model-independently and in concrete new physics models, namely the MSSM and two leptoquark scenarios, that a future precise measurement of the mass-eigenstate rate asymmetry in $B_s\to\mu^+\mu^-$ would allow to disentangle new physics scenarios that would be indistinguishable based on measurements of the branching ratio alone. We also highlight the complementarity between $B_s\to\mu^+\mu^-$ and direct searches in both model classes. Our numerics is based on the open source code flavio and is made publicly available.


I. INTRODUCTION
The decay B s → µ + µ − is the rarest b hadron decay ever observed. As flavour-changing neutral current based on the b → sµ + µ − transition, it is loop-and CKM-suppressed in the Standard Model (SM). Moreover, due to angular momentum conservation, the SM contribution to the amplitude is proportional to the muon mass squared, leading to a strong helicity suppression. This makes the decay extraordinarily sensitive to physics beyond the SM that gives rise to chirality-changing quark flavour violation, lifting the helicity suppression and potentially leading to sizable new physics (NP) effects. At the same time, it is theoretically exceptionally clean for a b hadron decay. In fact, the hadronic matrix element is fully determined in terms of the B s decay constant, which nowadays can be determined to a precision of two per cent with lattice QCD (see e.g. [1,2]). The short-distance contributions are known to next-to-next-to leading order in QCD [3,4] and next-to-leading order in the electroweak interactions [5], resulting in a remaining non-parametric uncertainty of only 1.5% on the branching ratio [6]. Thanks to the helicity structure, no photon-mediated hadronic FCNC contributions, that plague many of the semi-leptonic exclusive decays, are allowed.
After years of searches for this decay with stronger and stronger limits obtained by the Tevatron experiments CDF and D0, the huge amount of b quarks produced at the LHC lead to the first observation of the decay by the LHCb and CMS collaborations [7]. Recently, LHCb even presented the first single-experiment observation [8]. As we will discuss in Sec. III, in the next decade an improvement of the experimental precision by a factor of 7 and of the theoretical precision by a factor of 3 is not unreasonable. This will allow much more stringent constraints on NP, or has the potential to unveil significant deviations from the SM. Nevertheless, by just observing the branching ratio, there are NP scenarios which cannot be probed; for instance, a scenario leading to a different sign of the amplitude but to a similar magnitude will not affect the branching ratio. This is why the existence of an additional observable, the mass-eigenstate rate asymmetry A ∆Γ accessible from the untagged time-dependent decay rate, with complementary dependence on NP contributions, is crucial [9]. A main point of this work is to demonstrate the impact of a future measurement of A ∆Γ on the parameter space of specific NP models, complementary to the branching ratio as well as to direct searches.
The strong sensitivity of B s → µ + µ − to NP stems in particular from the fact that scalar operators of the form (s L b R )(μ R µ L ) or (s R b L )(μ L µ R ) can lift the helicity suppression present in the SM, and B s → µ + µ − is the only process sensitive to these operators in a significant way. This is why we will focus on NP scenarios that generate sizable contributions to these operators. 1 Such contributions typically arise in one of two ways. Either from the exchange of a coloured boson coupling to a lepton-quark current -a leptoquark; or from the exchange of a new, uncoloured scalar particle -i.e. a heavy Higgs boson -with tree-level or loop-induced flavour-changing couplings to quarks. In supersymmetric (SUSY) extensions of the SM, heavy Higgs exchange is particularly interesting as the corresponding scalar amplitudes can be enhanced by tan 3 β [14 -16]. Also extended Higgs sectors without SUSY can lead to sizable scalar amplitudes that are enhanced by tan 2 β [17][18][19][20].
The rest of this article is organized as follows. In Sec. II, we briefly review the observables accessible in the B s → µ + µ − decay. In Sec. III, we discuss the present status of experimental measurements and the prospects for experimental and theoretical improvements. These are then used in Sec. IV for a model-independent analysis of NP parametrized as Wilson coefficients of dimension-six operators, both from present data and from hypothetical future measurements. In Secs. V and VI, we analyze in detail the phenomenology of the B s → µ + µ − decay in the minimal supersymmetric standard model (MSSM), as an example of a model with scalar operators generated by heavy Higgs exchange, as well as in two leptoquark models. In both scenarios, we highlight the complementarity between the branching ratio and A ∆Γ as well as between B s → µ + µ − and direct searches. Sec. VII contains our conclusions.
II. OBSERVABLES IN B s → µ + µ − In this section let us briefly review the observables accessible in the B s → µ + µ − decay. For a more detailed discussion we refer e.g. to [21,22]. Since one cannot expect to measure the helicity of the muons in the foreseeable future, one considers the helicity-averaged decay rate [22], Depending on whether the initial flavour of the B s meson can be determined (tagged), one can then measure the untagged (CP-averaged) time-dependent rate or the CP asymmetry in the decay rate, Since the width difference in the B s system is sizable [23], the decay is characterized by three observables: • The mass-eigenstate rate asymmetry where B H s and B L s denote the heavy and light mass eigenstates of the B s system. As seen above, this asymmetry can be extracted from the untagged rate and it will be of primary interested in this work. For poor statistics, an experimental extraction of this observable is easier in terms of an effective lifetime [9], where The effective lifetime is connected to A ∆Γ via In the SM these observables take the values In general, A ∆Γ can only take values between −1 and +1.
• The time-integrated (and CP-averaged) branching ratio BR(B s → µ + µ − ). It is related to the "prompt" branching ratio, i.e. the branching ratio in the absence of B s -B s mixing, as [9] BR • The mixing-induced CP asymmetry S µµ . The measurement of S µµ requires flavour tagging, therefore large amounts of data are needed to overcome small tagging efficiencies at LHC. In the SM one has S SM µµ = 0. We will not consider S µµ in the following.

III. EXPERIMENTAL STATUS AND PROSPECTS
First evidence for B s → µ + µ − has been found by the LHCb [24] and CMS [25] experiments individually, who subsequently combined their measurements of B s → µ + µ − and their searches for B 0 → µ + µ − from LHC Run 1 data [7]. Recently, the LHCb collaboration presented their LHCb also performed a first measurement of the effective lifetime, albeit still with sizable uncertainties that do not allow significant constraints on A ∆Γ yet. Searches for B s → µ + µ − have also been performed by the CDF [26], D0 [27], and ATLAS [28] experiments, the latter being the most sensitive. For the branching ratio, ATLAS finds from the one-dimensional profile likelihood In our numerical analysis, we use a combination of the two measurements that have found evidence for the decay, namely CMS and LHCb. Since B s → µ + µ − and B 0 → µ + µ − are correlated, it would not be consistent to simply average the two numbers in (10) and (11), but the twodimensional likelihoods have to be combined. While we cannot perform a full combination taking into account correlations between the experiments, we expect a "naive" combination to be sufficient, given statistical uncertainties are still dominant. In Fig. 1, we show the two-dimensional likelihood contours in the plane of BR(B s → µ + µ − ) vs. BR(B 0 → µ + µ − ) for our combination and the individual measurements, showing also ATLAS for comparison. Profiling over BR(B 0 → µ + µ − ), we find the one-dimensional 1σ interval 3 This combination is in reasonable agreement with the SM prediction for the branching ratio, obtained with flavio v0.20.4 [29] using the calculation of [6] with updated hadronic input listed in Tab. I. To arrive at the average (13), we interpolated between the lines of constant likelihood in the experimental plots. We encourage the experimental collaborations to publish their two-dimensional likelihoods in numerical form in the future. We have implemented the interpolated numerical twodimensional likelihoods of the CMS and LHCb measurements as default constraints on BR(B s → µ + µ − ) and BR(B 0 → µ + µ − ) in the public flavio version 0.20.4.
The experimental uncertainty on the B s → µ + µ − branching ratio is presently dominated by statistics and can improve substantially with more data. The LHCb collaboration expects to measure it with a statistical precision of 0.19 × 10 −9 for an integrated luminosity of 50 fb −1 [33], corresponding to their expectation for LHC Run 4, while CMS expects a precision of 0.4 × 10 −9 for 3000 fb −1 [34], corresponding to their expectation for Run 5. An upgrade of LHCb for Run 5 that would allow the experiment to collect an integrated luminosity of 300 fb −1 could reduce the statistical precision, applying naive scaling, to 0.08 × 10 −9 , corresponding to around 2% of the SM prediction.
While the expectations for the branching ratio might be affected by overall systematic uncertainties stemming, for instance, from the B s production fraction or normalization modes, such overall uncertainties cancel in the determination of A ∆Γ . On the other hand, A ∆Γ is suppressed by the small B s -B s lifetime difference. A measurement of the time-dependent decay rate with a precision of 5% (2%), assuming the SM, thus translates into uncertainties on A ∆Γ of 0.8 (0.3). Optimistically assuming that the systematic uncertainties on the branching ratio can be sufficiently reduced, we thus arrive at the following future scenarios, σ exp (B s → µ + µ − ) = 0.19 × 10 −9 , σ exp (A ∆Γ ) = 0.8 , for 50 fb −1 ("Run 4"), σ exp (B s → µ + µ − ) = 0.08 × 10 −9 , σ exp (A ∆Γ ) = 0.3 , for 300 fb −1 ("Run 5"). (15b) On the theory side, the current relative uncertainty is dominated by the uncertainties on the CKM element V cb and the B s decay constant f Bs . On a time scale of ten years necessary for the projected experimental improvements, a lattice computation of f Bs to a precision of 1 MeV seems realistic [35,36]. Similarly a lattice computation of the B → D form factors to a precision of 0.7% might be possible, which would allow to extract V cb from B → D ν decays with this precision, given a sufficiently precise measurement at Belle-II. Assuming that also the remaining non-parametric uncertainties [6] can be reduced by a factor of two to three, a theoretical precision of 1.6% on the SM branching ratio could be reached, competitive with experiment. The theory uncertainty on A ∆Γ is instead negligible compared to the experimental expectations. We thus take for both future scenarios.

IV. MODEL-INDEPENDENT ANALYSIS
Before analyzing specific NP models, in this section we perform a model-independent analysis of new physics effects in B s → µ + µ − . The relevant effective Hamiltonian for the decay valid in any extension of the SM without new particles near or below the b quark mass scale reads with where only C 10 is non-vanishing in the SM. The normalization of O ( ) S,P is chosen to make their Wilson coefficients renormalization group invariant. We thus do not need to specify the renormalization scale for the Wilson coefficients in the following, but note that C ( ) S,P have dimensions of inverse mass in this convention.
The prompt branching ratio (as used in (9)) and A ∆Γ are given in terms of the Wilson coefficients as where the Wilson coefficients appear through the combinations with S = |S| exp(iφ S ), P = |P | exp(iφ P ), and φ NP s is a NP contribution to the B 0 s −B 0 s mixing phase. In the SM one has S = 0 and P = 1.
While the effective Hamiltonian (17) is appropriate for any NP that is heavy compared to the b quark mass scale, the absence of any signal for NP at LHC typically implies that the NP scale is much larger than the electroweak scale. If this is the case, one can express NP effects modelindependently in terms of Wilson coefficients of an effective Lagrangian invariant under the full SM gauge symmetry. At dimension six, the matching of this "SMEFT" [37,38] onto the above weak effective Hamiltonian leads to two relations among the Wilson coefficients [39], In models where the electroweak symmetry is realized non-linearly in the Higgs sector, these relations can be violated [40]. However, in realistic models with strongly-coupled electroweak symmetry breaking where the Higgs is a pseudo Goldstone boson, this breaking is suppressed by v 2 /f 2 , where f is the pseudo Goldstone decay constant. Although formally of order one, this ratio is usually required to be below the few percent level to comply with precision constraints in the electroweak and Higgs sectors. Therefore we expect the relations (22) to hold to good accuracy in most realistic models and we will assume them to hold for the rest of this section.
In models with NP effects in C 10 or C 10 , but not in the scalar operators, B s → µ + µ − plays an important role as it allows to measure the combination C 10 − C 10 with smaller theoretical uncertainties and without dependence on possible NP effects in electromagnetic penguins or semileptonic vector operators (usually called O ( ) 9 ). Here, we want to focus instead on NP effects in the scalar operators, since B s → µ + µ − has the unique property of being strongly helicity suppressed in the SM, but not in the presence of NP in scalar operators. In semi-leptonic decays, only two observables are sensitive to scalar operators with muons in principle: • The "flat term" F H in B → Kµ + µ − is tiny in the SM in the absence of scalar operators [41][42][43]. However it turns out that if relations (22) is satisfied, the observable is not complementary to B s → µ + µ − and not competitive. 4 • The angular observable S c 6 in B → K * µ + µ − [44]. Again however, the effects are too small to be competitive with B s → µ + µ − .
Thus we will restrict ourselves to B s → µ + µ − .
We now proceed to an analysis of NP effects in the two independent Wilson coefficients C S and C S , assuming C ( ) 10 to be SM-like. For the numerics, we used flavio v0.20.4 [29].

A. Constraints on a single Wilson coefficient
We start by considering NP effects in a single real Wilson coefficient C S = −C P . This applies to all models with heavy NP (and linearly realized electroweak symmetry in the Higgs sector) that satisfy the criterion of Minimal Flavour Violation [45], in particular the MFV MSSM to be analyzed in Sec. V. Within this framework, we have to an excellent approximation S 1 − P ≡ A, with a real A. The B s → µ + µ − branching ratio, BR(B s → µ + µ − ), and the mass-eigenstate rate asymmetry, A ∆Γ can therefore be written as There are two regions of parameter space that correspond to a SM-like branching ratio. The first one corresponds to a small new physics contribution A 1 and has A ∆Γ A SM ∆Γ = 1. The second region corresponds to a NP contribution that is comparable to the SM amplitude A 1. This second region of parameter space predicts A ∆Γ −1. While measurements of the branching ratio alone cannot distinguish the two regions, a measurement of A ∆Γ can.
We demonstrate this by performing the following fits of this Wilson coefficient: • A fit to the present branching ratio measurement, • A fit to a future measurement with the projected uncertainties in (15a) and (16) assuming the experimental central values for the branching ratio and A ∆Γ to equal the SM central values ("SM scenario" with A 1), • A similar "future" fit assuming the measured branching ratio to coincide with the SM expectation, but the central value of A ∆Γ to be −1 ("NP scenario" with A ≈ 1), • The previous two fits also for the future scenario in (15b).
For these fits, we have combined the correlated experimental and theoretical uncertainties into a likelihood function depending only on the Wilson coefficient, as described in [46] (and implemented in flavio as FastFit). The result is shown in Fig. 2. We observe that • the current measurement leaves two solutions with equal likelihood 5 , • future measurements of A ∆Γ will be able to completely exclude one of the two solutions, depending on the sign of A ∆Γ .

B. Constraints on a pair of Wilson coefficients
In the leptoquark models to be studied in Sec. VI, both coefficients C S and C S can be generated simultaneously. We thus proceed to analyze simultaneous NP effects in these two coefficients, first assuming them to be real.
Just as for the single Wilson coefficient, we performed five fits: to present data and to two future scenarios with SM-like or opposite-sign A ∆Γ . The results are shown in Fig. 3. We make the following observations.
• The degeneracy from the present branching ratio measurement leads to a ring in the C S -C S plane.
• Future measurements of A ∆Γ can break this degeneracy, but a two-fold ambiguity cannot be resolved by the branching ratio and A ∆Γ .
• The precision of the future scenario of (15b) ("Run 5") is required to obtain disjoint solutions at 2σ.
So far, we have assumed the Wilson coefficients to be real. In fact, being CP-averaged quantities, the branching ratio and A ∆Γ are not very sensitive to the imaginary parts of C S,P and C S,P , which do not interfere with the real SM coefficient. Nevertheless, for sizable complex NP contributions to the scalar Wilson coefficients, the qualitative picture shown in Fig. 3 changes. To demonstrate this, we have performed a full Bayesian fit to the current branching ratio data -using flavio in combination with emcee [47] -allowing real and imaginary NP contributions to C S and C S (still satisfying the SMEFT relations (22)), with flat priors, and varying also the relevant nuisance parameters such as f Bs and y s . The resulting marginal posterior distribution for the real parts is shown in Fig. 4 left. Compared to the case of real NP contributions shown in Fig. 3 and underlayed as dashed contours here for comparison, the region within the "ring" is now filled, as a too small branching ratio from the destructive interference between the SM coefficient and the real parts of the NP coefficients can be compensated by a contribution from the imaginary parts.
The Bayesian fit can also be used to derive posterior predictions for other observables. In Fig. 4 right, we show the posterior predictions for the flat term in B → Kµ + µ − at high and low q 2 , which is very small in the absence of scalar operators, as well as for the angular observable S c 6 5 We do not take into account the 2017 LHCb measurement of the effective lifetime that distinguishes the two solutions at less than half a standard deviation, [8]. The NP scenario is the same as in Fig. 2. in B → K * µ + µ − , which vanishes in the absence of scalar operators. Obviously, the effects still allowed for a completely general variation of the real and imaginary parts of C S and C S satisfying the SMEFT relations (22) are tiny in both cases, demonstrating that these observables cannot compete with B s → µ + µ − now or in the future, as anticipated.
All the plots in this section can be reproduced with flavio v0.20.4 using the scripts provided by us in a public repository [48]. The B s → µ + µ − decay is very well recognized as an important probe of supersymmetric extensions of the SM, as the decay rate can in principle be enhanced by orders of magnitude [14][15][16]. Therefore, the good agreement between SM prediction and the experimental results for BR(B s → µ + µ − ) leads to strong constraints on the parameter space of the minimal supersymmetric standard model. Here we will discuss the impact of the current and expected future B s → µ + µ − measurements in the context of the MSSM with minimal flavour violation, where the only sources of flavour violation are the SM Yukawa couplings. We will make the additional assumption that there are no sources of CP violation beyond the phase in the CKM matrix and that the squark masses are completely flavour blind, i.e. proportional to the unit matrix. In this case, FCNCs are proportional to Higgsino-stop loops. Despite the absence of any new sources of flavour violation, the B s → µ + µ − decay can provide powerful constraints in this scenario.
In the MSSM with MFV and no new sources of CP violation, the NP contributions can be described by a single, real amplitude A as in eqs. (23) and (24). The NP amplitude A is induced by the exchange of heavy scalar and pseudoscalar Higgs bosons, H and A. In the decoupling limit, these states are approximately degenerate, m H m A . Their combined contribution can be concisely written as [49] where Y 0 0.95 is a SM loop function. The above amplitude is proportional to the third power of t β = tan β = v 2 /v 1 , the ratio of the two Higgs vacuum expectation values, and decouples with the mass squared of the heavy Higgs bosons, m 2 A . The various parameters correspond to loop induced non-holomorphic Higgs couplings. We take into account gluino-squark loops, Wino-squark loops, and Higgsino-stop loops. General expressions for these contributions can be found in [49]. Here we will work in the limit where all SUSY masses are degenerate. In this simple benchmark case we have In the more general case, where SUSY masses are not degenerate, the above expressions are modified by O(1) loop functions that depend on the ratios of SUSY masses. The signs of the gluino, wino, and higgsino contributions are set by the relative signs of the µ term and the gluino mass mg, the µ term and the wino mass mw, and the µ term and the stop trilinear coupling A t , respectively. We follow the same sign conventions as in [49]. Note that the above expressions for the parameters do not depend directly on the overall SUSY mass scale. A mild dependence on the SUSY scale enters only through normalization group running, as the parameters α s , α 2 and m t in (27) correspond to MS values at the SUSY scale 6 . Therefore, bounds from direct searches for SUSY particles, in particular squarks and gluinos, can be easily avoided by decoupling, without significantly affecting the B s → µ + µ − phenomenology. The nullresults of existing searches for squarks and gluinos (see e.g. [54][55][56][57][58][59]) constrain the corresponding particle masses in some cases up to almost ∼ 2 TeV. Sensitivities will likely increase further for the future high luminosity scenarios. In the following, we will set the SUSY masses to 5 TeV, which should be safely outside the reach even of the high luminosity LHC.
While such heavy SUSY masses require a fine-tuning of the electro-weak scale at levels of below 1%, they are easily compatible with a 125 GeV Higgs mass, without having to resort to maximal stop mixing. Indeed, for stop masses of 5 TeV and trilinear couplings of the same order, the Higgs mass is in the ballpark of 125 GeV over a broad range of tan β values.
In the following we will concentrate on a specific benchmark scenario, to illustrate the complementarity of the branching ratio and the mass-eigenstate rate asymmetry in probing the MSSM parameter space. We choose sgn(µmg) = sgn(µmw) = −sgn(µA t ) = −1, corresponding to a positive SUSY amplitude A, and fix the absolute value of the stop trilinear coupling A t by demanding that m h = 125 GeV. Using the SusyHD code [60], we find that |A t | is almost independent of the heavy Higgs mass m A . It lies in the range |A t | ∈ (6.2, 7.5) TeV for tan β between 10 and 60, and grows to |A t | 10.5 TeV for tan β = 80. This leads to the following approximate values for the ε factors that enter the B s → µ + µ − amplitude The ranges for b and FC are due to the variation of A t for 10 < tan β < 80. 6 In our numerical analysis we use SM 2-loop RGEs [50][51][52][53] to run these parameters from the electro-weak scale to the SUSY scale. The sensitivity of the current branching ratio measurements to MSSM parameter space is illustrated in Fig. 5. The dark and light green regions correspond to the regions where BR(B s → µ + µ − ) is compatible with the measurements at the 1σ and 2σ level. The white region is excluded by BR(B s → µ + µ − ) by more than 2σ. We observe two distinct regions of parameter space. As expected, there is (i) a broad region for small tan β and large m A corresponding to a NP amplitude A 1, and (ii) a thin stripe for larger values of tan β where A 1 that also agrees well with the measured branching ratio.
In the plots of Fig. 6 we show the m A -tan β plane in the two future scenarios discussed above. While the size of the A 1 region and the A 1 stripe is shrinking with more precise data, the branching ratio measurement alone cannot exclude the A 1 scenario that corresponds to a sizable new physics contribution. The sensitivity of future measurements of the mass-eigenstate rate asymmetry A ∆Γ is also shown in the plots. The blue hatched regions correspond to A ∆Γ < −0.6 (left plot) and A ∆Γ < 0.4 (right plot). We can clearly see that that future measurements of A ∆Γ can cover unconstrained parameter space and fully probe the A 1 region.
Finally, we discuss the complementarity of the B s → µ + µ − observables and direct searches for the heavy Higgs bosons. The main production modes of heavy neutral Higgs bosons H and A in the MSSM are either gluon fusion or, at large tan β, production in association with b quarks. In the parameter regions that we are mainly interested in, namely multi-TeV Higgs bosons and large tan β, we find that the production in association with b quarks is by far dominant.
The corresponding production cross section can be easily obtained by rescaling known SM results where σ bb (H/A) SM is the production cross section of H/A with SM like couplings to b quarks, and the b parameter was already given above. The σ bb (H/A) SM cross section depends only on the mass of the neutral Higgs bosons and we compute it at NNLO using the public code bbh@nnlo [61]. Concerning the heavy Higgs decays, we note that multi-TeV Higgs bosons are sufficiently close to the decoupling limit, such that we can neglect decays of the scalar H to massive gauge bosons W W and ZZ and decays of the pseudoscalar into A → Zh. We also neglect decays into two light Higgs bosons H → hh (which is tan β suppressed) and A → hh (which is non-zero only in the presence of CP violation). In our setup, all other SUSY particles are sufficiently heavy such that "exotic" decays for example into neutralinos H → χ 0 χ 0 , or staus H →τ +τ − are not kinematically open. In this case, the main decay modes are H/A → tt, bb, τ + τ − . For low tan β, the decays to tops dominate. For large tan β one has roughly 90% branching ratio to bb and 10% branching ratio to τ + τ − . We approximate the total decay width as sum of the top, bottom and tau decay widths. The relevant expressions are In the decay to tt, we do not include higher-order non-holomorphic corrections. Those become relevant only for large tan β, where the tt width itself is negligibly small. For the decay widths with SM like couplings Γ(H/A → ff ) SM , we take into account NLO QCD corrections in the heavy Higgs mass limit from [62]. We find The fermion masses and α s in the above expressions are MS values at the scale of the heavy Higgs boson mass. The most sensitive searches for MSSM Higgs bosons look for the τ + τ − final state. Constraints are given on the cross section times τ + τ − branching ratio separately for gluon fusion production and production in association with b quarks. To obtain constraints on the heavy Higgs parameter space we add up the signal cross sections from H and A, σ bb (H) × BR(H → τ + τ − ) + σ bb (A) × BR(A → τ + τ − ) and compare to the latest 13 TeV bounds from [63,64]. For each value of the heavy Higgs mass, we use the stronger of the CMS and ATLAS bounds. The current constraints are shown in Fig. 5 as black hatched regions. We observe that for Higgs masses below ∼ 1.1 TeV, the constraints from the direct searches are stronger than the indirect ones from B s → µ + µ − . For larger Higgs masses, however, B s → µ + µ − covers unconstrained parameter space.
For the future scenarios we expect that the sensitivity of the direct searches will extend to considerably higher masses. We extrapolate the current expected sensitivities by scaling them with the square root of appropriate luminosity ratios. In the left plot of Fig. 6 we show the expected reach of the direct searches with 300 fb −1 of data, in the right plot the expected reach corresponds to 3000 fb −1 of data.
The plots clearly demonstrate the complementarity of the direct searches, the branching ratio measurement, and the measurement of A ∆Γ . All three probes are required to maximise the coverage of the m A -tan β parameter space.
Models with leptoquarks (LQs), i.e. scalar or vector states coupling to a quark-lepton current, have received lots of attention in recent years as they could explain various anomalies in B physics data, including angular observables in B → K * µ + µ − and hints for lepton flavour universality violation in B → Kµ + µ − and B → D ( * ) τ ν (see e.g. [65][66][67][68][69][70][71]). Assuming SM gauge invariance, the possible quantum numbers for LQs are restricted to ten possible representations (of which five with spin 0 and five with spin 1) [72]. In this work we are interested in scenarios that generate the scalar operators contributing to B s → µ + µ − . At tree level, this is only the case for two representations, Due to the strong enhancement of B s → µ + µ − in the presence of scalar operators, we have also investigated whether any of the other eight LQ representations could give rise to the scalar bsµµ operators at the one-loop level. We found that the only non-zero contribution comes from loop The bsµµ Wilson coefficients in the two LQ models in terms of mass-basis couplings. The superscript "NP" denotes the new physics contribution. The normalization factor N is defined in (36).
diagrams including Higgs exchange, making these contributions completely negligible. Thus we will restrict our discussion to the two vector LQs U 1 and V 2 and to tree-level effects in the following. Interestingly, the U 1 LQ is also able to explain the b → s anomalies [68,71]. For the two considered scenarios, U 1 and V 2 , the interaction Lagrangians are given in the G SMinvariant interaction basis by [72] c., (35b) where i, j are generational indices. After EWSB the quark fields have to be rotated into the mass basis in order to get to the effective operators in the weak Hamiltonian. This can be done by defining LQ couplings in the mass basis that one obtains by contracting the above gauge basis couplings with the quark rotation matrices. We will denote the mass basis couplings without the hat, e.g. λ bµ L . In terms of these couplings we find the Wilson coefficients 7 in Tab. II, where the normalization coefficient is given as where v = 246 GeV denotes the SM Higgs VEV. Since these couplings are invariant under G SM , the tree-level coefficients satisfy the SMEFT relations (22). We observe that both scenarios further satisfy the relation while U 1 fulfills C NP 9 = C NP 10 and C 9 = −C 10 and V 2 fulfills C NP 9 = −C NP 10 and C 9 = C 10 .

A. Direct searches for vector leptoquarks
Leptoquarks are color charged and therefore can be pair produced by the strong interaction at the LHC. Assuming that all other couplings of the LQs to SM particles are sufficiently small compared to the strong gauge coupling, pair production through QCD is the dominant production mode. While for scalar LQs, S, SU (3) gauge invariance completely fixes their couplings to gluons, for vector LQs, V µ , there is more freedom and the production cross section becomes more model dependent. Only if the vector LQs are gauge bosons of an extended gauge group, gauge invariance does completely determine the couplings to gluons. In the most general case the interactions read [75,76] with V µν = D µ V ν − D ν V µ and G µν ,G µν are the gluon field strength and the dual field strength, respectively. For vector LQs that are gauge bosons the anomalous couplings vanish, κ =κ = 0. The term proportional toκ breaks CP. We assume CP conservation and neglect this interaction. The anomalous coupling κ is a free parameter. Therefore, we will consider two scenarios: 1. the gauge boson case with κ = 0, 2. the maximally conservative choice where κ is fixed in such a way that it minimizes the production cross section.
The full expressions for the partonic cross sections for vector LQs can be found e.g. in [75]. In our numerical analysis we obtain the LQ production cross section at the LHC by convoluting the partonic cross section with CT14 NNLO parton distribution functions [77] using LHAPDF6 [78]. The decay modes of the LQ are model-dependent as they are determined by the couplings of the LQ to SM fields. We consider the minimal case, where we only allow for couplings that are relevant for bsµµ transitions. Including also other couplings would decrease the branching ratios and therefore weaken the constraints from direct searches.
In the case of the U 1 LQ we consider the decays U 1 → sµ and U 1 → bµ. Invariance under SU(2) L also implies the existence of the LQ coupling to up-type quarks and neutrinos (see eq. (35)). These couplings are related to the bµ, sµ couplings via the CKM matrix Therefore the relevant decay channels are where j denotes a jet coming from any quark but the top. The V 2 LQ is a SU (2) doublet and consists of a charge 4/3 and 1/3 state. Following a line of reasoning similar to the U 1 case, we find the relevant decay modes to be In principle, there can be a mass difference between the two V 2 components after EWSB, such that the heavier state can decay into the lighter emitting a W boson. Any mass splitting has to be proportional to the Higgs vev and therefore should be at most of the order of This implies that the W is off-shell for realistic LQ masses and we can neglect the branching ratio for such a process. For the V 2 LQ also a coupling to two quarks is allowed (rf. eq. (35)), such that a decay into two jets could be important. Such an interaction would violate baryon number and, therefore, induce proton decay which is highly constrained. In the following, we will simply assume baryon number conservation and neglect this type of coupling. Most experimental LQ searches target pair-produced LQs that decay into a quark and a lepton. The relevant signature is jj (or jj + / E T if one LQ decays into a neutrino). A common assumption in the experimental searches is that the LQ couples only to one generation of SM particles. However, the searches for first and second generation LQs usually do not apply a b-veto such that also the bµ final states relevant for this work are covered. A last caveat is that most searches implicitly assume scalar LQs. In the case of vector LQs, the angular distributions will be different, leading to slightly different acceptances. This effect has been found to be small in [79] and we will neglect it in the following. In our analysis we include the LQ searches from ATLAS [80,81] and CMS [82] using 8 TeV and 13 TeV data.

B. Leptoquark effects in
The LQ scenario U 1 can explain present anomalies in semi-leptonic b → sµµ transitions [68] if C 9 = −C 10 ≈ −0.5 and C 9 = C 10 ≈ 0 [46,83,84]. In the following we will focus on this case. We fix the left-handed couplings such that C 9 and C 10 take the above values and keep the right-handed couplings as free parameters. By (37) we are then forced to set either the unprimed or primed scalar operator to zero. We decide to consider the following benchmark scenario: 8 which results in C 9 = −C 10 = −0.5, C 9 = C 10 = 0, C S = −C P (free), C S = C P = 0. The LQ V 2 predicts C NP 9 = +C NP 10 and therefore cannot explain the b → sµµ anomalies. We assume degenerate masses for the two components with charges 1 4 and 4 3 , and consider the following scenario: This setup corresponds to a rather small value of C 9 = +C 10 = 0.1, vanishing primed Wilson coefficients, and free C S = −C P . In Fig. 7 we show the presently allowed parameter regions for both scenarios in the masscoupling plane. In the parameter space still allowed by experimental searches one finds two solutions that reproduce the current measurement of the branching ratio. One is SM-like while the other corresponds to large NP effects.
Projecting into the future, a measurement of A ∆Γ will be able to disentangle this situation. We find that the NP solution gives rise to a large negative A ∆Γ , such that already an estimation of its sign can rule out this scenario. In Fig. 8 we present future projections for the mass vs. coupling  (43) and (44). Inside the dark and light green bands, the present value of the experimental branching ratio (13) is reproduced at 1 and 2σ, respectively. The black //-hatched regions show the exclusions from present direct searches. The more densely hatched region corresponds to minimal LQ production, while the more coarsely hatched region is for YM-like production.
plane. We consider the Run 4 and 5 of the LHC. We estimate the direct constraints by rescaling the present exclusion limits by the square root of the luminosity ratios, L today /L future , where L future = 300 fb −1 , 3000 fb −1 for Run 4 and Run 5, respectively. In high mass ranges, for which presently there are no direct constraints, we conservatively extrapolate the current exclusion limits as constants. For the branching ratio we assume the current SM value (14) with uncertainties (15). These projections clearly show the power of an A ∆Γ measurement in eliminating a degenerate solution as well as the general impact of B s → µ + µ − on the LQ parameter space. The LQ representation V 2 also contributes to B → K ( * ) νν [85]. Using the notation of that reference we find that B s → µ + µ − is correlated to the b → sνν transition via m 2 As NP contributions to C 9 = −C 10 are strongly constrained, one cannot expect large contributions to B → K ( * ) νν. Taking a benchmark value of C 9 = −C 10 = 0.33 (which corresponds to the maximally allowed value at the 2σ-level found in [46]) and assuming degenerate masses for the two LQ components, we only find a modification of about 3% relative to the SM expectation. This, however, is only true if the LQ solely couples to muons. If also couplings to other lepton generations are allowed, then the effects can be much larger [85].
In summary, we have shown that a measurement of A ∆Γ can efficiently reduce degeneracies in the LQ parameter space that cannot be resolved by a measurement of the branching ratio alone. Interestingly, even the two-fold ambiguity remaining in the presence of simultaneous real NP effects in C S and C S as shown in Fig. 2 might be resolved in the LQ scenarios due to the relation (37) that predicts simultaneous effects in C 9 and C 9 as well as C 10 and C 10 for this solution, which could be tested in semi-leptonic B decays such as B → K * µ + µ − . Whether their sensitivity will be strong enough to really resolve this ambiguity is an interesting question but beyond the scope of our present study.  (43) and (44). The first row is for the "Run 4" scenario while the second row marks "Run 5". The green 1σ-and 2σ-regions correspond to the anticipated future experimental sensitivity of the branching ratio, assuming the SM central value (14). The black //-hatched regions show the extrapolated exclusions from direct searches. The more densely hatched region corresponds to minimal LQ production, while the more coarsely hatched region is for YM-like production. The blue \\-hatched region would be excluded at the 2σ level by a measurement of A ∆Γ with SM-like central value.

VII. CONCLUSIONS
In this work we discussed the new physics sensitivity of present and future measurements of the B s → µ + µ − decay.
The B s → µ + µ − decay is a rare flavour-changing neutral current process and widely recognized as one of the most important indirect probes of new physics at the LHC. The present experimental measurements of the branching ratio are compatible with the SM prediction and have uncertainties at the level of 20%. With data taken in future runs of the LHC, the experimental precision can improve significantly, potentially reaching ∼ 5% with 50 fb −1 in Run 4 at LHCb, and 2% with 300 fb −1 in Run 5 at LHCb, respectively. The SM prediction is expected to reach similar precision within one decade. In addition to high precision measurements of the branching ratio, the large statistics that will become available in future LHC runs will also allow the experiments to measure additional observables in the B s → µ + µ − decay. In particular, a measurement of the masseigenstate rate asymmetry A ∆Γ with a precision of 80% in Run 4 and 30% in Run 5 should become possible.
We analysed the interplay of the branching ratio and the mass-eigenstate rate asymmetry in constraining new physics both model-independently and in the context of new physics models, namely the minimal supersymmetric standard model and two leptoquark models. Our analysis shows that future measurements of the mass-eigenstate rate asymmetry would allow to disentangle new physics scenarios with scalar interactions, that would be indistinguishable based on measurements of the branching ratio alone.
Under mild assumptions, scalar contributions to the B s → µ + µ − decay can be modelindependently described by two complex parameters. Fits to branching ratio data alone leave a large degeneracy in parameter space. Assuming CP conservation, the degeneracy can be broken by A ∆Γ up to a two-fold ambiguity. We provide all the tools necessary to reproduce the plots from our model-independent analysis in Sec. IV with flavio in a public repository [48].
In the context of the MSSM, we use the measurements of BR(B s → µ + µ − ) to obtain constraints in the m A -tan β plane. We find that the branching ratio alone is not sensitive to a region of parameter space with a sizable SUSY contribution to B s → µ + µ − predicting opposite sign but similar magnitude for the transition amplitude. Future measurements of A ∆Γ can cover this region. We also highlight the complementarity with direct searches for heavy Higgs bosons. While the direct searches put the strongest constraints at low Higgs masses 1 TeV, the indirect constraints from B s → µ + µ − extend the coverage far into the multi-TeV region.
We make similar observations also in the considered leptoquark models. There are two leptoquark representations that can mediate scalar operators contributing to B s → µ + µ − . In both models, and similarly to the MSSM, there are regions of parameter space with sizable contributions to the B s → µ + µ − amplitude that lead to a SM-like branching ratio. Measurements of A ∆Γ can cover much of these regions, although ambiguities remain in the most general case. Direct searches for leptoquarks exclude most regions of parameter space for leptoquark masses of around 1 TeV and below. The B s → µ + µ − observables are able to indirectly probe much heavier leptoquarks.
Our work explicitly demonstrates the complementarity of the B s → µ + µ − branching ratio and mass-eigenstate rate asymmetry in probing new physics and stresses the relevance of the large statistics that can be expected from future high luminosity runs at the LHC.
The work of CN and DS was supported by the DFG cluster of excellence "Origin and Structure of the Universe".