Collider constraints on light pseudoscalars

We investigate the bounds on light pseudoscalars that arise from a variety of collider searches. Special attention is thereby devoted to the mass regions $[3, 5] \, {\rm GeV}$ and $[9,11] \, {\rm GeV}$, in which a meaningful theoretical description has to include estimates of non-perturbative effects such as the mixing of the pseudoscalar with QCD bound states. A compendium of formulas that allows to deal with the relevant corrections is provided. It should prove useful for the interpretation of future LHC searches for light CP-odd spin-0 states.


Introduction
The most significant achievement of the LHC Run-I physics programme has been the discovery of a new spin-0 resonance (h) with a mass of 125 GeV and with properties consistent with that of the standard model (SM) Higgs boson [1][2][3]. Besides precision measurements of processes involving a h, the LHC Higgs physics programme however also includes a wide spectrum of searches for additional Higgses (a summary of LHC Run-I results can be found in [4] for instance). Such states are predicted in many SM extensions such as supersymmetry or models where the Higgs is realised as a pseudo Nambu-Goldstone boson (PNGB) of a new approximate global symmetry.
In fact, if the extended electroweak (EW) symmetry breaking sector contains a PNGB, this state can be significantly lighter than the other spin-0 particles. A well-known example of a model that includes a light pseudoscalar (a) is provided by the next-to-minimal supersymmetric SM (NMSSM) where this state can arise as a result of an approximate global U(1) R symmetry [5]. Since in this case the amount of symmetry breaking turns out to be proportional to soft breaking trilinear terms, the mass of the a can naturally be less than half of the SM Higgs mass, if the trilinear terms are dialled to take values in the GeV range. Non-supersymmetric theories that can feature a light pseudoscalar are, to just name a few, simplified models where a complex singlet scalar is coupled to the Higgs potential of the SM or the two-Higgs doublet model (2HDM), Little Higgs models and hidden valley scenarios (see [6] and references therein for details).
Irrespectively of the precise ultraviolet (UV) realisation, a light pseudoscalar can lead to distinctive collider signatures. The most obvious consequence are exotic decays of the SM Higgs, namely h → aa for m a < m h /2 [7,8] and h → aZ for m a < m h − m Z [6,9]. Another feature that can have important phenomenological implications is that in the presence of the heavy-quark transition a → bb (a → cc) the pseudoscalar a can mix with bottomonium (charmonium) bound states with matching quantum numbers [10][11][12][13][14][15][16].
LHC searches for h → aa have been performed in the 4µ [17,18], 4τ [19,20], 2µ2τ [20,21], 2µ2b [20] and 2τ2b final states [22]. The obtained results have been used to set upper bounds on the h → aa branching ratio in 2HDMs with an extra complex singlet (2HDM+S) for pseudoscalar masses in the range of [1, 62.5] GeV. The analyses [17,[19][20][21][22] however all exclude m a values in the regions [3,5] GeV and [9,11] GeV for which a-η c and a-η b mixing effects as well as open flavour decays to D and B (s) meson pairs can be potentially important.
The main goal of this work is to extend the latter results to the cc and bb threshold regions by including effects that cannot be properly described in the partonic picture. In order to highlight the complementarity of different search strategies for a light a, we also compare our improved limits to other bounds on the 2HDM+S parameter space that derive from the LHC searches for h → Z d Z → 4 [23], h → Z d Z → 2µ2 [24], pp → a → µ + µ − [25,26], pp → abb followed by a → τ + τ − [27] or a → µ + µ − [28], pp → a → γγ [29,30], pp → a → τ + τ − [31], from the BaBar analyses of radiative Υ decays [32][33][34] and from the LHCb measurements of the production of Υ mesons [15,35] as well as the inclusive dimuon cross section [36,37].
This article is organised as follows. In Section 2 we briefly recall the structure of the 2HDM+S scenarios. Our recast of the results [17,[19][20][21][22] is presented in Section 3, where we also derive the constraints on the 2HDM+S parameter space that follow from the measurements and prosposals [15,[23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]. We conclude in Section 4. The formulas necessary to calculate the partial decay widths of the pseudoscalar a are collected in Appendix A, while Appendix B contains a concise discussion of the mixing formalism and of open flavour decays that are relevant in the vicinity of the bb and cc thresholds.

Theoretical framework
In the following section we will interpret various searches for light pseudoscalars in the context of 2HDM+S scenarios. In this class of models a complex scalar singlet S is added to the 2HDM Higgs potential (see e.g. [38,39] for 2HDM reviews). The field S couples only to the two Higgs doublets H 1,2 but has no direct Yukawa couplings, acquiring all of its couplings to SM fermions through its mixing with the Higgs doublets. A light pseudoscalar a can arise in such a setup from the admixture of the 2HDM pseudoscalar A and the imaginary part of the complex singlet S . The corresponding mixing angle will be denoted by θ, and defined such that for θ → 0 the mass eigenstate a becomes exactly singlet-like.
In order to eliminate phenomenologically dangerous tree-level flavour-changing neutral currents (FCNCs) the Yukawa interactions that involve the Higgs fields H 1,2 have to satisfy the natural flavour conservation hypothesis [40,41]. Depending on which fermions couple to which doublet, one can divide the resulting 2HDMs into four different types. In all four cases the Yukawa couplings between the pseudoscalar a and the SM fermions take the generic form Here y f = √ 2m f /v denote the SM Yukawa couplings and v 246 GeV is the EW vacuum expectation value. The parameters ξ M f encode the dependence on the 2HDM Yukawa sector and the factors relevant for the further discussion are given in Table 1. In this table the shorthand notations s θ = sin θ and t β = tan β have been used. Similar abbreviations will also be used in what follows. type   I  II  III  IV up-type quarks s θ /t β s θ /t β s θ /t β s θ /t β down-type quarks −s θ /t β s θ t β −s θ /t β s θ t β charged leptons −s θ /t β s θ t β s θ t β −s θ /t β Table 1.
Ratios ξ M f of the Yukawa couplings of the pseudoscalar a relative to those of the SM Higgs in the four types of 2HDM+S models without tree-level FCNCs.
In the presence of (2.1) the CP-odd scalar a can decay into fermions at tree level and into gluons, photons and EW gauge bosons at loop level. The expressions for the partial decay widths Γ(a → XX) that we employ in our study are given in Appendix A. Since in this work we will assume that the a is lighter than the W, Z, h and the other 2HDM Higgs mass eigenstates H, A, H ± , decays of the a into the latter states are kinematically forbidden.
If the a is sufficiently light, exotic decays of the SM Higgs into the two final states aZ and aa are however possible. The partial decay width Γ(h → aZ) is in 2HDM+S scenarios entirely fixed by the 2HDM parameters α, β and the mixing angle θ. Explicitly, one has at tree level and Notice that in the exact alignment/decoupling limit, i.e. α = β − π/2, in which the lighter CP-even spin-0 state h of the 2HDM becomes fully SM-like, the coupling g haZ and thus Γ(h → aZ) is precisely zero. However, given that the total decay width of the SM Higgs is only about 4 MeV, the process h → aZ can be important even if deviations from the alignment/decoupling limit are relatively small. Unlike g haZ , the triple Higgs coupling g haa depends not only on the physical Higgs masses and mixing angles but also on some of the trilinear couplings that appear in the full scalar potential. This feature makes the partial decay width Γ(h → aa) model dependent, and in consequence the two exotic branching ratios BR(h → aZ) and BR(h → aa) can be adjusted freely by an appropriate choice of parameters. Following this philosophy we will treat BR(h → aZ) and BR(h → aa) as free parameters in the remainder of this article.
To facilitate a comparison between the results obtained by the CMS collaboration and by us, we consider like [20] the following four 2HDM+S benchmark scenarios: the type I model with t β = 1, the type II model with t β = 2, the type III model with t β = 5 and the type IV model with t β = 0.5. The fermionic coupling factors ξ M f corresponding to each 2HDM+S type are reported in Table 1. It is important to realise that the s θ -dependence of ξ M f cancels in BR(a → XX) and it is thus possible to translate constraints on signal strengths such as The results of our recast are shown in the panels of Figure 1 and should be compared to the exclusion plots displayed in Figure 8 of [20]. The branching ratios BR(a → XX) used to interpret the results in the four particular 2HDM+S scenarios are calculated using the formulas given in Appendix A and include the mixing and threshold effects described in Appendix B. Notice that the inclusion of a-η c and a-η b mixing is crucial to obtain meaningful predictions in the m a regions [3,5] GeV and [9,11] GeV, which are left unexplored in the CMS analysis [20].
While overall we observe good agreement between the 95% confidence level (CL) exclusions set by CMS and by us, some differences in the derived limits are evident. Firstly, our analysis covers the mass region close to the cc (bb) threshold, where our limits display a resonance-like behaviour as a result of the mixing of the a with the three η c (six η b ) states included in our study. Second, in the m a range of [1,3] GeV our bounds on µ h BR(h → aa) tend to be somewhat weaker than those derived in [20]. The observed difference is again a consequence of the mixing of the a with QCD bound states. In fact, in the very low mass range the total decay width of the unmixed a is below 10 −3 MeV in the considered 2HDM+S scenarios, while that of the lightest η c state amounts to around 30 MeV [42]. Hence even a small η c -admixture in the mass eigenstate a can lead to an enhanced total decay width Γ a which in turn results in a suppression of BR(a → µ + µ − ) and a weakening of the bound on µ h BR(h → aa).
A light pseudoscalar a can also be searched for via the decay h → aZ. The only LHC analyses that presently can be used to set bounds on BR(h → aZ) are the ATLAS searches for new dark bosons Z d produced in h → Z d Z [23,24]. Notice that while the Z d decays democratically into electrons and muons in the case of the a one has Γ(a → e + e − )/Γ(a → µ + µ − ) = m 2 e /m 2 µ 2 · 10 −5 . As a result 4e and 2e2µ events originating from h → aZ → 4e and h → aZ → 2e2µ give essentially no contribution to the signal strength in pp → h → aZ → 4 . The 8 TeV ATLAS study [23] however only provides exclusion bounds on BR(h → Z d Z → 4 ) from a combination of final states. To correct for this mismatch we have calculated r Aε = X=4µ,2µ2e Aε X / X=4µ,4e,2e2µ,2µ2e Aε X , where Aε X denotes the product of acceptance and reconstruction efficiency in the final state X -the values for Aε X can be found in the auxiliary material of [23]. We find that r Aε has only a mild dependence on m a and amounts to around 60%. The actual limits are then obtained by and solving for BR(h → aZ). To improve upon this naive recast one would need individual bounds for the different combinations of final-state lepton flavours. In fact, the very recent 13 TeV ATLAS analysis [24] provides Aε 2µ2 as well as limits on the relevant fiducial cross section. Our recast of the latter results thus only has to rely on the assumption that the product Aε 2µ2 is roughly the same for the Z d model and the 2HDM+S scenario, which we indeed believe to be the case. type IV, t β = 0.5 Figure 1. Limits on µ h BR(h → aa) in the 2HDM+S of type I with t β = 1 (top left), type II with t β = 2 (top right), type III with t β = 5 (bottom left) and type IV with t β = 0.5 (bottom right). The purple, blue, orange, red, cyan, green and dark red exclusions correspond to the search for h → aa → 4µ [17], h → aa → 4τ [19], h → aa → 4τ [20], h → aa → 2µ2τ [20], h → aa → 2µ2τ [21], h → aa → 2µ2b [20] and h → aa → 2τ2b [22], respectively. The dashed black lines indicate µ h BR(h → aa) = 1 and all coloured regions are excluded at 95% CL.
The exclusion limits on µ h BR(h → aZ) corresponding to the four 2HDM+S benchmark scenarios discussed earlier are presented in Figure 2. From the panels it is evident that, apart from pseudoscalar masses around 25 GeV where the data [24] has a local deficit, the constraints that derive from the 13 TeV analysis [24] are significantly stronger than those that one obtains from the 8 TeV data [23]. One also observes that the constraints in the first and second benchmark are weak as they just start to probe the region µ h BR(h → aZ) 1, whereas in the third and fourth 2HDM+S scenario already values of µ h BR(h → aZ) 0.1 can be probed with the available LHC data sets. . Limits on µ h BR(h → aZ) in the 2HDM+S of type I with t β = 1 (top left), type II with t β = 2 (top right), type III with t β = 5 (bottom left) and type IV with t β = 0.5 (bottom right). The red and green bounds correspond to the ATLAS search for pp → h → Z d Z → 4 [23] and pp → h → Z d Z → 2µ2 [24], respectively. The dashed black lines indicate µ h BR(h → aZ) = 1 and all coloured regions are excluded at 95% CL.
Since the asymmetry between electron and muon final states from h → aZ decays is a striking signature of a light pseudoscalar, we strongly encourage our experimental colleagues to provide as in [24] separate bounds for the 2e2 and 2µ2 final states in future searches for signatures of the type h → Z d Z → 4 .
For concreteness we study the same four 2HDM+S scenarios that we have already considered before. The most stringent limits on |s θ | that can be derived at present are displayed in Figure 3.
In order to recast the results of the CMS searches for a → µ + µ − [26], pp → abb → τ + τ − bb [27], pp → a → γγ [29], pp → a → τ + τ − [31], the LHCb measurements of Υ production [15,35] and the inclusive dimuon cross section [37], one needs to know the production cross sections of a light a in gluon-fusion and in association with bb pairs. Our predictions for gg → a production are obtained at next-to-next-to-leading order in QCD using HIGLU [43], while the pp → abb cross sections are calculated at next-to-leading order (NLO) in QCD in the four-flavour scheme with MadGraph5_aMCNLO [44] employing an UFO implementation [45] of the 2HDM model discussed in the publication [46].
Our recast of the results of the LHCb search for dark photons A [37] proceeds as follows. We calculate the inclusive pp → A production cross section at NLO in QCD with the help of MadGraph5_aMCNLO [44], while we extract BR(A → µ + µ − ) from the well-measured cross section ratio R = σ(e + e − → hadrons)/σ(e + e − → µ + µ − ) [42]. Following [36,37], model-dependent A -Z mixing effects are included in our calculation employing the formulas given in [47]. We have also taken into account detector acceptance differences between pp → A → µ + µ − and pp → a → µ + µ − by computing the ratio r A = A a /A A of signal acceptances. We find that r A amounts to around 2.0, 1.3, 1.0 at m a = 5 GeV, 15 GeV, 70 GeV and scales approximately linear between the quoted m a values. Concerning the detection efficiencies ε A and ε a we assume that they are identical for A → µ + µ − and a → µ + µ − , which should be a good approximation when the dimuon signal is prompt [37]. We finally add that in our recast of the LHCb dark photon results, we only consider the mass region m a > 4.5 GeV to avoid a-η c mixing contributions to the pp → a cross section associated to pp → η c production. The mass region m a ∈ [9.1, 10.6] GeV is also not covered by our recast, because in [37] the LHCb collaboration does not present bounds on the kinetic mixing of the A close to the bb threshold.
The main conclusion that can be drawn from the results presented in Figure 3 is that only in the 2HDM+S scenario of type IV with t β = 0.5 it is possible to set physical meaningful bounds on the sine of the mixing angle θ, i.e. |s θ | < 1, over the entire range of studied pseudoscalar masses. One furthermore observes that solely the BaBar search for the radiative decay Υ(1S ) → aγ → µ + µ − γ [33] allows to probe parameter regions with |s θ | < 0.1. This search is however kinematically limited to m a < m Υ(1S ) 9.5 GeV. Improvements in the existing LHC search strategies (and/or new approaches) are needed to reach the same sensitivity on |s θ | for pseudoscalar masses above approximately 10 GeV in the examined 2HDM+S benchmark models. Measurements of the inclusive dimuon cross section [36,37] seem to be quite promising in this context.

Conclusions
Beyond the SM theories with an extended Higgs sector can naturally lead to pseudoscalar resonances with masses significantly below the EW scale if these states serve as PNGBs of an approximate global U(1) symmetry. The R-symmetry limit in the NMSSM and the case of spontaneously broken U(1) subgroups in Little Higgs models are just two working examples of this general idea. Searches for light CP-odd spin-0 states are thus theoretically well-motivated and in the case of a detection could help to illuminate the structure and dynamics of the underlying UV model.
The existing collider searches for pseudoscalars with masses of approximately [1, 100] GeV fall into two different classes. Firstly, searches that look for the presence of a light a in the decay of a SM particle. Searches for h → aa and h → aZ, but also the radiative decays Υ → aγ belong to this category. In the case of the exotic Higgs decays the resulting signature that the ATLAS and CMS experiments have explored are four-fermion final states containing at least two opposite-sign leptons [17][18][19][20][21][22][23][24], while what concerns the radiative Υ decays, BaBar has considered the hadronic, dimuon and ditau decays of pseudoscalars [32][33][34]. The second type of searches instead relies on the direct production of the a in pp collisions and its subsequent decays to either charged lepton or photon pairs. Both the gluon-fusion channel [15,25,26,[29][30][31]35] and abb production [27,28] have so far been exploited to look for light pseudoscalars at the LHC in this way.
In this work, we have performed a global analysis of the present collider constraints on light pseudoscalar states. To facilitate a comparison with the recent CMS study [20], we have considered the class of 2HDM+S models, treating the parameters t β and s θ as well as the branching ratios BR(h → aa) and BR(h → aZ) as free parameters -see Section 2 for a concise introduction to the 2HDM+S setup. A complication that arises in our analysis is that in the mass regions [3,5] GeV and [9,11] GeV, non-perturbative effects such as the mixing of the pseudoscalar with QCD bound states have to be taken into account to allow for a meaningful interpretation of the experimental data. We have worked out the theoretical formalism necessary to calculate the most relevant shortdistance and long-distance effects and provide a collection of the corresponding formulas in the two Appendices A and B.
Our numerical analysis consists of three parts. In the first part, we have derive 95% CL exclusion limits on the signal strength µ h BR(h → aa) that follow from the latest CMS searches for the exotic h → aa decay [17][18][19][20][21][22], while in the second part we present the limits on µ h BR(h → aZ) that stem from the ATLAS searches for h → Z d Z → 4 [23] and h → Z d Z → 2µ2 [24]. The exclusion bounds on |s θ | that arise from the searches [15,26,27,29,31,33,35,37] are finally derived in the third part of our numerical study. In all three cases, we have considered four specific 2HDM+S benchmark scenarios that differ in the choice of Yukawa sector and t β . We have found that the inclusion of a-η c a-η b mixing effects as well as open flavour decays to D B (s) meson pairs has a visible impact on the obtained limits only in the mass region of approximately [1,4] GeV [10,15] GeV , while perturbative calculations are perfectly adequate for m a values away from the cc and bb thresholds.
The main conclusion that can be drawn from the results presented in Figures 1, 2 and 3 is that existing collider constraints on the parameter space of 2HDM+S models are in general not very strong. Exceptions are the [1,3] GeV region in which µ h BR(h → aa) is well-constrained by the CMS search for h → aa → 4µ [17] and the [1, 9.5] GeV range where the Υ(1S ) → aγ → µ + µ − γ search of BaBar [33] provides stringent limits on |s θ |. Much to the opposite, the 2HDM+S parameter space turns out to be least constrained for m a values in the range of approximately [15,70] GeV. The development of improved or new search techniques (such as for instance dedicated searches for h → aZ [24] and inclusive diphoton [30] or dimuon [36,37] cross section measurements) that specifically focus on the latter mass region therefore seems to be a worthwhile scientific goal.
Mike Williams for their interest in our work, constructive feedback and their useful suggestions. UH appreciates the continued hospitality and support of the CERN Theoretical Physics Department. JFK acknowledges the financial support from the Slovenian Research Agency (research core funding No. P1-0035 and J1-8137). MS would like to thank the organisers of Les Houches for the great and fruitful atmosphere of the workshop.

A Decay width formulas
In the calculation of the total decay width Γ a of the unmixed pseudoscalar a, we employ the following expressions for the partial decay widths (see the reviews [49][50][51][52] for instance) where MS masses are indicated by a bar while masses without a bar are evaluated in the pole scheme. We have furthermore defined τ f /a = 4m 2 f /m 2 a and β f /a = 1 − τ f /a and used the symbol Q q to denote the electric charge of the quark in question. All MS masses as well as the coupling constants α s and α are renormalised at the scale µ R = m a . Table 1 finally contains the coupling assignments ξ M f that we consider in our work. The QCD corrections to the partial decay width into light quarks (A.2) that are included in our numerical analysis read [10,[53][54][55][56][57][58][59][60][61][62][63] and [64,65] The symbol N f introduced above denotes the number of light quark flavours that are active at the scale m a . For pseudoscalar masses far above the threshold, i.e. m a 2m q , the results (A.6) and (A.7) represent at the moment the most accurate predictions for the QCD corrections to Γ(a → qq). In our numerical analysis, we hence use them to calculate the partonic rate of a → ss.
In the case of the partial decay width into heavy-quark pairs (A.3) the QCD corrections are given to first order in α s by [10,[53][54][55][56][57] Here we have introduced the abbreviation x β Q/a = (1 − β Q/a )/(1 + β Q/a ) and the one-loop function entering (A.8) takes the form with Li 2 (z) denoting the usual dilogarithm. In the threshold region, i.e. m a 2m Q , mass effects are important and as a result the QCD corrections (A.8) should be used to describe them. Following the prescription implemented in HDECAY [66,67], the transition between the region close to threshold to that far above threshold is achieved by a smooth where for analytic continuation it is understood that τ → τ − i0.
The multiplicative factor K g entering (A.4) takes the following form where the second term encodes the virtual two-loop QCD corrections, while the third term corresponds to the finite part of the real QCD corrections in the heavy-quark limit [49,68]. We have verified that quark mass effects of the real corrections not included in (A.11) amount to no more than 5%. The virtual corrections can be written as where y q/a = −x q/a with τ q/a → τ q/a +i0 for analytic continuation and the prime denotes a derivative with respect to τ q/a . To reproduce the position of the a → qq threshold correctly, we set µ q = m a /2 in our study. The two-loop function appearing in (A.12) reads [68,69] G(y) = y (1 − y) 2 48H(1, 0, −1, 0; y) + 4 ln(1 − y) ln 3 y − 24ζ 2 Li 2 (y) − 24ζ 2 ln(1 − y) ln y Here H(1, 0, −1, 0; z) is a harmonic polylogarithm of weight four with two indices different from zero, which we evaluate numerically with the help of the program HPL [70]. The polylogarithm of order three (four) is denoted by Li 3 (z) Li 4 (z) , while ζ 2 = π 2 /6, ζ 3 1.20206 and ζ 4 = π 4 /90 are the relevant Riemann's zeta values.
In the case of (A.5) we decompose the relevant QCD corrections as with [68,69,71]

B Mixing and threshold effects
Even though the decay a → bb (a → cc) is kinematically forbidden below the open-flavour threshold, the presence of heavy quarks can become relevant through mixing between the pseudoscalar a and bottomonium (charmonium) bound states with the same quantum numbers [10][11][12][13][14][15][16]. Such mixings can effectively be described through off-diagonal contributions δm 2 aη b (n) to the pseudoscalar mass matrices squared. In the case of a-η b mixing, we employ with The masses and radial wave functions of the η b (n) states are denoted by m η b (n) and R η b (n) , respectively. The latter quantities can be extracted from the Υ(n) leptonic decay widths (see [72] for instance) which are measured rather precisely [42]. In the case of a-η c mixing, we only include the first three states in the pseudoscalar mass matrix squared (B.1) and rely on the potential model calculations of [73] to determine the radial wave functions R η c (n) . The values of the η b (n) and η c (n) masses and radial wave functions that are used in our numerical analysis are collected in Table 2 for convenience.
To be able to determine the eigenvalues and eigenvectors of (B.1) one also needs to know the total decay widths of the η b (n) and η c (n) states. The digluon decay widths of the η b (n) states are given to leading order in α s by (see [10] for example) and an analogous formula holds in the case of the charmonium resonances. The partial decay widths (B.3) essentially saturate Γ η b (n) with n 5, 6. For η b (5) and η b (6), however, also decays to final states involving π and B (s) mesons are relevant. In the case of the decays to pion final states, we employ [42] Γ(η b (5) → π mesons) = 1. where Γ(ψ(4040) → D mesons) Γ ψ(4040) = 80 MeV [42]. We furthermore emphasise that the branching ratios η b (n) → µ + µ − are all below the 10 −10 level [15] and therefore can be safely ignored in the mixing formalism. The effects of the ditau decays of the bottomonium bound states are negligible as well and so are the dilepton decays of the η c (n) mesons. Effects of a-η b mixing in h → aa such as for instance h → 2η b (n) → aa are part of BR(h → aa) and thus effectively included in our numerical analysis. The same is true for contributions of intermediate η c (n) states to the exotic decay h → aa of the SM Higgs.
Above the bb (cc) threshold a perturbative description of the production and the decay of the pseudoscalar a breaks down. In this region one can however approximate the bb (cc) contributions to the total decay width Γ a through a heuristic model that is inspired by QCD sum rules [10,14,15] and interpolates to the continuum sufficiently above threshold. The interpolations take the form with m B = 5.28 GeV, m B * = 5.33 GeV, m D = 1.86 GeV and m D * = 2.01 GeV [42]. In our analysis, the interpolation is achieved by simply multiplying the partonic decay width Γ(a → bb) and Γ(a → cc) by the factor N b a and N c a , respectively. For m a > 2m K decays into kaons become kinematically allowed. The decay a → KK however violates CP, and as a result a can in practice only decay into three-body final states such as KKπ. Following [48], we estimate the hadronic width Γ(a → ss → KKπ) by multiplying Γ(a → ss) by the suppression factor with λ (x, y, z) defined in (2.4).