A fresh look into m_{c,b} and precise f_{D_(s),B_(s)} from heavy-light QCD spectral sum rules

Using recent values of the QCD (non-) perturbative parameters given in Table 1 and an estimate of the N3LO QCD perturbative contributions based on the geometric growth of the PT series, we re-use QCD spectral sum rules (QSSR) known to N2LO PT series and including all dimension-six NP condensate contributions in the full QCD theory, for improving the existing estimates of {m}_{c,b} and f_{D_(s)}, f_{B_(s)} from the open charm and beauty systems. We especially study the effects of the subtraction point on"different QSSR data"and use (for the first time) the Renormalization Group Invariant (RGI) scale independent quark masses in the analysis. The estimates [rigourous model-independent upper bounds within the SVZ framework] reported in Table 8: f_D/f_\pi=1.56(5)[<1.68(1)], f_B/f_\pi=1.58(5)[<1.80(3)] and f_{D_s}/f_K= 1.58(4) [<1.63(1)], f_{B_s}/f_K=1.50(3)[<1.61(3.5)], which improve previous QSSR estimates, are in perfect agreement (in values and precisions) with some of the experimental data on f_{D,D_s} and on recent lattice simulations within dynamical quarks. These remarkable agreements confirm both the success of the QSSR semi-approximate approach based on the OPE in terms of the quark and gluon condensates and of the Minimal Duality Ansatz (MDA) for parametrizing the hadronic spectral function which we have tested from the complete data of the J/\psi and \Upsilon systems. The values of the running quark masses m_c(m_c)=1286(66) MeV and m_b(m_b)= 4236(69) MeV from M_{D,B} are in good agreement though less accurate than the ones from recent J/\psi and \Upsilon sum rules.


Introduction and a short historical review
The (pseudo)scalar meson decay constants f P are of prime interests for understanding the realizations of chiral symmetry in QCD. In addition to the well-known values of f π = (130.4 (2) MeV and f K = 156.1(9) MeV [2] which control the light flavour chiral symmetries, it is also desirable to extract the ones of the heavy-light charm and bottom quark systems with highaccuracy. These decay constants are normalized through the matrix element: where: is the local heavy-light pseudoscalar current; q ≡ d, s; Q ≡ c, b; P ≡ D (s) , B (s) and where f P is related to the leptonic width: where m l is the lepton mass and |V Qq | the CKM mixing angle. Besides some earlier attempts based on non-relativistic potential models to extract these quantities (which are however not applicable for the heavy-light systems), the first bounds on f D and f B from QCD spectral sum rules (QSSR) [3] 1 were derived by NSV2Z [8], which have been improved four years later in [9][10][11]. Since then, but long time before the lattice results, different QSSR papers have been published in the literature for estimating f D,B 2 . These results look, at first sight, in disagreement among each others and some of them, claimed the observation of the scaling f P ∼ 1/ √ M P expected in the large M P limit [12]. These different papers have been scrutinized in [6,13], where Narison found that the apparent discrepancies between the different results can be solved if one applies carefully the stability criteria (also called sum rule window) of the results versus the external QSSR Laplace/Moments sum rules variables and continuum threshold t c . In this way, and for given values of m c,b , he obtained the values: which are independent of the forms of the sum rules used. However, the result has been quite surprising as it indicates a large violation of the heavy quark symmetry scaling predictions,where 1/M Q corrections have been estimated in [14]. This "unexpected result" has been confirmed few years later by lattice calculations [15]. Since then, some progresses have been done for improving the QCD expression of the 2-point correlator. It starts from a confirmation of the SVZ original expression of the LO perturbative and non-perturbative contributions. Then, Broadhurst and Generalis [10,16] have provided 1 For reviews, see e.g: [4][5][6][7]. 2 For reviews and more complete references, see e.g: [5,6].
the complete PT α s NLO including light quark mass corrections. It has been completed by the PT α 2 s N2LO corrections of Chetyrkin and Steinhauser [17] in the case of one heavy and one massless quarks. This result has been completed by the inclusion of the NP contributions up to dimension-six [11] and of the light quark mass corrections to LO by [11,18]. All of these previous QCD expressions have been given in terms of the onshell quark mass. In [19], Narison has used (for the first time) the running c, b quark masses in the QSSR analysis, by using its known relation with the on-shell mass known at present to NLO [20][21][22], N2LO [10,16] and N3LO [24] where it has been noticed that the QSSR PT expressions converge faster. It has also been noticed that the values of f D,B are very sensitive to the value of m c,b motivating him to extract m c,b (for the first time) from the known values of M D and M B . Recent analysis, including the α 2 s corrections have been presented in the literature, in the full theory where the running MS mass has been used [25][26][27] and in HQET [28] where the radiative corrections are large due to the uses of the on-shell mass 3 . In the following, we shall present analysis based on the full QCD theory where we use as inputs the most recent values of the (non-)perturbative QCD parameters given in Table 1. We assume the geometric growth of the PT series [30] as a dual to the effect of a 1/q 2 term [31,32] for an estimate of the N3LO perturbative contributions. We shall also study systematically the effect of the substraction points on each "QSSR data" and use (for the first time) in the analysis, the Renormalization Group Invariant (RGI) s, c, b quark masses introduced by [33] and which are scale and (massless) scheme independent.

QCD spectral sum rules (QSSR)
• The Laplace sum rules (LSR) We shall be concerned with the two-point correlator : where Jq Q (x) is the local current defined in Eq. (2). The associated Laplace sum rules (LSR) Lq Q (τ) and its ratio Rq Q (τ) read [3] 4 : where µ is the subtraction point which appears in the approximate QCD series when radiative corrections are included. The ratio of sum rules Rq Q (τ, µ) is useful, as it is equal to the resonance mass squared, in the Minimal Duality Ansatz (MDA) parametrization of the spectral function: (8) where f P is the decay constant defined in Eq. (1) and the higher states contributions are smeared by the "QCD continuum" coming from the discontinuity of the QCD diagrams and starting from a constant threshold t c .
• The Q 2 = 0 moment sum rules (MSR) We shall also use for the B-meson, the moments obtained after deriving n + 1-times the two-point function and evaluated at Q 2 = 0 [3], where an expansion in terms of the on-shell mass M b can be used. They read: and the associated ratio: • Test of the Minimal Duality Ansatz (MDA) from J/ψ and Υ We have checked explicitly in [6] that the MDA presented in Eq. (8), when applied to the ρ-meson reproduces within 15% accuracy the ratio Rd d measured from the total cross-section e + e − → I = 1 hadrons data ( Fig. 5.6 of [6]). In the case of charmonium, we have also compared M 2 ψ from R (n) cc with the one from complete data and find a remarkable agreement for higher n ≥ 4 values ( Fig. 9.1 of [6]), indicating that for heavy quark systems the rôle of the QCD continuum will be smaller than in the case of light quarks and the exponential weight or high number of derivatives suppresses efficiently the QCD continuum contribution but enhances the one of the lowest ground state in the spectral integral. We redo the test done for charmonium in Fig. 9.1 of [6] and analyze the bottomium channel for the LSR and MSR. We show in Fig. (1a) the τ-behaviour of the ratio of L exp cc normalized to L dual cc where we have used the simplest QCD continuum expression for massless quarks to order α 3 s from the threshold t c 5 : We show in Fig. (1b) the τ-behaviour of M ψ , where the continuous (oliva) curve corresponds to √ t c ≃ M ψ(2S ) − 0.15 GeV. We show a similar analysis for the bottomium sum rules in Fig. (2) for the LSR and in Fig. (3) for the MSR where we have taken √ t c ≃ M Υ(2S ) − 0.15 GeV. One can see that the MDA, with a value of √ t c around the value of the 1st radial excitation mass, describes quite well the complete data in the region of τ and n where the corresponding sum rules present τ or n stability [35]: as we shall see later on. This good description of the data by the MDA shows the efficient rôle of the exponential weight or high number of derivatives for suppressing the higher mass states and QCD continuum contribution in the analysis. This nice feature prevents the introduction of some more involved models bringing new parameters in the analysis where some of them cannot be understood from QCD 1st principles. Moreover, MDA has been also used in [36] (called Minimal Hadronic Ansatz in this paper) in the context of large N c QCD, where the restriction of an infinite set of large N c narrow states to a Minimal Hadronic Ansatz which is needed to satisfy the leading short-and long-distance behaviours o the relevant Green's functions, provides a very good approximation to the observables one compute.  of moments Rd Q for the harmonic oscillator as a function of the imaginary time variable τ, where one knows the exact and approximate results, one can find [37] that the exact energy E 0 of the ground state can be approached from above by the approximate series (see Fig. 4). At the minimum or inflexion point (stability) of the curves, one has a ground state dominance. For small time (large Q 2 ), all level contributes, while for large time (small Q 2 ) the series breakdown. We shall apply this stability criterion inspired from quantum mechanics in our analysis. In principle, the continuum threshold √ t c in Eq. (8) is a free parameter, though one expects its value to be around the mass of the 1st radial excitation because the QCD spectral function is supposed to smear all the higher state contributions in the spectral integral as explicitly shown previously in Section 2. In order to avoid the model-dependence on the results, Refs. [5,6,13,14,19,25] have considered the conservative range of t c -values where one starts to have τ-or n-stability until which one reaches a t c -stability where the contribution of the lowest ground state to the spectral integral completely dominates. For the D and B mesons, this range is [5,6,13,14,19,25]:

The QCD input parameters
The QCD parameters which shall appear in the following analysis will be the strange, charm and bottom quark masses m s,c,b (we shall neglect the light quark masses q ≡ u, d), the light quark condensate qq , the gluon condensates g 2 G 2 ≡ g 2 G a µν G µν a and g 3 G 3 ≡ g 3 f abc G a µν G b νρ G c ρµ , the mixed condensate qgσGq ≡ qgσ µν (λ a /2)G a µν q = M 2 0 qq and the fourquark condensate ρ qq 2 , where ρ ≃ 2 indicates the deviation from the four-quark vacuum saturation. Their values are given in Table 1 and we shall work with the running light quark parameters known to order α 3 s [5,6,38]. They read: where β 1 = −(1/2)(11 − 2n f /3) is the first coefficient of the β function for n f flavours; a s ≡ α s (τ)/π;m q,Q is the RGI quark mass,μ q is spontaneous RGI light quark condensate [33]. The QCD correction factor C(a s ) in the previous expressions is numerically: which shows a good convergence. We shall use: from τ-decays [39,40], which agree perfectly with the world average 2012 [41,42]: We shall also use the value of the running strange quark mass obtained in [43] 6 given in Table 1. The value of the running qq condensate is deduced from the value of (m u + m d )(2) = (7.9±0.6) MeV obtained in [43] and the well-known GMOR relation: The values of the running MS mass m Q (M Q ) recently obtained in Ref. [35] from charmonium and bottomium sum rules, will also be used 7 . Their average is given in Table 1. From which, we deduce the RGI invariant heavy quark masses to order α 2 s , in units of MeV: For the light quarks, we shall use the value of the RGI mass and spontaneous mass to order α s for consistency with the known α s m s and qq condensate corrections of the two-point correlator. They read, in units of MeV:  [39,48,51]

QCD expressions of the sum rules
• The LSR To order α 2 s , the QCD theoretical side of the sum rule reads, in terms of the on-shell heavy quark mass M Q and for m d = 0: where: µ is an arbitrary subtraction point; R 2s is the α 2 s -term obtained semi -analytically in [17] and is available as a Mathematica package program Rvs.m. Neglecting m d , the PT NLO terms read [10]: The contribution up to the d = 4 gluon condensate and up to d = 6 quark condensates have been obtained originally by NSV2Z [8]. The contribution of the d = 6 g 3 f abc G 3 and j 2 gluon condensates have been deduced from the expressions given by [11] (Eqs. II.4.28 and Table II.8) where: after the use of the equation of motion. ρ ≃ (2 ± 0.2) measures the deviation from the vacuum saturation estimate of the d = 6 quark condensates [39,48,51]. The α s correction to d d , in the MS -scheme, comes from [26], where the running heavy quark mass m Q enters into this expression. Using the known relation between the runningm Q (µ) and on-shell mass M Q in the MS -scheme to order α 2 s [20][21][22][23][24]: for n l light flavours, one can express all terms of the previous sum rules with the running mass m Q (µ). It is clear that, for some non-perturbative terms which are known to leading order of perturbation theory, one can use either the running or the pole mass. However, we shall see that this distinction does not affect, in a visible way, the present result, within the accuracy of our estimate, as the non-perturbative contributions are relatively small though vital in the analysis.

• The MSR
The moments read for m d = 0: where:

Estimates of f P andm
After inspection, one finds that f P and the RGI massm Q can only be simultaneously determined from Ld Q (τ, µ) and . This particular value of µ = τ −1/2 is also interesting because the subtraction scale moves with the sum rule variable τ in the analysis.
• Analysis of the τ-and t c -stabilities of Ld c and Rd c Using the central values of the QCD input parameters in Table  1 and in Eqs. (16), (18) and (19), one can show in Fig. 5 the influences of τ and t c on the value of f D and M D for a given value of the subtraction point µ = τ −1/2 , where, the τ-stability for f D is reached for: When extracting the RGI massm c from Rd c by requiring that it reproduces the experimental mass squared M 2 D , one can notice in Fig. (5) that, unlike f D , M D present τ-stability for larger range of t c -values: The existence of τ-stability at values of t c below 5.3 GeV 2 depends on the heavy quark mass value and disappears when we require the sum rule to reproduce the value of M D , such that we shall not consider a such region. The values of t c ≃ (6.5 ∼ 9.5) GeV 2 given in Eqs. (27) and (28)  contributions. (see also similar cases of the harmonic oscillator in Fig. 4 and of the Laplace sum rules for charmonium and bottomium in [35,37]). Like in earlier versions of this work [13,14,19,25], we consider this large range of t c -values in the aim to extract the most conservative result from the analysis and to avoid, in the same way, any (ad hoc) external input for fixing the exact value of t c . This procedure implies a larger error in our result than often quoted in the literature where (to my personal opinion) the systematics have been underestimated. A similar procedure will be done in the following and for the B-meson channel.
• Analysis of the convergence of the QCD series We study the convergence of the QCD series in the case of the charm quark at a such low value of the subtraction point µ = τ −1/2 and taking t c = 6 GeV 2 . We work in the MS -scheme as we know from previous works [19] that the PT series converge better than using the on-shell subtraction. In so doing, we estimate the α 3 s N3LO contribution using a geometric PT series as advocated in [30] which is dual to the effect of the 1/q 2 term when large order PT series are resummed. We show the τ-behaviour of f D in Fig. (6). One can notice that, all corrections act in a positive way. The prediction increases by about 17% from LO to NLO and another 14% from NLO to N2LO but remains unaffected by the inclusion of the N3LO contribution estimated above. These features indicate that the PT series converge quite well at this low scale, while the size of each PT corrections are reasonably small and will be even smaller for higher values of the subtraction point µ and for the B-meson. Therefore, a confirmation of this N3LO estimate requires an explicit evaluation of this contribution. As far as the non-perturbative contributions are concerned, their effects are relatively small.

• QCD and systematic error estimates
Using the previous QCD input parameters and their corresponding errors, we deduce the different errors on f P andm Q given in Table 2, where the optimal results have been taken at the τ-and t c -stability regions mentioned in the previous subsection: As mentioned earlier, we consider a such large range of t cvalues in the aim to extract the most conservative result from the analysis. However, this procedure induces a larger error in the analysis than the one quoted in the literature using some other models or using some other criteria. In fact, the range of values of our result includes most of the predictions given in the literature which are often quoted with smaller errors. Therefore, we expect that, within this procedure, we take properly into account most of the systematics of the sum rule approach.
In so doing, we take the central value of f D in Table 2 as coming from an arithmetic average of its values from the different t c given in the legend of Fig. (5) inside the range given by Eq. (27). We may have improved the accuracy of our predictions by introducing more model-dependent new parameters for parametrizing the continuum contribution, which we would not do as, in addition to the test performed in Section 2, we also want to check the degree of accuracy of the MDA parametrization for the heavy-light systems by confronting the results obtained in this paper with the some known data on f P or from lattice simulations. Indeed, such tests are important as the MDA model is widely used in the literature for predicting some not yet measured masses of new exotic hadrons like four-quark, molecules [57] and hybrid [58] states. However, we do not also try to fix more precisely t c by e.g. using Finite Energy Sum Rule [53] like did the authors in Ref. [59] as we want to have more conservative results. •

Results for f D andm c
Considering the common range of t c -values for f D and M D given in Eq. (27), we obtain the results quoted in Table 2 which we consider as improvement of the result obtained from the same sum rule and at the same subtraction point by [19] 9 : The smaller errors in the present analysis, come from more precise input parameters, more complete NP-corrections included into the OPE and more constrained range of t c -values.The value 8 Using the larger range of t c -values, we would have obtained a slightly different value:m c ≃ 1492(102) tc (82) qcd MeV, where the errors come respectively from the choice of t c and QCD parameters given in details in Table 2. 9 An extended discussion about the value of f D at different subtraction points will be done in the next section. obtained in Eq. (30) also agrees within errors with the accurate determination from charmonium systems quoted in Table 1 though less accurate. The main sources of errors from the present determination can be found in Table 2. One can notice that the contributions of the d = 6 condensates are negligible for f D (less than 0.3 MeV) and small for m c ( d d 2 and G 3 which contribute respectively to 17 and 6 MeV). Table 2: Central values and corresponding errors for f P andm Q in units of MeV from the LSR at the subtraction point µ = τ −1/2 . We have usedm Q in Eq. (18) for getting f P . The +(resp. -) sign means that the values of f P ,m Q increase (resp. decrease) when the input increases (resp. decreases). The relative change of sign from c to b in some errors is due to the effects of τm 2 Q appearing the OPE. Notice that the error in G 2 also affects the G 3 contribution. The Total error comes from a quadratic sum. show in Fig. (7) the τ-behaviours of f B and M B for different values of t c . One can see, that in this channel, τ-stability for f B is reached at 10 : 10 The apparent minima at τ ≤ 0.1 GeV 2 obtained for lower values of tc corresponds to the region where the continuum contribution to the spectral integral is dominant and should not be considered.
which we again consider as improvement of the result from [19]: obtained from the same sum rule.

Effects of the subtraction point on f D,B from LSR
The choice of subtraction points is also one large source of errors and discrepancies in the existing literature. In order to cure these weak points, we extract the values of f D,B and the corresponding errors at a given value of the subtraction point µ. We show in Fig. (9 10) and (14), the set of "QSSR data points " obtained in this way for different values of µ.

• Final results for f D and f B from LSR
Using the fact that the "physical observable" is independent of µ, we average (fit horizontally) the different data points of f D from LSR in Tables 2 and 3 and Fig. (10) (red triangle). The average is represented by the horizontal band in Fig. (10). The narrower (grey) domain corresponds to the resulting averaged error, while the larger one corresponds to the case where the error from the most precise determination has been taken. A similar analysis for f B from LSR has been done using the data  in Tables 2 and 4 and Fig. (14). We deduce from this analysis, the final results: where the quoted errors are the averaged errors. The previous errors are multiplied by about 1.8 for f D and 1.65 for f B if one keeps the errors from the most precise determinations.

• Final value ofm b from LSR
In addition to the sum rule Rd b subtracted at µ = τ −1/2 , we also notice that the sum rule Rd b subtracted at µ = M b , where the log (µ/M b )-term disappears in the QCD expression, presents τstability [see Fig. (11)] and can then provide another estimate ofm b . The result is given in Table 5. Taking the average of this result with the previous one in Table 2   when the N3LO term is included.The convergence of the PT series is comparable with the one of LSR shown in Fig. (8). •

Optimal values of f B andm b from MSR
Using a similar procedure as for the LSR, we study, in the case of MSR, the n-and t c -stabilities of f B andm b for different values of the subtraction point µ. The analysis is illustrated in Fig.  (13). The results are shown in Tables 4 and 5. One can notice that the sum rule does not stabilize for µ < 2 GeV, while for other values of µ, the range of values t c at which the n-stability is reached depends on the value of the subtraction point µ and are inside the range 32-42 GeV 2 . We show the results in Table  4 an in Fig. (14) from which we deduce the result from the moments in units of MeV:

Final values of f D , f B andm c,b
As a final result of the present analysis, we take the average of the results from LSR for f D andm c . This result is given in Eq.  which are: where we have used the more precise value of m b (m b ) given in Table 1 for getting f B . One can notice that f D ≃ f B , confirming previous results quoted in Eq. (4). This (almost) equality instead of the 1/ √ m b behaviour expected from HQET has been qualitatively interpreted in [60] using semi-local duality, while in [14] large mass corrections to the HQET lowest order expression have been found. These results are also confirmed by recent lattice calculations (see Table 8).

SU(3) breaking and estimates of f D s and f B s
We extend the previous analysis for extracting f D s ,B s by including the m s -corrections and by taking into account the S U (3) breaking of the quark condensate ss / d d given in Table 1.

• f D s from LSR
In so doing, we use the complete PT expression in m s of the QCD spectral function given to order α s by [10].  expressions for N2LO and N3LO have been used. The nonperturbative contributions come from the expressions given by [11,18,26] Table 6 and Fig (17).

• f B s from LSR and MSR
In this case, we only use the PT expression to order α s of the QCD spectral function expanded up to order m 2 s which is given by [26]. The non-perturbative contributions are the same as in the case of f D s . We show in Fig. (19) the τ-behaviour and nbehaviour of f B s from LSR and MSR for different values of t c at given µ = 4 GeV. The results for different values of µ are given in Table 7 and Fig. (18).  •

Results for f D s and f B s
From the previous analysis , we deduce: which, with the help of the results in Eqs. (36) and (39) lead to: These results agree within the errors with the ones obtained by using the semi-analytic expressions of the correlator to order α s [61]: with data when available [2,62] and with recent lattice simulations (see Table 8).

Rigourous model-independent upper bounds on f D (s) ,B (s)
Upper bounds on f D has been originally derived by NSV2Z [8] and improved four years later in [9][10][11] and more recently in [27]. In this paper, we shall use LSR and the positivity of the continuum contributions to the spectral integral for obtaining the upper bounds on the decay constants. The procedure will be similar to the estimate done in previous sections where the optimal bound will be obtained at the minimum or inflexion and: These bounds are stronger than earlier results in [9][10][11], while the results for f D,D s agree (within the large errors quoted there) with the ones in [27]. These large errors come mainly from m c , µ and¯ dd . The previous bounds can be used for excluding some experimental data and some theoretical estimates.
In deriving these bounds, we have only used the positivity of the spectral function and we have checked that the SVZexpansion converges quite well both for the PT radiative and non-perturbative corrections such that the approximate series is expected to reproduce with a good precision the exact solution. This fact can be (a posteriori) indicated by the remarkable agreement of our estimates with the lattice results. In this sense, we may state that the upper bound obtained previously is rigourous (at least within the SVZ framework).

Summary and conclusions
We have re-extracted the decay constants f D,D s and f B,B s and the running quark masses m c,b (m c,b ) using QCD spectral sum rules (QSSR). We have used as inputs, the recent values of the QCD (non-)perturbative parameters given in Table 1 and (for the first time) the renormalization group invariant quark and spontaneous masses in Eqs. (18) and (19). The results given in Eqs. (36), (39), (43) and (44) agree and improve existing QSSR results in the literature. Along the analysis, we have noticed that the values of the decay constants are very sensitive to the heavy quark mass and decrease when the heavy quark masses increase. Here we have used (for the first time) the scale independent Renormalization Group Invariant (RGI) heavy quark masses in the analysis. We have translated the on-shell mass expressions of the PT spectral function known to N2LO into the MS one where (as has been already noticed in previous works [19]) the PT series converge faster. We have also remarked that f P and m Q are affected by the choice of the continuum threshold t c which gives the largest errors. Here, like in our previous works [13,14,19,25], we have taken the conservative range of t c -values where the τ-or n-stability starts until the one where ones starts to have t c -stability. We have also seen that the subtraction point µ affects the truncated results within the OPE which has been the sources of apparent discrepancies and large errors of the results in the literature. Here, we have considered carefully the results at each subtraction point and deduced, from  these "QSSR data", the final results which should be independent on this arbitrary choice. In view of previous comments, we consider our results as improvements of the most recent ones to N2LO and using MDA in [25][26][27].
The results on f D and f D s agree within the errors with the data compiled in [2,62], while the upper bound on f D s can already exclude some existing data and theoretical estimates 12 .
As one can see in Table 8, our results are comparable (in values and precisions) with recent lattice simulations including dynamical quarks [63][64][65][66] 13 . These agreements are not surprising as both methods start from the same observables (the two-point correlator though evaluated in two different space-times) and use the 1st principles of QCD (here is the OPE in terms of the quarks and gluon condensates which semi-approximate the confinement regime). These agreements also confirm the accuracy of the MDA for describing the spectral function in the absence of a complete data, which has been tested earlier [5,6] and in this paper from the charmonium and bottomium systems.