CP Violation in Leptonic Rare $B^0_s$ Decays as a Probe of New Physics

The decay $B^0_s\to\mu^+\mu^-$ is a key probe for the search of physics beyond the Standard Model. While the current measurements of the corresponding branching ratio agree with the Standard Model within the uncertainties, significant New-Physics effects may still be hiding in $B^0_s\to\mu^+\mu^-$. In order to reveal them, the observable $\mathcal{A}^{\mu\mu}_{\Delta \Gamma_s}$, which is provided by the decay width difference $\Delta\Gamma_s$ of the $B^0_s$-meson system, plays a central role. We point out that a measurement of a CP-violating observable ${\cal S}_{\rm \mu\mu}$, which is induced through interference between $B^0_s$-$\bar B^0_s$ mixing and $B_s\to\mu^+\mu^-$ decay processes, is essential to obtain the full picture, in particular to establish new scalar contributions and CP-violating phases. We illustrate these findings with future scenarios for the upgrade(s) of the LHC, exploiting also relations which emerge within an effective field theory description of the Standard Model, complemented with New Physics entering significantly beyond the electroweak scale.


Introduction
The decay B 0 s → µ + µ − is one of the most interesting processes offered by Nature, allowing us to test the Standard Model (SM) and probe New Physics (NP). In the SM, this channel has no contributions at the tree level and shows a helicity suppression [1]. Consequently, the SM branching ratio is enormously suppressed, and only about three out of one billion B 0 s mesons decay into the µ + µ − final state. Another key feature of B 0 s → µ + µ − is related to the impact of strong interactions. As gluons do not couple to the leptonic final state, only the B 0 s decay constant f Bs enters the theoretical description, which can be calculated by means of lattice QCD [2].
As NP effects may enhance the branching ratio of B 0 s → µ + µ − significantly, experiments have searched for this channel for decades [3]. It has been a highlight of the results of the Large Hadron Collider (LHC) that B 0 s → µ + µ − could eventually be observed by the CMS and LHCb collaborations and is now experimentally well established [4], with a measured branching ratio in the ballpark of the SM prediction. In addition to the branching ratio, B 0 s → µ + µ − offers another observable, A µµ ∆Γs , which is accessible thanks to the sizeable decay width difference ∆Γ s of the mass eigenstates of the B 0 s -meson system [5]. This observable is theoretically clean and plays an important role in the search for NP effects [6][7][8]. A pioneering measurement of A µµ ∆Γs has recently been reported by the LHCb collaboration [9]. This analysis requires, in contrast to the measurement of the branching ratio, time information for untagged B s data samples.
If also tagging information is available, a CP-violating observable S µµ can be measured which arises from the interference between B 0 s -B 0 s mixing and decay processes. Should it be possible to determine the helicity of the final-state muons, yet another CP asymmetry C µµ can be measured, as discussed in detail in Refs. [5,6]. It is not independent from A µµ ∆Γs and S µµ , as the observables satisfy the following relation: A µµ ∆Γs 2 + (S µµ ) 2 + (C µµ ) 2 = 1 .
In these observables, as in the case of A µµ ∆Γs , the decay constant f Bs cancels. Consequently, they are theoretically clean. Within the SM, the CP asymmetries vanish. However, in the presence of physics beyond the SM, we may in general encounter new sources of CP violation, generating non-vanishing CP asymmetries and affecting also the observable A µµ ∆Γs . In analyses of rare B (s) decays, it is usually -for simplicity -assumed that CPviolating NP phases vanish. Within specific models, such assumptions can be made, where an important example is given by scenarios with "Minimal Flavour Violation" [10]. However, we would rather like to learn from experimental data whether new CP-violating phases enter the dynamics of the decay B 0 s → µ + µ − . In this paper, we explore this question. Interestingly, we find that S µµ is an essential observable to reveal the nature of possible NP effects. The sign of the CP asymmetry C µµ would allow us to resolve certain ambiguities. We shall illustrate these findings with various examples, showing in particular how we may establish new (pseudo)-scalar contributions to B 0 s → µ + µ − and further resolve their structure and dynamics. These considerations are completely general and can also be applied to the rare B 0 s → τ + τ − and B 0 s → e + e − decays [8]. The outline of this paper is as follows: in Section 2, we discuss the theoretical description of B 0 s → µ + µ − and introduce the corresponding observables. In Section 3, we explore then the situation with general CP-violating NP contributions. Assuming relations between short-distance coefficients, which are motivated by considerations within effective field theory, we analyze the interplay between the B 0 s → µ + µ − observables in Section 4. In Section 5, we shall address experimental aspects by discussing scenarios and illustrating their physics reach by making assumptions about the experimental precision. Finally, we summarize our key results and give a brief outlook in Section 6.

Decay Amplitude
The theoretical framework to describe the decayB 0 s → µ + µ − is given by effective quantum field theory, which allows the calculation of a low-energy effective Hamiltonian of the following general structure [1,5,7]: Here G F is Fermi's constant, V * ts V tb is a factor with elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix, and α denotes the QED fine structure constant. The Wilson coefficients C ( ) 10 , C ( ) P and C ( ) S describe heavy degrees of freedom, which have been integrated out from appearing as explicit fields, and are associated with the four-fermion operators O 10 = (sγ µ P L b)(μγ µ γ 5 µ), with m b denoting the b-quark mass, and In general, the Wilson coefficients are different for b → s and b → d transitions, and depend on the flavour of the final-state leptons [8]. For simplicity, we do not give the corresponding labels explicitly in the following discussion. In the SM, we have only to deal with the O 10 operator, having a real coefficient C SM 10 . Introducing the combinations of Wilson coefficients where M Bs , m µ , m b , m s are the corresponding particle masses and ϕ P , ϕ S denote CPviolating phases, we obtain the following expression for the decay amplitude [5]: Here λ = L, R describes the helicity of the final-state leptons with η L = +1 and η R = −1.
In the SM, we have and the relevant Wilson coefficient is given as [6] where η Y describes QCD corrections, θ W is the weak mixing angle, Y (x t ) represents one of the Inami-Lim functions, and x t ≡ m 2 t /M 2 W parametrizes the top-quark and W mass dependence [11]. We would like to emphasize that, by convention, C SM 10 does not have a complex phase. However, it takes a negative value, such that In the following discussion, the CP-violating phases ϕ P and ϕ S play a central role. While the latter is directly related to the phase of the short-distance coefficient C S − C S of new scalar contributions, the former may get contributions both from C 10 − C 10 and from the coefficient C P − C P , which arises from new pseudo-scalar operators.

Branching Ratio and Effective Lifetime
Due to B 0 s -B 0 s mixing, an initially, i.e. at time t = 0, present B 0 s meson evolves into a time-dependent linear combination of B 0 s andB 0 s states. For the "untagged" rate no "tagging" of the initially present B s meson is needed. This quantity depends only on two exponentials and involves the parameter which characterizes the decay width difference of the B s mass eigenstates, with τ Bs ≡ 1/Γ s denoting the B s mean lifetime [12,13]; for the experimental value, see Ref. [14]. The decay dynamics enters through the following observable [5,6]: which is independent of the muon helicity, as reflected by the definition of A µµ ∆Γs . Within the SM, we have A µµ ∆Γs | SM = +1.
The phase φ NP s originates from possible CP-violating NP contributions to the B 0 which is already strongly constrained by experimental data for CP-violating effects in B 0 s → J/ψφ and decays with similar dynamics, yielding the following results [14][15][16]: where we have used the SM value φ SM s = −2β s = −(2.12 ± 0.04) • in Eq. (17). Since it is challenging to measure the muon helicity, we consider the helicity-summed rates and use them to define an untagged rate Γ(B s (t) → + − ) in analogy to Eq. (11). The branching ratio reported by experiments actually corresponds to the following timeintegrated untagged rate [5,12]: Combining the CMS result from 2013 [17] with the most recent LHCb analysis [9] yields This average was calculated by means of the Particle Data Group (PDG) procedure [18]. For comparison, we give also the constraint B(B s → µ + µ − ) ATLAS'16 = (0.9 +1.1 −0.8 ) × 10 −9 reported by the ATLAS collaboration [19].
In the SM, we have the following expression [1]: where special care has to be taken concerning the use of renormalization schemes to properly include next-to-leading-order electroweak corrections (for details, see Ref. [1]). Using current state-of-the-art input parameters yields the following result [8]: In a very recent analysis [20], QED corrections from dynamics below the renormalization scale µ = m b were calculated, affecting the branching ratio by almost 1%. In order to search for NP effects by means of the branching ratio of B 0 s → µ + µ − , the following ratio plays the key role [5,6]: taking by definition the SM value R| SM = 1.
Using the expressions given above yields with The numerical results in Eqs. (21) and (23) give The effective lifetime of the decay B 0 s → µ + µ − , which is defined through contains the same physics information as the observable A µµ ∆Γs [5]: A pioneering measurement of the effective lifetime of B 0 s → µ + µ − was recently reported by the LHCb collaboration [9]: Using Eq. (30), this result can be converted into where the error is fully dominated by the uncertainty of τ s µµ . In view of the general model-independent range − 1 ≤ A µµ ∆Γs ≤ +1, it will be crucial to improve the experimental precision for this observable at the LHC upgrade(s) in order to use this quantity for testing the flavour sector of the SM.

CP Asymmetries
In contrast to the untagged B s rate in Eq. (11), the tagged, time-dependent rates involve oscillatory sin(∆M s t) and cos(∆M s t) terms, where ∆M s is the mass difference between the heavy and light B s mass eigenstates. We obtain a CP-violating rate asymmetry of the following form [5,6]: with the observables where η L = +1 and η R = −1 for left-and right-handed muon helicity, respectively. It should be noted that the CP asymmetry S λ µµ , which is caused by interference between B 0 s -B 0 s mixing and B s → µ + µ − decay processes, does actually not depend on the muon helicity, just as the observable A µµ ∆Γs ≡ A λ ∆Γs . Using the helicity-summed rates introduced above yields where the C λ µµ terms cancel because of the η λ factor. It should be noted that a nonvanishing C µµ would be a smoking-gun signal for a new scalar contribution S. CPviolating asymmetries of this kind in B s,d → + − decays were also considered for various NP scenarios in Refs. [21][22][23], neglecting the effects of ∆Γ s and the associated observable A µµ ∆Γs . For a more recent study, including the untagged observable, see Ref. [6]. It should be stressed that the non-perturbative decay constant f Bs cancels in A µµ ∆Γs as well as in S µµ and C µµ , thereby making these observables theoretically clean probes for the search of NP signals [5,6]. In the SM, a tiny residual uncertainty arises from QED corrections, which lead to effects at the 10 −5 and 10 −3 levels for A µµ ∆Γs and the CP asymmetries S µµ , C µµ , respectively [20].
In the following discussion, we will explore the interplay of A µµ ∆Γs and S µµ with the observable R to search for NP and reveal its nature, in particular whether it involves new (pseudo)-scalar contributions. Experimental feasibility studies of measurements of the CP asymmetry in Eq. (37) have not yet been performed to the best of our knowledge. However, we envision that an effort should be made to perform such a measurement at the LHC upgrade(s). In view of the relation in Eq. (1), a measurement of C µµ would not provide independent information. As such an analysis would require the reconstruction of the muon helicity, it is much more challenging than the asymmetry in Eq. (37) involving the helicity-averaged rates. However, we will show that already information on just the sign of C µµ would be sufficient to resolve certain ambiguities affecting the determination of P and S. We encourage experimentalists to explore avenues to eventually measure the sign of the C µµ observable.
3 General CP-Violating New Physics

Theoretical Description
Let us start the general discussion of the CP-violating coefficients P and S in Eqs. (5) and (6), respectively, with the ratio R in Eq. (24). Using the expression in Eq. (26), we obtain If we had a precise measurement of A µµ ∆Γs , we could straightforwardly convert R into r. In view of the large uncertainty in Eq. (32), we use the general range in Eq. (33) to calculate 0.69 ≤ r ≤ 1.13, where we have also taken into account the 1σ uncertainty of R, given in Eq. (28). This observable fixes a circular band with radius √ r in the |P |-|S| plane, which we show in where we have introduced the abbreviations If we assume that the CP-violating phases ϕ P and ϕ S take trivial values, i.e. 0 • or 180 • , R allows us to fix a circle in the |P |-|S| plane through Eq. (26), and the intersection with the straight line following from fixes |P | and |S|, as discussed in detail in Refs. [5][6][7][8]; note that we use the result for φ NP s in Eq. (17). However, if we allow for general CP-violating phases, any point on the circle with radius √ r is allowed since we obtain |S| = 0 for cos Φ P = A µµ ∆Γs and |P | = 0 for cos Φ S = −A µµ ∆Γs .
The measurement of a non-vanishing CP asymmetry S µµ would immediately establish the presence of non-trivial CP-violating phases. This observable fixes another straight line in the |P |-|S| plane: Figure 2: Flowchart to illustrate the general strategy to determine |P | and |S| as functions of the the CP-violating phase ϕ S from the B 0 s → µ + µ − observables.
However, as the CP-violating phases are in general unknown, the slope of this straight line is not determined, in analogy to the constraint following from A µµ ∆Γs . We have three independent observables at our disposal, r as well as A µµ ∆Γs and S µµ , which depend on the four unknown parameters |P |, Φ P and |S|, Φ S . Using the general expressions for A µµ ∆Γs and S µµ in Eqs. (13) and (36), respectively, yields This equation allows us to determine Φ P as a function of Φ S with the help of The expression under the square root is actually factorizable, thereby yielding Using then the observables r and A µµ ∆Γs , we may determine as functions of the CP-violating phase Φ S . Using instead of A µµ ∆Γs the CP asymmetry S µµ yields The expression in Eq. (49) leaves us with a twofold ambiguity for ϕ P for every value of ϕ S . Information on the sign of C µµ allows us to determine the correct branch and thus obtain a single solution for ϕ P as a function of ϕ S . However, both branches have the same dependence of |P | and |S| on ϕ S , so a single solution for |P | and |S| as a function of ϕ S can be obtained even when no information on the sign of C µµ is available. In the flowchart in Fig. 2, we illustrate this general method for analyzing the observables provided by the B 0 s → µ + µ − decay, and we will provide an example of this formalism in the next subsection.

Vanishing Mixing-Induced CP Violation
An interesting situation arises for S µµ = 0. Although one may naively conclude that the CP-violating phases take then simply trivial values, this is actually not the case because of the structure of the expression in Eq. (36). In fact, we obtain the following extremal values on the circle with radius √ r in the |P |-|S| plane: where the region between these points can be accessed by varying Φ S . In the case of A particularly interesting situation arises for A µµ ∆Γs = 0, corresponding to the following point in the |P |-|S| plane:

Sizeable Mixing-Induced CP Violation
Let us now turn to mixing-induced CP violation in B 0 s → µ + µ − , and discuss a scenario with a large value of S µµ , which requires significant CP-violating phases originating from physics beyond the SM. In order to illustrate this situation and the formalism discussed in Subsection 3.1, we consider an example which is characterized by Assuming furthermore These values of |P | and |S| fall well within the currently allowed region in the |P |-|S| plane shown in Fig. 1. We obtain the following set of observables: and assume that they were measured at a future experiment.
Let us now illustrate how we may obtain insights into NP effects using these observables. The deviation of A µµ ∆Γs from the SM prediction +1 would indicate NP effects. Having the measured A µµ ∆Γs at hand, we may use Eq. (38) to convert R into r, yielding Moreover, the precision of the measured B 0 s → µ + µ − branching ratio will then have significantly increased (see Section 5 for a more detailed discussion), allowing us to reduce the width of the circular band in Fig. 1. However, without any information on S µµ , we could not narrow down further |S| and |P | in a model-independent way, i.e. we would still be left with the whole circular region, and could in particular not establish a non-vanishing scalar contribution S.
The measurement of the observable S µµ different from zero would signal new sources of CP violation. Using then Eq. (49), we could determine ϕ P as a function of ϕ S , as illustrated in Fig. 3. The information on the sign of C µµ would allow us to resolve the ambiguity, as indicated in the figure. Note that the points (ϕ S , ϕ P ) = (0 • , 0 • ) and (180 • , 180 • ) would be excluded through the contours. Using Eqs. (51) or (52), we obtain |S| and |P | as functions of ϕ S , as shown in Fig. 4. Here, information about the sign of C µµ plays no further role. Interestingly, we would now be able to put a lower bound on |S|, i.e. could conclude that we have new scalar contributions. We insist on the fact that in order to obtain this highly non-trivial information, a measurement of the CP asymmetry S µµ is required. Although we can only determine the B 0 s → µ + µ − parameters as functions of ϕ S , this analysis would have profound implications, establishing in particular new scalar and pseudo-scalar contributions with CP-violating phases. In order to obtain further insights, more information is needed and assumptions about short-distance coefficients have to be made.

General Framework
The effects of new particles enter the coefficients in Eqs. (5) and (6) through the shortdistance coefficients C P , C P and C S , C S , which describe new pseudo-scalar and scalar contributions, respectively, and C 10 , C 10 . As the constraints from the ATLAS and CMS experiments at the LHC for direct searches of new particles support the picture of a NP scale Λ NP which is much larger than the electroweak scale Λ EW , the corresponding NP effects can be described in a model-independent way through an effective Lagrangian where the heavy degrees of freedom, i.e. the NP particles, have been integrated out at Λ NP . If we require then invariance under the SM gauge group SU (2) L × U (1) Y for the renormalization group evolution between Λ NP and Λ EW , a "SM Effective Field Theory" (SMEFT) can be set up [24,25] and matched to the effective Hamiltonian in Eq. (2) describing B 0 s → µ + µ − decays. Following these lines and applying the machinery of effective quantum field theory, the following relations among the corresponding shortdistance coefficients can be derived [26]: A further application of these relations -assuming no new sources of CP violation -can be found in Ref. [7], while a fit of data to the SMEFT scenario with complex coefficients was performed in Ref. [27]. For a discussion within specific models, see Ref. [6].
In this section, we explore the implication of Eqs. (63) and (64) for the general analysis of CP violation discussed in Section 3. To this end, we express the relevant quantities in terms of the scalar short-distance coefficients which yields and We will refer to this notation as the SMEFT parametrization. It it useful to write Eq. (66) in the following form: In the Appendix, we present expressions that allow us to obtain the B 0 s → µ + µ − observables in terms of the parametrization introduced above. As P requires input for C 10 , C 10 , we shall now first discuss these coefficients.

Closer Look at C 10 and C 10
The Wilson coefficients C 10 and C 10 enter in P through the following combination: where ϕ 10 is a CP-violating phase and parametrizes NP effects. The relations allow us to express C 10 in terms of the -in general -complex NP coefficient C NP 10 . In order to reveal the substructure of P , information on C 10 is required. In specific models, we may calculate C NP 10 (see, for instance, Ref. [6]). Alternatively, using experimental data for B → K ( * ) + − decays, we may determine C 10 − C 10 from experiment (see Ref. [28] and references therein). In practice, the corresponding NP contributions are extracted through involved global fits to sets of large numbers of observables. We use the results from Ref. [28], where different scenarios for NP in real Wilson coefficients are discussed. Considering NP in individual Wilson coefficients, the authors find that the data is best explained by a contribution to the short-distance coefficient C 9 of the four-fermion operator O 9 = (sγ µ P L b)(μγ µ µ), which does not contribute to B 0 s → µ + µ − , yielding C NP 10 = 0 and thus C 10 = 1. However, a similarly good fit is obtained by assuming the relation C NP 9 = −C NP 10 for real coefficients, which appears in models with new particles that couple only to left-handed leptons. In this case, we find where the minus sign follows from C SM 10 taking a negative value, as given in Eq. (9), resulting in In Ref. [28], CP-violating phases are neglected. However, the short-distance coefficients are in general complex, and the phases can be included in the fit. In Ref. [29], such an analysis is performed. The results are presented as 2D confidence contours in the complex plane of the coefficients C 10 and C 10 . To probe for the possible size of ϕ 10 and |C 10 |, we assume that C 10 = 0 and convert the 1σ allowed regions for the complex Wilson coefficient C 10 shown in Ref. [29] into C 10 using Eq. (70), yielding Due to the structure of Eq. (74), we obtain a rather constrained range for the CPviolating phase ϕ 10 . It is also interesting to note that the range for the absolute value |C 10 | is consistent with the result in Eq. (76).
In the future, analyses of CP-violating effects in B → K ( * ) + − and B s → φµ + µ − observables, as introduced in Refs. [30,31], will allow us to get a much sharper picture of |C 10 | and a possible complex phase ϕ 10 . It would be very useful to add the complex coefficient C 10 as a default output to the corresponding sophisticated fits to the semileptonic rare B (s) decay data.
For the numerical illustrations below, we will either use the range in Eq. (76) for real Wilson coefficients C 10 and C 10 , or we will consider the case |C 10 | = 1, ϕ 10 = 0 • , where NP effects would enter exclusively through (pseudo)-scalar contributions.
An interesting situation arises if we consider a scenario where NP effects enter only through C 10 , with vanishing coefficients C P , C P and C S , C S , yielding P = C 10 and S = 0. Specific examples are given by models with extra Z bosons (see, for instance, Ref. [6]) and scenarios with modified Z couplings (such as in models with vector-like quarks [32]). We would then have the simple expressions with C µµ = 0. Consequently, the observables would lie on a circle with radius one in the A µµ ∆Γs -S µµ plane.

Extraction of |x| and ∆
Applying the method presented in Subsection 3.1, we may determine |S|, |P | and ϕ P from the B 0 s → µ + µ − observables as functions ϕ S . Using Eq. (66), we may convert these parameters into the ratio x of the -in general -complex scalar short-distance coefficients: Sign of C µµ

SMEFT relations
General analysis with and yielding The quantities entering these expression can be expressed in terms of the absolute values and phases of the relevant complex coefficients as It is instructive to consider the example in Subsection 3.2.2, where |S| = 0.30 and ϕ S = 20 • . Using the expressions given above, we can convert the corresponding values of |P | = 0.89 and ϕ P = 30 • into where we have assumed no NP in C 10 , so |C 10 | = 1 and ϕ 10 = 0 • . In Fig. 5, we give a flowchart for this strategy, and show in Fig. 6 the situation corresponding to Eq. (87). Using information on the sign of C µµ , we would only be left with the red contours. We observe that |x|e i∆ could be constrained in a very non-trivial way. The resulting contours depend strongly on the associated B 0 s → µ + µ − decay observables. In order to constrain the parameters more stringently, it is useful to make assumptions about scenarios, as we will illustrate in the next section. Following these lines, we may rule out a given scenario or confirm it, allowing us then to extract the corresponding parameters. By the time we may have measurements of CP violation in B 0 s → µ + µ − available, we should have a much better picture of the physics beyond the SM, thanks to The grey contours could be excluded through sign information for the observable C µµ . the interplay between model building and data coming both from the high-energy and the high-precision frontiers. In particular, we should then also have some preferred scenarios, including specific patterns for the CP-violating phases, which could be confronted with experimental data and the new strategies presented in this paper.

Illustration
As experimental data have already constrained the NP contribution φ NP s to the B 0 s -B 0 s mixing phase to be tiny, as given in Eq. (17), we may simplify the discussion by neglecting this quantity. Moreover, for the decay B 0 s → µ + µ − , we have with excellent precision w = 1. Let us now illustrate the formalism and strategy discussed above through various examples. Here we shall choose values for the input parameters to calculate the decay observables. Assuming then that these quantities have been measured at the future LHC upgrade(s), we discuss the pictures emerging from the strategy discussed above. For simplicity, we do not consider experimental aspects in this section but will illustrate scenarios assuming uncertainties of future measurements in Section 5.
The number of allowed solutions for a given angle ϕ S depends on the value of the Wilson coefficient C 10 . In order illustrate this feature, we consider two scenarios for C 10 . Let us first assume that there is a vanishing NP contribution C NP 10 = 0, which yields |C 10 | = 1, ϕ 10 = 0 • . In this case, Eq. (93) results in two solutions for |S| as a function of ϕ S , as can be seen in the top-left plot in Fig. 7. Using Eqs. (90), (91) and (92), we can determine the observables A µµ ∆Γs , S µµ and C µµ as functions of ϕ S , respectively, as shown in Fig. 7. In particular, once A µµ ∆Γs has been measured, the value of S µµ can be predicted. Should this CP asymmetry be measured correspondingly, this scenario would be confirmed, allowing us to determine the corresponding NP parameters. On the other hand, should the measurement of S µµ be in conflict with the prediction, the NP scenario would be ruled out by experimental data.
Let us now consider a scenario with NP contributions to C 10 . If we follow the analysis of Ref. [28] and use the central value of C 10 in Eq. (76), we obtain the functional dependence of |S| and the corresponding observables on ϕ S shown Fig. 8. Interestingly, for a given value of ϕ S , Eq. (93) gives now a single solution for |S|. Consequently, unlike their counterparts in Fig. 7, the contours no longer form closed loops, thereby indicating that the degeneracy with respect to ϕ S has disappeared. In Fig. 9, we illustrate this strategy, which is actually more general, i.e. does not only apply to the case of x = 0.

Further experimental data
|S|, ϕ S Sign information Figure 9: Flowchart to illustrate the use of the relations in Subsection 4.1 with information on C 10 to convert the measured value of R into predictions of the B 0 s → µ + µ − observables. Once these are measured in accordance with the pattern characterizing the NP scenario, |S| and ϕ S can be extracted from the data.
In the case of C µµ , the symmetry is broken by an overall minus sign: More explicitly, we have As we will see below, this feature has interesting phenomenological implications.
In order to illustrate the expressions given above, we consider two examples with different values of the coefficient C 10 :

Example (a):
We first assume a situation with vanishing NP contributions C NP 10 = 0, and employ the following setup: Using Eq. (93), we determine |S| as a function of ϕ S . As discussed above, for |C 10 The corresponding values for the observables read as follows: Let us now assume that these observables have been measured, and discuss how we may then -with the help of the strategy discussed above -reveal the dynamics of the B 0 s → µ + µ − decay and distinguish between the x = 0 and |x| → ∞ cases: • It is plausible to expect that A µµ ∆Γs is the next observable to be measured. With the help of the top-right plot in Fig. 7 • This ambiguity can be resolved through information on the sign of C µµ , which is given by C µµ = −0.16 and C µµ = +0.16 for |x| → ∞ and x = 0, respectively, as can be seen in Fig. 7. Consequently, the fact that C µµ breaks the symmetry in Eq. (97) gives us a powerful tool to distinguish between x = 0 and |x| → ∞.

Example (b):
Now we have a look at a scenario with NP contributions to C 10 , which is characterized as follows: Here the value of C 10 follows from Eq. (76), and is discussed in more detail in Subsection 4.2. In contrast to Example (a), we obtain now a single solution for |S| from Eq. (93), which is given by Using Eq. (66), we find |P | = 0.81, ϕ P = 33 • , resulting in the following values of the observables: In analogy to Example (a), using the plots in Fig. 8, we may again show the compatibility of the "measured" observables with the scenario x = 0, and rule out the case of |x| → ∞ through the sign of the C µµ asymmetry. For the convenience of the reader, we summarize the main features of these examples in Table 1.
In Fig. 10, we show the correlation between A µµ ∆Γs and S µµ through the CP-violating phase ϕ S . It should be noted that the corresponding regions for |C 10 | = 0.84, ϕ 10 = 0 • and |C 10 | = 1, ϕ 10 = 0 • do not differ substantially and are included in a single plot. Due to the symmetry transformation in Eq. (97), the scenarios x = 0 and |x| → ∞ cover the same region once we make a scan over the full range of ϕ S . The allowed region in  Table 1: Summary of the strategy followed in Examples (a) and (b) to disentangle the scenario x = 0 from |x| → ∞ and determine the value of ϕ S .
1. The currently available measurement of R implies a remarkably constrained circular region in the A µµ ∆Γs -S µµ plane for CP-violating NP scenarios characterized by x = 0 and |x| → ∞.

A future measurement of the observable combination A µµ
∆Γs and S µµ lying outside the allowed region would rule out the x = 0 and |x| → ∞ scenarios.

The allowed region in the A µµ
∆Γs -S µµ plane is close to the unit circle. Consequently, due to Eq. (1), the observable C µµ is constrained to take a smallish value.
4. The allowed region is similar to the one arising for the scenario described in Section 4.2. While here ϕ 10 = 0 • would imply the SM results A µµ ∆Γs = 1 and S µµ = 0, in the case of x = 0 or |x| → ∞ we may still deviate substantially from the SM even in spite of having a vanishing phase ϕ 10 .
In a complementary way, if we can obtain the value of the phase ϕ S from external information or theoretical considerations, we will be able to predict the observables A µµ ∆Γs and S µµ compatible with vanishing short distance contributions C P,S or C P,S . Strong deviations from these determinations will indicate that the corresponding scenarios are not realized in Nature. A discussion of NP scenarios characterized by the relations P ± S = 1 (see Eqs. (88) and (96)) can be found in Ref. [6].

∆ = 0 •
Another interesting case arises if C S and C S have the same CP-violating phases, i.e. ∆ = 0 • , which yields leading to while the symmetry is again broken by the observable C µµ through an overall sign change: The three observables r, A µµ ∆Γs and S µµ in Eqs. (112)-(114) depend on the three unknowns x, |S| and ϕ S . Consequently, if the observables are measured, we may determine these parameters. The twofold ambiguity following from the symmetry transformation in Eq. (116) can be resolved through the measurement of the sign of C µµ . Unfortunately, in view of the highly non-linear structure of the equations, we cannot give simple analytic solutions. However, the parameters can be determined numerically. In Section 5, we will illustrate this determination through fits to scenarios of future measurements.
Alternatively, we can apply the strategy depicted in the flowchart in Fig. 9. We start with the experimental value of R given in Eq. (28). Furthermore, we assume that |x| = 0.5 and |C 10 | = 1, ϕ 10 = 0 • . This allows us to solve for |S| as a function of ϕ S , and to subsequently determine A µµ ∆Γs , S µµ and C µµ as functions of ϕ S . The results are shown as the blue contours in Fig. 11. Here, also the symmetric situation with |x| = 2 is shown in red, illustrating nicely how C µµ breaks the symmetry.
Finally, in Fig. 12, we show the correlation between A µµ ∆Γs and S µµ for |x| = 0.5 and |x| = 3. Contrary to the situation for x = 0, |x| → ∞, we are not constrained to a contour close to the unit circle, but can also obtain values in the interior region. For the scenario |x| = 3, the relations |S|, A µµ ∆Γs , S µµ and C µµ as functions of ϕ S are similar to the ones shown in Fig. 11.   Table 2: Summary of the scenarios described in Subsection 4.4. In all the cases we have assumed ϕ 10 = 0 • .
In the expression for r given in Eq. (112), a pole seems to arise for |x| = 1, which corresponds to However, this is a spurious divergence, which is cancelled by the C S − C S term in the expression for S in Eq. (6), implying Using the relations in Eqs. (63) and (64), we obtain Consequently, Eq. (5) yields and shows that also the divergence in Eq. (66) for |x| = 1, ∆ = 0 • is spurious. If we neglect, for simplicity, again the tiny NP contribution φ NP s to the B 0 s -B 0 s mixing phase, we obtain For a discussion of NP models describing this situation, see Ref. [6]. Obviously, also extensions of the SM with scalars, which couple in a left-right-symmetric way to quarks (see the operators in Eq. (3) and the relations in Eq. (119)), fall into this category. As in the case given by Eq. (79), the correlation between the observables A µµ ∆Γs and S µµ describes a circle with radius one. The overall phase ϕ P includes effects from the, in general, complex quantities C 10 and C S . This is particularly interesting if future measurements reveal (A µµ ∆Γs ) 2 + (S µµ ) 2 compatible with the unit circle and if we have bounds available on the phase ϕ 10 from other processes. Then, results incompatible with Eq. (79) will indicate the potential presence of a scalar or pseudo-scalar contribution.

∆ = 180 •
In the case of ∆ = 180 • , we obtain the following expressions for the B 0 s → µ + µ − observables: (130) These equations could be solved numerically to determine |x|, |S| and ϕ S , in analogy to the discussion of ∆ = 0 • .
It is interesting to have a closer look at x = −1, i.e. |x| = 1 and ∆ = 180 • . In terms of the short-distance coefficients, this case corresponds to Using the relations in Eqs. (63) and (64), we obtain furthermore implying Using Eqs. (26) and (27), we obtain Special care should be paid when using Eq. (134), since the expression on the right-hand side has to be greater than or equal to zero. This feature implies the following upper bound: where we have used that 1 + y s cos(2ϕ 10 ) ≥ 1 − y s , with y s given in Eq. (12). With the current experimental value of R in Eq. (28), the corresponding bound is given by The different scenarios described in the previous subsections are presented in Table 2, where we show the connection between the standard parametrization used for the short distance contributions and the SMEFT one introduced in Sec. 4.1.

Experimental Aspects
Up to now we have not considered experimental uncertainties in the observables A µµ ∆Γs , S µµ and C µµ when studying the different scenarios. Nevertheless, we would like to demonstrate the potential for the determination of the underlying parameters at future experiments. Since the asymmetries are not independent, due to the relation in Eq. (1), it is not possible to determine all four parameters |S|, ϕ S , |x| and ∆. However, as discussed in Subsection 4.3, we expect to have a better picture of physics beyond the SM by the time the CP asymmetries of B 0 s → µ + µ − have been measured. Therefore, we consider some of the examples discussed in Subsection 4.4, which correspond to specific values of |x| or ∆. We assume uncertainties for the observables, allowing us to extract the NP parameters through fits.
Unless specified otherwise, within this section we use a future measurement of where we have assumed a relative uncertainty of 10% for B(B s → µ + µ − ), which is achievable at the LHCb upgrade [33], while keeping the current central value fixed.
Notice that the relative uncertainty in our "measurement" forR in Eq. (137) leads to a 2σ tension with the SM. Thus the statistical significance will not be high enough to claim for the discovery of NP effects. The major limiting factor of the precision is the ratio f d /f s of the fragmentation functions of the B 0 d and B 0 s mesons [34], which is required for normalization purposes. To the best of our knowledge, no information about the expected precision of future measurements of A µµ ∆Γs , S µµ and C µµ is available. The key question we want to address is the precision of the measurement of these observables that is required to establish in particular new (pseudo)-scalar contributions at the 5 σ confidence level.

x = 0 and |x| → ∞
To begin with, we evaluate the impact of experimental errors for the observables in Example (a) of Subsection 4.4.1, corresponding to a scenario where x = 0. An absolute uncertainty of ±0.2 for the asymmetries leads to the following set of observables: A µµ ∆Γs = 0.58 ± 0.20, S µµ = −0.80 ± 0.20, C µµ = 0.16 ± 0.20.
In such a situation, S µµ would indicate CP-violating NP effects at the 4σ level, while A µµ ∆Γs and C µµ would deviate from the SM picture at the 2σ and 1σ levels, respectively. Let us assume that the values above have been measured at a future experiment, and that there are strong reasons to consider models characterized by x = 0. We will now illustrate through a χ 2 fit how well we can reveal the underlying decay dynamics.
Let us first obtain the regions allowed for |S| and ϕ S if we only include R and A µµ ∆Γs in the statistical analysis. Thus, using the expression in Eqs. (89) and (90) with the "data" in Eqs. (137) and (138), we perform a χ 2 fit to these two observables and obtain the blue contours in the left panel of Fig. 13, which correspond to 1σ allowed regions. We indicate the input parameters used to determine our observables in Eq. (138) with the green dot. This plot allows us to establish a non-zero value for |S| at the 3σ level. If we include also the "measurement" for S µµ indicated in Eq. (138), along with Eq. (91), and repeat the χ 2 fit, we can eliminate the dashed contour in the left panel and obtain the right plot of Fig. 13.  Figure 13: Illustration of the determination of |S| and ϕ S from the observables in Eq. (138) for a scenario with x = 0, which is degenerate with |x| → ∞. The contours correspond to the 1σ allowed regions obtained from χ 2 fits. In the left panel, we show the result of the fit to only R and A µµ ∆Γs , while in the right panel we have also included S µµ . The blue contours were obtained by assuming x = 0, whereas the red contours follow for |x| → ∞. A measurement of the sign of C µµ would allow us to distinguish these cases, excluding the |x| → ∞ scenario.
As we have pointed out in Subsection 4.4.1, there is a symmetry between x = 0 and |x| → ∞, implying the same values of A µµ ∆Γs and S µµ for these two cases. Conversely, we could not distinguish x = 0 and |x| → ∞ at the phenomenological level having only measurements of these observables available. Indeed, repeating the χ 2 fits assuming |x| → ∞ for the same set of input observables leads to the red contours in Fig. 13. Although we use x = 0 as our favoured model in this illustration, it would certainly be desirable to rule out the degenerate |x| → ∞ scenario. As C µµ breaks the symmetry by an overall minus sign, we could actually exclude the |x| → ∞ case through experimental information on the sign on this CP asymmetry. If we add C µµ to the analysis, a solution only arises in case of the x = 0 scenario, thereby singling out the blue contour. We would then find |S| = 0.43 +0.07 −0.08 , where C µµ has a minor impact on the numerical values themselves, apart from excluding |x| → ∞. In this scenario, the assumed experimental uncertainties in Eq. (138) would allow us to establish non-zero values of |S| and ϕ S at the 5σ and 7σ levels, respectively, which would provide highly non-trivial insights into the underlying dynamics.

∆ = 0 •
Let us now have a closer look at another interesting scenario: ∆ = 0 • , where C S and C S have the same CP-violating phases. The expressions in Eqs. (112)-(115) form a system of three independent equations which allows us to determine |S|, ϕ S and |x|. Due to the highly non-linear structure of the mathematical expressions, we cannot provide analytical solutions in general. Instead we give an example of how to solve the system through a χ 2 fit. We consider the input parameters |x| = 0.5, ϕ S = 20 • , |C 10 | = 1, where we have considered the same absolute uncertainties as in Subsection 5.1. Assuming that these observables have been measured correspondingly at a future experiment, A µµ ∆Γs would indicate NP at the 6σ level, while S µµ and C µµ would differ from the SM at the 2σ and 4σ levels, respectively. The latter observable would require a non-vanishing scalar contribution S. Performing a χ 2 fit to these quantities, we can determine the underlying decay parameters |x|, |S| and ϕ S simultaneously from the best fit point. We start our statistical analysis by considering only R, A µµ ∆Γs and S µµ . In the left and right panels of Fig. 14, we show the corresponding 1 σ confidence regions in the ϕ S -|S| and ϕ S -|x| planes, respectively. We obtain two solutions, given by the blue and red contours, as we expect based on the symmetry relations in Eq. (117). Our input parameters are indicated by the green dot. Consequently, non-zero values of |S| and |x| at the 4σ and 6σ levels, respectively, could then be established.
If we include also C µµ in the analysis, we can eliminate the solution corresponding to the red contours, since C µµ breaks the symmetry relation in Eq. (116) by an overall minus sign. The resulting 1σ regions are shown in Fig. 15, corresponding to the results As a matter of fact, non-zero values of these parameters could be pinned down at the 5σ, 4σ and 7σ levels, respectively. In general, the precision for the CP asymmetries required to determine |S|, ϕ S and |x| with a given confidence level depends on the situation in parameter space. Moreover, we may end up with an ambiguity even after including C µµ in the χ 2 fit. Nevertheless, this example nicely complements the one in Subsection 5.1 and shows the potential of the CP asymmetries to determine the (pseudo)-scalar contributions, and even to discriminate between the corresponding primed and unprimed Wilson coefficients.

Conclusions and Outlook
The rare decay B 0 s → µ + µ − has been in the focus of particle physics for decades, offering one of the theoretically cleanest probes for physics beyond the SM, in particular for new (pseudo)-scalar contributions, which are still largely unconstrained. Finally, this channel could be observed by the CMS and LHCb collaborations and is now an experimentally well established process, exhibiting a branching ratio encoded in R in the ballpark of the SM. The observable A µµ ∆Γs , which is accessible thanks to the decay width difference ∆Γ s and requires an untagged -but time-dependent -analysis, will play an important role to shed light on possible NP contributions to B 0 s → µ + µ − . In general, these effects involve also CP-violating phases, which are usually neglected in theoretical analyses for simplicity.
In this paper, we have presented a comprehensive strategy for the future LHC upgrade(s), allowing us to reveal the presence of new sources of CP violation. The key role in this endeavour is played by the mixing-induced CP asymmetry S µµ , which requiresin contrast to A µµ ∆Γs -also tagging information for the experimental analysis. Another observable, C µµ , would become accessible if the helicity of the final-state muons could be determined; already sign information for this CP asymmetry would be very valuable information. These three observables do not depend on the decay constant f Bs and are not affected by theoretical uncertainties.
Interestingly, the interplay of R with A µµ ∆Γs and S µµ allows us to establish new (pseudo)-scalar contributions and new sources of CP violation. In general, we can only obtain constraints as we do not have sufficient independent observables to determine the short-distance coefficients |S|, |P | and their phases ϕ S , ϕ P . To obtain further insights, additional information is required. This could either be obtained by assuming specific NP models, or in a model-independent way through relations between the shortdistance coefficients C ( ) P , C ( ) S , which can be derived within the SMEFT approach. We have followed the latter avenue, discussing a variety of scenarios to illustrate how the corresponding parameters can be determined from the measured observables.
Since the pseudo-scalar coefficient P involves C 10 , we need information on this quantity. By the time precise measurements of the observables A µµ ∆Γs and S µµ are available, we expect to have a detailed picture of C 10 , following from analyses of semileptonic rare B → K ( * ) µ + µ − and B s → φµ + µ − decays. Current anomalies in the data for the former and B → K ( * ) e + e − decays indicate NP effects in C 10 , which we have also considered in our explorations. It will be important to utilize CP violation in the corresponding observables in the future.
To the best of our knowledge, experimental feasibility studies for the measurement of S µµ at the LHC upgrade(s) are not yet available. Performing fits to the observables for given future scenarios, we find that an absolute precision at the 0.2 level for A µµ ∆Γs and S µµ could have a dramatic impact on our search for new (pseudo)-scalar contributions in leptonic rare B s decays, allowing us to reveal the underlying dynamics. We urge the LHC collaborations to add studies of CP violation in rare B 0 s → + − decays to their physics agenda for the long-term future and super-high-precision era of B physics.
In view of the complexity of the resulting general expressions, we refrain from listing them for the observables r, A µµ ∆Γs , S µµ and C µµ . However, we have given formulae for specific examples discussed in Subsection 4.4.