On the effective lifetime of $B_s \to \mu \mu \gamma$

We consider the $B^0_s \to \mu^{+} \mu^{-} \gamma$ effective lifetime, and the related CP-phase sensitive quantity $A_{\Delta\Gamma_s}^{\mu \mu \gamma}$, as a way to obtain qualitatively new insights on the current $B$-decay discrepancies. Through a fit comparing pre- to post-Moriond-2021 data we identify a few theory benchmark scenarios addressing these discrepancies, and featuring large CP violation in addition. We then explore the possibility of telling apart these scenarios with $A_{\Delta\Gamma_s}^{\mu \mu \gamma}$, once resonance-modeling and form-factor uncertainties are taken into account. We do so in both regions of low and high invariant di-lepton mass-squared $q^2$. For low $q^2$, we show how to shape the integration range in order to reduce the impact of the $\phi$-resonance modelling on the $A_{\Delta\Gamma_s}^{\mu \mu \gamma}$ prediction. For high $q^2$, we find that the corresponding pollution from broad-charmonium resonances has a surprisingly small effect on $A_{\Delta\Gamma_s}^{\mu \mu \gamma}$. This is due to a number of cancellations, that can be traced back to the complete dominance of semi-leptonic operator contributions for high $q^2$ -- at variance with low $q^2$ -- and to $A_{\Delta\Gamma_s}^{\mu \mu \gamma}$ behaving like a ratio-of-amplitudes observable. Our study suggests that $A_{\Delta\Gamma_s}^{\mu \mu \gamma}$ is -- especially at high $q^2$ -- a potentially valuable probe of short-distance CP-violating effects in the very same Wilson coefficients that are associated to current $b \to s$ discrepancies. Its discriminating power, however, relies on progress in form-factor uncertainties. Interestingly, high $q^2$ is the region where $B^0_s \to \mu^{+} \mu^{-} \gamma$ is already being accessed experimentally, and the region where form factors are more accessible through non-perturbative QCD methods.


Introduction
The decay B 0 s → µµγ has recently attracted new attention in connection with, but also independently from, the discrepancies currently observed in b → s transitions [1][2][3][4][5].One reason is the fact that, once sizeable new-physics effects on B 0 s → µµ have been excluded by recent measurements of its branching fraction, the decay B 0 s → µµγ becomes interesting, as it allows to probe a larger set of Wilson coefficients than B 0 s → µµ, in particular all those of current interest in connection with the mentioned discrepancies [6][7][8][9][10][11][12][13][14].Moreover, B 0 s → µµγ enjoys an enhancement due to the lifting of the chiral suppression by the photon coupling, which translates into a branching ratio of the order of 10 −8 [15][16][17][18].
Admittedly, from a theory point of view B 0 s → µµγ is not nearly as clean as B 0 s → µµ [19][20][21] because of the required B 0 s → γ form factors (f.f.'s), and the limited knowledge thereof [16,18,[22][23][24][25][26].However, this difficulty can be circumvented in selected regions of the measurable phase space [27], and/or by defining ratio observables whereby the f.f.uncertainties cancel to a large extent [17,18].Besides, from an experimental point of view the high-q2 B 0 s → µµγ spectrum may be accessed from the very same, abundant dataset as B 0 s → µµ [27], i.e. without direct detection of the photon.In fact, a search with this method has been performed recently by the LHCb experiment [28,29], with the full statistics collected so far, for dimuon masses above 4.9 GeV, yielding the first experimental limit on the branching fraction of this decay of 2.0 × 10 −9 at 95% CL for the mass region considered.This confirms that prospects of detection or of improved limits at LHCb are favourable with future datasets.
One under-explored feature of B 0 s → µµγ is its capacity to probe new CP violation, i.e. complex Wilson coefficients.Actually, and as well-known [30], in most constraining b → ssector branching ratios and CP-averaged angular observables, only NP contributions aligned in phase with the SM can interfere with the SM contributions.As a consequence, NP with non-standard CP violation is in fact constrained more weakly than NP where CP violation stems only from the CKM phase.As a further consequence, for the coefficients present in the SM, i.e.C 7 , C 9 and C 10 , the constraints on the imaginary part of the NP contributions end up being looser than on the real part, as we will discuss later.
The effective lifetime of B 0 s → µµγ allows to access the quantity known as A µµγ ∆Γs , which offers a sensitive probe of new CP-violating effects in b → s-sector Wilson coefficients.We introduce and calculate explicitly this quantity. 1A µµγ ∆Γs is naturally a ratio-of-amplitudes observable, so that sensitivity to hadronic uncertainties may be accordingly reduced.We explore this question in detail.Specifically, this quantity may be studied in two separate regions, to be indicated as low-and high-q 2 , located beneath and respectively above the J/ψ and ψ(2S) resonances.In the low-q 2 region, a calculation of the long-distance dynamics based on factorisation methods was recently made available [22].Besides, for both low and high q 2 , a further recent parameterization based on light-cone sum rules (LCSRs) has recently appeared in Ref. [33]. 2 For low q 2 we are thus able to perform an explicit comparison between the two parameterizations, and thereby explore the sensitivity of A µµγ ∆Γs to the f.f.choice.Although a similar comparison is not yet possible for high q 2 , we are however able to address the question of the A µµγ ∆Γs sensitivity to broad-charmonium resonances, which we find to be reassuringly small.This region is the most sensitive to the operators O 9 and O 10 allegedly responsible for b → s discrepancies, and is the region where B 0 s → µ + µ − γ may be accessible in the short/medium term with the method in Ref. [27].
The plan of the paper is as follows.In Sec. 2 we collect the basic facts about effective lifetimes and rederive the related CP-violating observable for the case of B 0 s → µµγ, as an example of a decay with more than two particles in the final state.We specifically discuss the necessity of including phase-space integrals, separately in the numerator and denominator of the A µµγ ∆Γs definition.In Sec. 3 we discuss to what extent b → s and related data as of Moriond 2021 allow for complex contributions to the Wilson coefficients C 7,9,10 .We identify several benchmark new-physics scenarios.Our aim is to compare their effect on A µµγ ∆Γs with the effect of hadronic uncertainties, including f.f.parameterisations, as well as resonant effects from the φ and from broad charmonium.We perform such comparison in Sec. 4, devoted to low q 2 , and Sec. 5, on high q 2 .Finally, in Sec.6 we provide a discussion of the experimental prospects as well as our conclusions.We relegate to appendices: notational details on the necessary amplitudes and f.f.'s (Appendix A); a discussion of f.f.'s within the BBW approach [22] and a comparison with the LCSR parameterisation in Ref. [33] (Appendix B); explicit A µµγ ∆Γs formulae (Appendix C); a detailed analytic argument on why the uncertainty induced by broad-charmonium resonances is small (Appendix D); our input table (Appendix E).

Basics and main formulae
Given a final state f common to both the B 0 s and the B0 s , the most 'natural' experimental observable assuming equal production rates for a B 0 s and a B0 s is the 'untagged' rate [32,36] where with dΦ f an element of n-body phase space for the final state f [37], and where Af (t) denote the amplitudes of decay to f for a B s -meson that was a B 0 s at t = 0.The explicit time dependence of the two amplitudes squared appearing on the r.h.s. of eq. ( 1) is well-known.One introduces Af (t)| 2 can be calculated as follows [36] where we omitted terms of O(a) in the | Āf (t)| 2 case.Clearly, the sin and cos terms in eq. ( 4) will cancel in the sum on the r.h.s. of eq. ( 1), and one gets [32] The A f ∆Γs expression usually found in the literature is the r.h.s. of eq. ( 6), without the phase-space integrals.We note explicitly that these phase-space integrals are separate in the numerator and the denominator of A f ∆Γs , see eq. ( 6).Therefore, for decays to more than two particles in the final state, they cannot be accounted as the overall normalization factor that is usually understood in the literature.In fact, the products Āf A * f , |A f | 2 and | Āf | 2 will depend on kinematic invariants if f consists of ≥ 3 particles, and this dependence is different for the different Wilson-coefficient combinations that these amplitude products depend on.The integrals allow to integrate out such kinematic dependence, in accord with the l.h.s. of eq. ( 5), which depends on t only.
From eq. ( 1) one may define the general relation between the experimental and the theoretical branching ratio as [32,36,38] where The parameter A f ∆Γs relating the two branching ratios is final-state as well as model dependent.Interestingly however, A f ∆Γs can be extracted directly from another observable, the effective lifetime [31,32], defined as In the above relations we introduced, as customary, the average B s -system lifetime τ s and the fractional lifetime difference y s ≡ ∆Γ s /(2Γ s ).

A µµγ ∆Γ s calculation
The calculation of A µµγ ∆Γs is a special case of eq. ( 6).The integral over the 3-body phase space reads where, after defining the decay kinematics and a hatted quantity denotes that it has been normalized to the appropriate power of M Bs to make it dimensionless.Further using We obtain where, to ease notation, we used the abbreviation Aµµγ , (14) took into account that | Ā| = |A|, and again neglected terms of O(a) (see eq. ( 3)).An explicit formula for A and the relation between A and Ā are reported in Appendix A. The corresponding formula for A µµγ ∆Γs is provided in App. C. A few comments are in order.First, the CKM phase in the weak-Hamiltonian coupling in ĀA * exactly cancels the q/p phase, as expected.More specifically, after introducing the CP transformation with φ CP a convention-dependent phase (φ CP = 0 in FLAG [39] and = π in [36]), the phase φ M appearing in the q/p ratio (see eq. ( 3)) is given by (see e.g.[40]) Making the flavour of the initial state explicit, we then have Āf f , (q = d, s), a convention-independent quantity that depends on the initial state B 0 q as well as on the final state f .The quantity ξ (q) f can be determined in terms of two of the three observables where, for q = s, the first relation coincides with eq. ( 6).Since the φ M dependence is the same in Āf /A f and in Āf A * f , the latter appearing in eq. ( 13), we see that A f ∆Γs is phase-convention independent, as well-known [36,40].
We observe that, given the very A µµγ ∆Γs definition, complex phases in any of the Wilsoncoefficient combinations contribute to 'misaligning' numerator and denominator in eq. ( 13), in turn causing A µµγ ∆Γs to depart from unity.3This makes A µµγ ∆Γs a very sensitive probe of new, short-distance, CP-violating effects; besides, the fact that A µµγ ∆Γs is a 'ratio-of-amplitudessquared' quantity helps reducing its sensitivity to certain hadronic effects for high q 2 , as we will see in detail in Sec.5.2.1 and Appendix D. Of course, the above mentioned 'misalignment' depends on the theory scenario assumed, e.g. the SM vs. a given new-physics shift to the Wilson coefficients, in particular on the imaginary component of such shifts.Hence, to address the hadronic sensitivity of A µµγ ∆Γs in detail, we have to first establish theory scenarios to use as benchmarks.We discuss the latter in the next section.

Parenthesis: a CPV Global Fit in the Light of Recent Data
This section lies somewhat outside the main line of discussion of this paper, which is centered around A µµγ ∆Γs .This section, however, emphasises that large CP-violating effects on certain b → s Wilson coefficients [14,30,41] are still possible with updated data, and that A µµγ ∆Γs offers a new, theoretically clean observable to put them to the test.
To make our point, we consider the b → sµµ effective Hamiltonian (for notation see Appendix A), and we compare real vs.complex deviations on one-Wilson-coefficient combinations including C .Blissfully, this is still the case after the most recent LHCb analyses including the full Run-2 dataset [28,29,42].These updates, recently presented at Moriond EW 2021, are included in our analysis, which is summarised in table 1. Existing studies discussing the complex-Wilsoncoefficient case include [14,30], as well as the very recent [41] with a maximum-likelihood approach, as implemented in the Python packages smelli and flavio [43][44][45].We scan the likelihood of each scenario in the Re(C NP i )/ Im(C NP i ) plane, using all observables directly relevant to the b → s sector, which are summarised in table 2. In particular, we include the R K and B 0 s → µ + µ − updates as of Moriond 2021. 4  The important qualitative message from table 1 is that, for a given Wilson coefficient or Wilson-coefficient combination, sizeable imaginary contributions are entirely compatible with data.Representative shifts can be read off from the second column in table 1, which refers to the NP contribution only.These shifts show that the wealth of available b → s data is still relatively under-constraining for Wilson-coefficient shifts with a large phase not aligned to the corresponding one in the SM contribution.As we discuss in later sections, A µµγ ∆Γs is a novel, efficient probe of precisely this case.
In fig. 1 we also show the 1, 2, 3σ CL contours for the global fit in four selected scenarios, as solid black lines.These scenarios are C NP 7 , C NP 9 , C NP 10 , and 5  We also plot individual constraints coming from the subsets of observables that are the most constraining among all those included in both smelli and flavio.These subsets are indicated in boldface in table 2. We make the following comments about these results: 4 The recent B 0 s → µ + µ − update includes the first (ever) limit on B 0 s → µ + µ − γ in a limited portion of the q 2 spectrum, following the method in Ref. [27].This decay is sensitive to all Wilson coefficients relevant to the current discrepancies [23,46].The new bound is not yet constraining for our analysis. 5These scenarios will henceforth be referred to by dropping the NP superscript, i.e. as C7, C9, C10 and CLL, respectively.In addition, the 'SM scenario' will denote the case with all C NP i = 0.
• The fit improvement with respect to the purely SM solution (Re(C NP i ) = Im(C NP i ) = 0) is quite significant for the real part of the Wilson coefficient in the C 9 , C 10 and C LL scenarios, as quantified by the pull, and as well-known [6][7][8][9][10][11][12].On the other hand, in each considered scenario there is no significant pull on the sheer imaginary part of the NP Wilson coefficient considered: in each scenario, the SM solution is consistent with the available data.More CP-sensitive data would therefore be welcome in this respect.
• The allowed region in the C 9 and C LL scenarios is smaller in our case as compared to the results of [14].The difference can be traced back to both updated measurements and the addition of new observables to the global fit.In particular, the recent measurements of the branching ratios and angular observables in B → K * µµ decays increases the constraining power along the imaginary axis; the addition of Λ b → Λµµ branching ratio and angular observables pulls the allowed region towards the SM prediction; the updates in R K/K * and in the B s → µµ branching ratio modifies the constraint on the real axis.
• The imaginary part of C 10 and C LL is mostly constrained by the updated measurements in B → K ( * ) µµ observables.In particular, the addition of B 0 → K * µµ CP asymmetric angular observables [4] (yielding the contours displayed in darker blue) shows a preference for Im(C NP 10 ) > 0 and Im(C NP LL ) < 0. On the other hand, R K ( * ) and B(B s → µµ) constrain more the real part.
• C 7 is mostly constrained by B → K * µµ and b → sγ observables. 6We do not find a significant deviation from the SM prediction for either the real or imaginary parts.Among the b → sγ observables implemented in flavio and smelli, the most constraining ones are the branching ratios of B 0,+ → K * 0,+ γ, B s → φγ and B → X s γ.On the contrary, A ∆Γ (B 0 d,s → φγ), S K * γ and S φγ do not yield competitive constraints on the real and imaginary parts of C 7 .
From the above fits, we infer the following benchmark scenarios to be considered in the sections to follow: These benchmark points are chosen to lie within the 1σ region of the fit, and to have an imaginary part as large as possible.

A µµγ ∆Γ s at low q 2
The scenarios collected in eq. ( 18) serve us to establish reference shifts to A µµγ ∆Γs with respect to its SM value.We can then investigate to what extent these shifts are resolved once we keep into account hadronic uncertainties due to f.f.'s and to resonance modelling.
We first consider the region of low q 2 , located between the lower kinematic limit 4m 2 µ and an upper bound around 6 GeV 2 , determined in order to minimise possible pollution from the  J/ψ-resonance sideband.The low-q 2 region is of interest because of the sizeable enhancement of the C 7 contribution due to the 'nearby' photon pole, and because a determination of the f.f.'s using factorisation methods has recently become available [22].Away from resonant regions, this determination allows for a rigorous assessment of theory errors.Besides the analysis in Ref. [22] (BBW in what follows), another Bs → γ f.f.determination, based on Light-Cone Sum Rules (LCSR) has recently appeared in Ref. [33], and will henceforth be denoted as JPZ.The JPZ parameterization also holds for the high-q 2 region to be discussed in Sec. 5. Before the BBW and JPZ parameterizations, the only other existing f.f.parameterisation, Ref. [18], was based on a phenomenological model.
We are thus in a position to perform a direct comparison between the BBW and JPZ f.f.'s in A µµγ ∆Γs .To this end, we first re-express the full amplitude calculated in [22] in terms of the F V,A,T V,T A f.f.'s used in [18] as well as in related literature [24,35].We thereby obtain 'effective' BBW f.f.'s, F BBW V,A,T V,T A , in the sense that the amplitudes in eqs.( 29)- (30), with the effective f.f.'s F BBW V,A,T V,T A , reproduce the BBW amplitude.A detailed discussion, along with explicit formulae, and a numerical comparison of F BBW V,A,T V,T A with the JPZ f.f.'s [33], is presented in Appendix B.
Here we would like to study the question how the f.f.choice, BBW vs. JPZ, impacts the prediction on A µµγ ∆Γs in the different scenarios of eq. ( 18) as well as in the SM.Before presenting such comparison, an important qualification is in order.The low-q 2 region we are considering includes a narrow resonance, the φ, which escapes at present a description in terms of fully rigorous QCD methods, and on the other hand accounts for a substantial fraction of the low-q 2 signal.One has therefore to establish a compromise between minimising the theoretical error and maximising the experimental statistics.In our numerical study, we use sliding s 1 and s 2 values such that 4m 2 µ < s 1 < M 2 φ < s 2 < 6 GeV 2 , and we identify the GeV 2 ] with the above compromise in mind.The φ region is at present accounted through phenomenological approaches, whereby the relevant amplitude is shifted by a Breit-Wigner(BW)-like shape, and the latter is suitably parameterised.We follow the approach in Ref. [17], that we very briefly summarise here.One first identifies the relevant reduced amplitude where T ⊥, (q 2 ), G ⊥, (q 2 ) and L i (q 2 ) denote long-distance contributions, and takes into account diagrams where the e.m.-penguin photon emits the lepton pair and diagrams where the e.m.-penguin photon is the final-state one.Note that in more common notation T B0 s →φ 2 (0), see e.g.[69].One then rewrites this amplitude as an n-times subtracted dispersion relation, such description being necessary because of the resonant behaviour.Taking n = 1 yields This relation corresponds to the one also used in [23] in the approximation Â B0 [23] (see also [24,35]).The f.f. is known most precisely from LCSR,7 yielding [69-71] Finally, the O 8 and four-quark contributions in eq. ( 19) are known in the 1/m b -limit [72,73] and in LCSR [74,75].As pointed out in Ref. [17], an alternative and possibly more effective strategy towards improving the prediction for the amplitudes Â B0 s →φγ ⊥, appearing in eq. ( 21) is to extract them from experiment.This approach is promising since the branching ratio [37] is known to 10% accuracy.An update including the entire dataset is in progress.Hence the statistical component of the error on this measurement-about half of the error quoted in eq. ( 23)-will decrease steadily.As a consequence one can expect to extract the amplitudes at the 5% level, which is well below a theory error above 10%. 8In other words, rather than taking T B0 s →φ 1 (0) from eq. ( 22)-which translates into an error around 15% on the low-q 2 B 0 s → µ + µ − γ branching ratio-one may trade T B0 s →φ 1 (0) for eq.( 23). 9 This translates into which is over 1 standard deviation above eq.( 22).The above discussion allows to identify the central value in eq. ( 22) as reference, as well as the two values (T 1 ) min = 0.282 and (T 1 ) max = 0.393, which are respectively the 1σ lower bound from eq. ( 22) and the 1σ upper bound from eq. ( 24).These two values determine a realistic range for T , that in particular takes into account the 'tension' between the determinations ( 22) and (24).
A first task is thus to identify the in the [(T 1 ) min , (T 1 ) max ] range affects the A µµγ ∆Γs prediction negligibly with respect to the rest of the error components.We display such dependence in fig. 2 for the SM case as well as for all the NP scenarios in eq. ( 18).
The figure shows that the choice of T B0 s →φ 1 in the quite generous range discussed has little impact on the prediction of A µµγ ∆Γs even for s 1 as large as 0.035 and s 2 as small as 0.037-for reference, q 2 = M 2 φ corresponds to s = 0.036.This conclusion holds generally true in all considered theory scenarios.The figure suggests however to choose ŝ1 ≤ 0.03 and ŝ2 ≥ 0.05 Bs , see text below eq.( 10)), given the rapid variation of the A µµγ ∆Γs prediction above and respectively below such values.With these ŝ1 and ŝ2 values in mind, we choose as the central value in eq. ( 22) for definiteness.We next discuss the comparison between the BBW and the JPZ parameterisation on the A µµγ ∆Γs prediction in the low-q 2 region of integration.Of course, such comparison depends on the errors associated to either parameterisation.
The error on the BBW parameterisation is determined as the envelope of the errors discussed in detail in Appendix B, and in most of the integration interval is entirely dominated by the uncertainty on the inverse moment of the B s -meson distribution amplitude, λ Bs .The JPZ parameterization comes with an estimation of the associated errors, as well as of the correlations between the different form factors, and we adhere to these estimations.In the range ŝ ∈ [0, 1], f.f.errors are around 10% for F V and F T V , and comprised between 9% and 17% (28%) for F A (F T A ).
The A µµγ ∆Γs prediction integrated in the low- GeV 2 ] with BBW vs. JPZ f.f.'s is displayed in fig.3, the two columns of panels denoting the region below and respectively above the φ peak, and the rows referring to the different scenarios, including the SM and those in eq. ( 18).The figure shows that, whatever the choice of s 1 and s 2 , the two parameterisations yield generally consistent results, within the large errors.There are a few notable exceptions though.The clearest is the C 7 scenario in the full [4m 2 µ , ŝ1 ] 8 As discussed in Ref. [17], one may apply the above approach to other resonances, in particular the φ .In spite of a potentially sizeable T B 0 s →φ 1 coupling, the large φ width turns out to suppress the φ contribution to the B 0 s → µ + µ − γ spectrum to be a below-1% correction to the total branching ratio.Finally, the above approach may be applied to narrow charmonium as well.However, the required radiative branching ratios are, again, not yet measured.Besides, short-distance dynamics in this region is dominated by the O GeV 2 ] (right panels), within the scenarios indicated.We use JPZ f.f.'s [33] as well as three choices of T

B0
s →φ 1 -see below eq.(24).range.The JPZ parameterization tends to predict lower values for Re(F T A ) and especially Re(F T V ) than the BBW counterparts, see fig. 6 in Appendix B. This difference translates into lower A µµγ ∆Γs predictions in the full [4m 2 µ , ŝ1 ] range for any scenario (see left column of fig.3).These predictions become sizeably lower in the C 7 scenario, where the NP shift multiplies the F T V,T A f.f.'s.The other appreciable difference occurs in the C LL scenario for ŝ2 0.10 and in the C 9 scenario for ŝ2 0.07.We note in this respect that the JPZ shape visibly decreases for decreasing ŝ2 , whereas the BBW shape is nearly constant throughout the ŝ2 range.This difference holds true in the full ŝ2 interval considered, and for any scenario, hence it is not attributable to features specifically related to the φ-resonance region.On top of this difference, the BBW error in A µµγ ∆Γs tends to squeeze for decreasing ŝ2 in the SM, C 9 , and C LL scenarios, due to an accidental cancellation between λ Bs -and r LP -induced components of the uncertainty, namely the two dominant ones.The fact that, aside from these somewhat accidental instances, the two parameterizations give consistent results on A µµγ ∆Γs is, we believe, a non-trivial finding.

A µµγ
∆Γ s at high q 2 In this section we discuss the high-q 2 region, the one most sensitive to O 9,10 , and the region where B s → µµγ may be accessed in the short/medium term with the method in Ref. [27].
Throughout the numerics in this section, we use JPZ f.f.'s [33] only, as BBW f.f.'s are deemed valid only for q 2 values below the narrow-charmonium threshold [22].

NP sensitivity at high q 2
We would like to first display the sensitivity of A µµγ ∆Γs to the new-physics scenarios summarised in (18).To this end, we fix for definiteness the ŝ integration range as The lower bound, corresponding to √ s 4.1 GeV, has been chosen to minimise the possible contamination by broad-charmonium resonances, while maximising statistics.This value will be justified better in Sec.5.2.1.Conversely, we fix the upper bound to be the endpoint ŝ = 1.More generally, the upper bound may be expressed as ŝmax = 1 − 2E cut /M Bs , where E cut may be interpreted as the smallest photon energy detected by the experimental apparatus.This quantity is actually inferred from the apparatus' resolution in s, and is of the order of 50 MeV. 10Concretely, the region of the B s → µµγ spectrum very close to the endpoint (i.e. with 2E cut /M Bs 1) is completely dominated by bremsstrahlung (a.k.a.final-state radiation) contributions [77,78], that experimentally are routinely subtracted by a MonteCarlo [79].Throughout this paper, we accordingly consider the B s → µµγ spectrum with the bremsstrahlung contribution set to zero, whereas we keep interference terms, that the MonteCarlo does not subtract.
With these qualifications, a first interesting question is the A µµγ ∆Γs central-value prediction as a function of Re(C NP i ) vs. Im(C NP i ), with i = 7, 9, 10, LL, where we recall that the i = LL case denotes C NP 9 = −C NP 10 .We display such dependence in the panels of fig. 4. We also show as solid, dashed or dotted red contours the 1, 2 and 3σ regions allowed by the global fit described in Sec. 3. GeV 2 ] (right panels), within the scenarios indicated.We superimpose BBW [22] (yellow) and JPZ [33] (blue) f.f.'s.See text for further details.As expected, in the chosen integration region the dependence of A µµγ ∆Γs on C NP 7 is fairly weak, although shifts in Im(C NP  7 ) allowed by present global fits would lead to departures with respect to (A µµγ ∆Γs ) SM of O(10%), thereby in principle detectable.Effects in the C 9 , C 10 scenarios, and especially in the C LL one are more sizable, because in this region the dependence of A µµγ ∆Γs on these Wilson-coefficient combinations is stronger than it is on C 7 , similarly as for the B s → µµγ differential width.
The effects just discussed have to be compared with the theory error associated to A µµγ ∆Γs , which we do next.

A µµγ ∆Γ s theory error at high q 2
The two outstanding sources of theory error for A µµγ ∆Γs in this region are the possible pollution from broad-charmonium resonances and the f.f.uncertainty.We address these two sources in turn.
Following a standard approach to account for these effects, we note that they enter our decay of interest via the subprocess B 0 s → V (→ )γ, where V denotes any of the aforementioned resonances.The associated long-distance contributions may be modelled as a sum over Breit-Wigner (BW) poles [24], i.e.
In this relation, In order to address the theory uncertainty due to this modelling, we subsequently scan simultaneously over the BW normalisation factors |η V | ∈ [1, 3] as well as on the phases δ V ∈ [0, 2π) with independent uniform distributions. 11Such procedure is expected to provide a conservative way to measure the deviation from naive factorisation (|η V | = 1 and δ V = 0). 12We perform the above scan for ŝmin ∈ [0.50, 0.70], and for all the NP scenarios in eq. ( 18), as well as for the SM case.The result is shown in fig. 5 as yellow bands.For each given value of ŝmin the central value and the values in the upper and lower bands are calculated as, respectively, the mean and ±1σ of the distribution obtained through the corresponding scan.The panels in fig. 5 show that the broad-charmonium modelling uncertainty due to eq. ( 26) (yellow 'bands') has no appreciable impact on the prediction of A µµγ ∆Γs , compared with the uncertainty (blue bands) induced by f.f.'s, to be discussed in Sec.5.2.2.Since the shift induced by eq. ( 26) is typically of O(5%), and its phase is completely misaligned with the C 9SM phase, it is not obvious why most of this uncertainty cancels out in A µµγ ∆Γs .On closer inspection, one can understand this cancellation as the result of specific features, in particular the complete dominance, in this kinematic region, of quadratic contributions in C 9 and C 10 ; the fact that the multiplying form factors F V and F A are real; the fact that the broad-charmonium shift in eq. ( 26) can be treated as a 'small' modification of the numerically large SM value for C 9 .These features ease cancellations between the numerator and the denominator of A µµγ ∆Γs -suggesting that, with respect to this class of long-distance contributions, A µµγ ∆Γs behaves well like a ratio observable.We present a more detailed analytic argument in Appendix D.
Because of the above conclusions, the C 9 and C 10 dominance in this region may be regarded as one key advantage with respect to the low-q 2 region, where the C 7 contribution becomes important, without being dominant in the full kinematic range.It is quite encouraging that, as a result, A µµγ ∆Γs is an observable largely immune to broad-charmonium uncertainties, given that they escape a rigorous description.The control of the theory prediction rests thus entirely in the f.f.error, to which we turn next.

Impact of form-factor modelling
A further source of uncertainty that needs be addressed is the error attached to the f.f.'s F V,A,T V,T A that parameterise the B0 s → γ amplitudes, see eq. ( 31).As discussed in Sec. 4, the JPZ parameterization includes an estimate of the error on each of the F V,A,T V,T A f.f.'s [33]. 13We accordingly vary these f.f.'s within uncorrelated normal distributions around their respective errors.The impact of these errors on A µµγ ∆Γs is shown in fig.5, as blue bands.In calculating these bands, the narrow-charmonium contributions are modelled with eq. ( 26), with |η V | = 1 and the δ V phases as in table 3, although we verified that a simultaneous variation of broad-charmonium and form-factor parameters leads to tiny differences.Fig. 5 allows to visually compare the impact of the two sources of uncertainty, broad charmonium vs. f.f.'s, for each theory scenario.
We see that, in the whole ŝmin range considered and for any theory scenario, the f.f.uncertainty is always much larger than the uncertainty induced by broad-charmonium modelling.The basic difference between the two sources of error may be appreciated as follows.On the one hand, the broad-charmonium correction enters 'only' as a shift to C 9 , hence efficient cancellations can take place, as discussed in detail in Appendix D. On the other hand, f.f.'s enter in different ways, all numerically relevant, for the different Wilson-coefficients combinations, as already mentioned at the end of Sec.2.1.As a consequence, cancellations of f.f.effects are much less efficient.The A µµγ ∆Γs error induced by f.f.'s is still too important to allow to tell apart the chosen theory scenarios from each other.However, it is noteworthy that, in this kinematic region, the dominance of, jointly, C 9 and C 10 gives rise to a particularly high sensitivity to the C LL scenario, which could be resolvable with an improvement on f.f.'s uncertainties.

Experimental Outlook and Conclusions
In this paper, we have discussed the theoretical interest of pursuing effective-lifetime measurements for B 0 s → µ + µ − γ.With the method of Ref. [27] these measurements, in the high-q 2 part of the spectrum, can be performed on the very same event sample as for B 0 s → µ + µ − , i.e. di-muon events with no detected photon.The case of B 0 s → µ + µ − γ is actually one of many possible examples.In fact, the underlying idea-the idea that if one can estimate all other B 0 s → µ + µ − backgrounds as q 2 departs downwards from the M 2 Bs peak, one may extract B 0 s → µ + µ − γ 'parasitically' as one particular background-is in principle applicable to any event of the kind B 0 s → µ + µ − X.Such decay allows to reconstruct a primary and a secondary vertex (the latter through the di-muon), and thereby to perform time-dependent measurements, in turn usable to reconstruct the corresponding effective lifetime. 14he specific task of measuring the effective lifetime for B 0 s → µ + µ − γ introduces a number of daunting challenges-although not insurmountable.The measurement of the effective .The blue vs. yellow bands refer to the f.f.error, and respectively on the uncertainty associated to the modelling of broad-charmonium resonances.See text for details.lifetime of a B 0 s decay requires a clean signal sample, a task that for B 0 s → µ + µ − γ decays is non-trivial, and a well measured secondary vertex, which is typically the case for displaced di-muons at LHC experiments.At present, two possible, different approaches exist to reconstruct these decays: a full reconstruction, more suited for the low-q 2 part of the spectrum, where the photon energy is larger and, as such, easier to reconstruct; and a method relying on partial reconstruction [27], mostly relevant for the high-q 2 range, where the photon energy is small, making photon reconstruction inefficient.The two methods present different challenges as concerns possible backgrounds.
In the high-q 2 region, the partial-reconstruction method allows high efficiency, comparable to the B 0 s → µ + µ − decay for which, as mentioned, a first measurement of the effective lifetime has been performed by the LHCb [64] and CMS [66] experiments, combined in Ref. [67]. 15ssuming this method, the closer the di-muon mass to the B 0 s mass the more similar the experimental performances in B 0 s → µ + µ − γ will be to the B 0 s → µ + µ − case.Alas, the backgrounds are considerably different, due to the absence of a clear narrow peak on which to perform a statistical background subtraction from the invariant mass sidebands.In particular the semileptonic B → hµν, with the hadron h = π, K misidentified as a muon, constitute a large background, together with the rarer B → hµµ decays, where the additional hadron is not reconstructed.The latter in particular will constitute almost an irreducible background, although one possible handle to tame it may be the shape of the di-muon mass distribution.A measurement of the B 0 s → µ + µ − γ effective lifetime thus requires both a precise knowledge of the yield and of the di-muon mass distribution of the listed backgrounds.This is not unrealistic, for example a large sample of B 0 s → K − µ + ν µ decays has been recently seen at LHCb [84] and can be used to constrain its equivalent misidentified yield.The lifetime of this decay, being a flavour-specific decay, is precisely known.Similar backgrounds from B 0 and B + decays will also have precisely known lifetimes.In this partially reconstructed method, a small uncertainty can arise from the absence of the photon momentum in the full B reconstructed momentum used to calculate the boost (and hence the lifetime from the decay distance): the higher the q 2 considered the smaller is this uncertainty.This would add a small component to the experimental uncertainty on the lifetime on an event-by-event basis, which will average out with statistics.In addition, kinematic fits exploiting the known flight direction of the B can further improve this resolution and relative measurements can cancel out these effects, and indeed partially reconstructed semileptonic decays are used for precision lifetime measurements [85].
Conversely, at low di-muon mass, in order to distinguish B 0 s → µ + µ − γ signal decays from background it will be crucial to efficiently reconstruct the photon.While this reduces the overall signal efficiency, it allows to statistically separate it from background through the invariant mass peak at the B 0 s mass.In this region, performances can be benchmarked against the B 0 s → φγ decay, that can also be reconstructed in the muonic final state, and that is well studied in its kaonic final state.
Clearly, a quantitative study is needed to translate these considerations into a luminosity vs.A µµγ ∆Γs -accuracy plot.In particular, in either of the above two cases detailed studies of the backgrounds are required to assess the performance on this new observable.Such survey requires a full experimental simulation and is thus beyond the scope of this paper.But it is clear that the measurement of this observable will require high statistics and is thus feasible in the upgrade phase of the LHC experiments.
In synthesis, the present paper introduces the B 0 s → µ + µ − γ effective lifetime.This observable offers the possibility of probing all Wilson coefficients that are currently of interest in view of the B anomalies, and, in particular, of novel CP-violating effects.Interestingly, scenarios with new weak phases not aligned with the corresponding ones in the SM contributions are under-constrained.Besides, the relevant CP-sensitive quantity-A µµγ ∆Γs -encoded in the effective lifetime is per definition a ratio-of-amplitudes (squared) observable.One can therefore expect a partial cancellation of the hadronic uncertainties inherent in the B 0 s → µ + µ − γ amplitude.We find that this cancellation is surprisingly efficient at high q 2 for broad-charmonium-modeling uncertainties, whereas errors induced by form factors are still sizeable and represent the main limitation at present.However, the high-q 2 region is especially favorable for non-perturbative, first-principle approaches to the calculation of the concerned form factors, e.g.within lattice QCD.Also interestingly, this is the same region where B 0 s → µ + µ − γ is already being accessed from the B 0 s → µ + µ − dataset.Refs.[33,34].FD would like to thank Flavio Archilli, Matteo Rama and the whole LHCb B 0 s → µ + µ − group for useful discussions.We acknowledge Peter Stangl for help on various aspects of flavio.This work is supported by the ANR under contract n. 202650 (PRC 'GammaRare').

A. Amplitudes and form factors: basic notation
The dynamics of the B0 s → µµγ amplitude may be parameterised by the following b → s effective Hamiltonian [87][88][89] where λ i ≡ V * is V ib and V the CKM matrix.The operators relevant to our discussion read Bs 0 → µµγ) entering eq. ( 39) may be expressed as the sum of a 'direct-emission' (DE) and a bremsstrahlung component [18,23], which read The O ( ) 7,9,10 matrix elements are parameterised by f.f.'s of dilepton momentum transfer q 2 [18, 23,24] and f Bs is defined through The amplitudes in eqs.( 29)-( 30) can be generalized to include primed coefficients through the substitutions where the +(−) refers to terms proportional to F V,T V (F A,T A or f B 0 s ).Eq. ( 29) exactly reproduces the result in [23]; eq. ( 30) matches the result in [23] provided X f = +1.We note however that the PDG [37] and FLAG [39] imply X f = −1.This point is accounted for in Refs.[18,46].
The calculation of A µµγ ∆Γs involves also the amplitudes for an initial B 0 s .The B 0 s → γ hadronic matrix elements are related to the B0 s → γ ones through a CP transformation, hence they depend on the phase φ CP appearing in eq. ( 15).The amplitude for B 0 s → µµγ is related to the amplitude for B0 s → µµγ by replacing CKM matrix entries and Wilson coefficients with their complex-conjugates, and by the following substitutions where denote any of F A , F T A , and f Bs in eqs.( 31)- (32), whereas F (B 0 s ) V,A denote the f.f.'s and decay constant for B 0 s matrix elements.17As concerns the dependence on weak phases, we note that a B0 s (B 0 s ) initial state corresponds to a b → s ( b → s) Hamiltonian, with Wilson coefficients C 7,9,10 (C * 7,9,10 ).As a consequence, both of Ā and A * are proportional to C 7,9,10 and the numerator of A µµγ ∆Γs depends on C 2 i .Hence the dependence on the Wilson-coefficients' phases does not cancel in A µµγ ∆Γs as it would in, e.g., |A| 2 .However, and as remarked in the main text, the CKM-phase dependence in Āµµγ A * µµγ cancels exactly with the analogous dependence in q/p, so that the r.h.s. of eq. ( 40) is proportional to |N | 2 .As a consequence, A µµγ ∆Γs is sensitive to CP-violating phases coming from Wilson coefficients only.

B. BBW vs. JPZ form factors
For the hadronic matrix elements relevant to our study we use two alternative approaches.The first one relies on parameterising the required matrix elements as in eq. ( 31).This parameterisation [18,23,24] defines the f.f.'s F V,A,T V,T A .These f.f.'s are subsequently estimated using a relativistic quark model (see [90][91][92] for details).We use the latest determination in [18] as reference for this approach.Importantly, F T V (q 2 , 0) and F T A (q 2 , 0) denote only the contributions where the e.m. penguin emits the di-lepton pair.In order to also include topologies where the e.m. penguin emits the real, final-state photon, as well as so-called weak-annihilation contributions [93] one replaces in eqs.(31), with the FT V,T A functions defined in Ref. [18].Recently, Beneke, Bobeth and Wang [22] have performed the first evaluation of the Bs → µµγ amplitude at low q 2 < 6 GeV 2 using rigorous factorisation methods.The amplitude is expressed as the sum of initial-plus final-state-radiation terms.The final-state-radiation terms correspond to B0 s -to-vacuum matrix elements proportional to O 10 , hence helicitysuppressed and completely negligible [22,27,94].On the other hand, from the initial-stateradiation terms one can read off the 'equivalent' of the F V , F A , FT V and FT A functions, to be denoted as where i = 9, 10 and C i stands for C eff 9 (q 2 ) and C 10 [95,96], respectively, and In the above expressions: 19• For the , where C i = C eff 7 , C eff 9 , C 10 , with µ h the hard scale, we use expressions in [22].We take the necessary F j(u) i functions from [97], the F j(c) i functions from [98](see also [99,100]), with both sets validated through [86].
• The other main source of error comes from the next-to-leading-power f.f.'s ξ and ξ.Following [104], Ref. [22] includes in the ξ and ξ definition a term that subtracts the leading-power contribution in the tree approximation, multiplied by r LP = 0.2 ± 0.2.The r LP range yields the second largest source of error in our case.
• For the m b and m c masses, we use the pole-mass values in table 3. 21 The distinction with respect to MS masses is dropped in the next-to-leading-power terms in eqs.( 36)- • Any other function or symbol not commented upon can be retrieved from [22].
It is useful to explicitly compare the JPZ f.f.'s in the F V , F A , FT V , FT A nomenclature, with the BBW 'effective' f.f.'s in eqs.( 36)- (37) in the low-q 2 region.We do so in figs.6 and 7 for the real and imaginary parts, respectively.We show separately the BBW-f.f.range of variation due to the λ Bs and to the r LP uncertainties discussed in the previous paragraph.Note that the JPZ parameterization in Ref. [33] yields real f.f.'s.The only source of an imaginary contribution is the shift due to the φ resonance, which affects only the tensor f.f.'s.The 'dictionary', provided in Ref. [33], between JPZ and KMN f.f.'s allows to include such shift in a way akin to the KMN parameterization.We include this contribution only for the sake of the indicative comparison in fig.7, keeping in mind that the φ region is basically entirely excluded in the considered kinematic range.The figures show that the overall agreement between the JPZ and BBW f.f.'s is mostly acceptable within the quite large errors inherent especially to the BBW parameterization.

∆Γ s
Starting from eq. ( 13) we can perform the integration 1 −1 d cos θ.We obtain Figure 6: Comparison of the real parts of JPZ f.f.'s F V , F A (first row left vs. right), FT V (second row), FT A (third row), with the BBW 'effective' f.f.counterparts in eqs.( 36)- (37).Note that the first row shows vector and axial f.f.'s with an index i = 9, 10, as in eqs.(36).The BBW-f.f.range of variation due to the λ Bs and to the r LP uncertainties, both discussed in the text, are shown separately (see legend entries).For rendering reasons, in the instance of tensor f.f.'s we show separately the cases of ŝ below (left panels) and respectively above (right panels) the φ peak. and where the functions f ij are given by
We estimate the complex shift δ cc from a simultaneous scan to the ten resonance parameters |η V | ∈ [1, 3] and δ V ∈ [0, 2π) (see Sec. 5.2.1 for details).For the sake of the present discussion, we need an absolute upper bound on the size of this shift.To this end we choose ŝ = 0.49, which yields the 1σ-range |δ cc | ∈ (4.3 ± 1.6)% × C 9,SD , (56) with an arbitrary phase.In the discussion to follow, we will conservatively assume eq. ( 56) as the size of the δ cc complex correction throughout the integration range, ŝ ∈ [0.5, 1], considered for high q 2 .Such assumption will prove to be enough for the sake of our argument.We note, however, that this is a truly conservative assumption.In fact, fig.8  correction steadily decreases for ŝ > 0.49, and that the correction is below 1% for ŝ > 0.75, i.e. in half of the actual integration range.With this in mind, we will consider δ cc as constant over this range, and given, in size, by eq.(56) where, in view of eq. ( 56), we neglect terms quadratic in the δ cc shift, and we use the fact that f 99 (see eq. ( 44) and fig.7) is real.We note that the δ cc shift also affects C 7 × C 9 (eqs.( 48) and ( 49)) as well as C 9 × C 10 (eq.( 53)) interference terms, that are bundled in N rest and D rest .We will come back to these terms later on.Denoting Re δ cc / C9 = /2, and treating C9 and δ cc as constant over the q 2 range of interest (see footnote 22), we can rewrite eq. ( 57) as Because of eq. ( 56), we can expand eq. ( 58) for small , obtaining The N and D coefficients in the above expression can now be read off from relations ( 42)- (53).
In particular, for N 99 and D 99 we have Besides, we can express N rest and D rest as In the N rest and D rest relations we have not written explicitly terms that either involve C 7 , or are bilinear in C 9 C 10 , that are also concerned by the δ cc shift.These terms, indicated in eqs.( 62)-( 63) by δ N,D (δ cc ), are suppressed by the small size of the C 7 coefficient in comparison with |C 9,10 | (see footnote 22), or by the small ratio f Bs /M Bs 0.04.We find that the δ N,D (δ cc ) contributions are around 15% of those explicitly written.Hence, recalling eq. ( 56), the broad-charmonium shift within these contributions may be estimated as a ≈ 4% × 15% effect, i.e. well below 1%.

Figure 4 : 9 = −C NP 10 .
Figure 4: A µµγ ∆Γs contours in the Re(C NP i ) vs. Im(C NP i ) plane, for i = 7, 9, 10, and for the scenario Re vs. Im of C NP 9 = −C NP 10 .The solid or dashed or dotted contours denote respectively the 1, 2 and 3σ regions allowed by the global fit discussed in Sec. 3. See text for further details.

Figure 5 :
Figure 5: A µµγ ∆Γs prediction in the kinematic interval [ŝ min , 1].The blue vs. yellow bands refer to the f.f.error, and respectively on the uncertainty associated to the modelling of broad-charmonium resonances.See text for details.
with α, β color indices and the primed operators obtained from eq. (28) by the replacements {L → R, m b → m s }.We use the covariant-derivative definition D µ = ∂ µ + ieQ f A µ + ig s G µ (where e.g.Q e = −1).This definition, and the chosen normalisations of O 7,8 yield C SM 7,8 < 0. With the above definitions, the amplitudes

Table 1 :
Comparison table of real vs.complex one-Wilson-coefficient scenarios for the b → , C 10 or C LL are no better than no NP at all.