1 Introduction

Ever since the discovery of first charmonium state (J/Ψ) at SLAC and BNL, quarkonia spectroscopy has always been a subject of continuous research, both experimentally and theoretically [1, 2]. Thanks to the experimental facilities at various colliders, new states have been constantly discovered [2, 3]. Recently discovered conventional quarkonium states include ηb(1S), ηb(2S), hb(1P), hb(2P), χb1(3P) and χb2(3P) [3]. There are also lot of states that are discovered (named X, Y, Z) that do not fit into the conventional quarkonia spectra and are assumed to be strong candidates for non-conventional hadronic states [2,3,4,5,6]. Comprehensive reviews on the status of heavy quarkonium physics can be found in References [5,6,7,8,9,10,11,12,13].

Even though quantum chromodynamics (QCD) is the theory of strong interactions, there is no direct method to extract the properties of observed hadrons directly from QCD Lagrangian. Lattice QCD provides great hope in this direction and significant results are being reported [1, 2]. Various other approaches widely used in the literature to study hadronic systems include QCD sum rules [14,15,16,17,18,19], nonrelativistic and relativistic potential models [20,21,22,23,24,25,26,27,28,29,30, 32,33,34,35,36], effective field theory [37,38,39] and formalisms based on Bethe-Salpeter equation and Dyson-Schwinger equation [40,41,42,43,44,45,46,47,48,49]. Among theses QCD inspired potential models have been largely successful in explaining hadronic properties. Importance of potential models stems from the fact that masses (upto rather high excitations) and decay properties can be investigated in an unified way.

One of the important quantities that can be obtained through potential models are the vector and pseudoscalar decay constants of heavy mesons. Decay constants are important because it can be used to estimate the hadronic matrix elements and provide information about the short distance structure of hadrons and non-perturbative QCD dynamics [19, 27]. Precise knowledge of their values play very important roles, for example in determination of CKM matrix elements, study of weak decays, meson mixing, etc. [44, 50]. They are input parameters which constrains various decays and processes and their precise knowledge is important in various experimental measurements [51, 52]. Vector (fV) and pseudoscalar (fP) decay constants are defined through the matrix elements [27]:

$$ \begin{array}{@{}rcl@{}} m_{V}f_{V}\epsilon^{\mu}&=&\left<0|\bar{{\varPsi}}\gamma^{\mu}{\varPsi}|V\right>\\ p^{\mu} f_{P}&=&i\left<0|\bar{{\varPsi}}\gamma^{\mu}\gamma^{5}{\varPsi}|P\right> \end{array} $$

The quark model definition of decay constants are through [25, 27]:

$$ \begin{array}{@{}rcl@{}} f_{V}&=&\sqrt{\frac{3}{M}}\int \frac{d^{3}k}{(2\pi)^{3}}\sqrt{1+\frac{m_{Q}}{E_{k}}} \sqrt{1+\frac{m_{\bar{Q}}}{E_{\bar{k}}}}\left( 1+\frac{k^{2}}{3(E_{k}+m_{Q})(E_{\bar{k}}+m_{\bar{Q}})}\right) \phi(\mathbf{k}) \end{array} $$
(1)
$$ \begin{array}{@{}rcl@{}} f_{P}&=&\sqrt{\frac{3}{M}}\int \frac{d^{3}k}{(2\pi)^{3}}\sqrt{1+\frac{m_{Q}}{E_{k}}} \sqrt{1+\frac{m_{\bar{Q}}}{E_{\bar{k}}}}\left( 1-\frac{k^{2}}{(E_{k}+m_{Q})(E_{\bar{k}}+m_{\bar{Q}})}\right) \phi(\mathbf{k}) \end{array} $$
(2)

The above expressions are mostly used in potential models which use the momentum space representation. In the non-relativistic limit \((k\rightarrow 0)\), one obtains the well-known Van Royen-Weisskopf formula which relates the decay constants to the value of the meson wave-function at the origin [53]:

$$ f_{P/V}^{2}=12 \frac{|\psi_{P/V}(0)|^{2}}{M_{P/V}} $$
(3)

In (3), ψP/V(0) is the wave function of the pseudoscalar/vector meson calculated at the origin (r = 0) and MP/V is the mass of the pseudoscalar/vector meson. This formula is widely used in potential model calculations [22,23,24,25, 27, 29, 32,33,34,35]. One reason that it is widely used is that the coordinate space wave function is sufficient to calculate the decay constants. In potential models, in which the vector and pseudoscalar states differ only by the spin-dependent corrections, one obtains fPfV [22,23,24,25, 32, 35, 54]. The relativistic effects which enters as the v2/c2 corrections in hyperfine spltting cannot be neglected in quarkonia. They produce significant difference in the wavefunction at origin for pseudoscalar and vector mesons [55], and hence fP becomes very different from fV. Ebert et al. [25] have shown that inclusion of relativistic effects produces considerable (\(\sim \)70 MeV) difference between vector and pseudoscalar decay constants. Also the present experimental estimates (fJ/ψ = 416 ± 6 MeV [3] and \(f_{\eta _{c}}=335\pm 75\) MeV [56]), lattice [19, 52] and other theoretical studies [19, 44, 45, 48, 50, 57, 58] predict that fV is greater than fP. Hence the non-relativistic Van Royen-Weisskopf formula fails to provide exact results in accordance to experimental and lattice results. We suggest that the higher order terms (\(\mathcal {O} (k^{2}/c^{2})\)) in (1) and (2) are important. In this paper, we derive expressions for the decay constants of S-wave quarkonia by approximating the quark model definition by including terms upto k2/c2.

This paper is organised as follows: in Section 2 we derive expressions for pseudoscalar and vector decay constants for S-wave quarkonia including the k2/c2 corrections. In Section 3 we discuss the theoretical model used in the present analysis. Results and discussions of our present analysis are given in Section 4.

2 Decay Constants

For quarkonia \(m_{Q}=m_{\bar {Q}}=m\) and hence \(E_{k}=E_{\bar {k}}=E\). Also using the relation, k2 = (E + m)(Em), (1) becomes

$$ \begin{array}{@{}rcl@{}} f_{V}&=&\sqrt{\frac{3}{M}}\int\frac{d^{3}k}{(2\pi)^{3}}\left( 1+\frac{m}{E}\right)\left( 1+\frac{(E+m)(E-m)}{3(E+m)(E+m)}\right) \phi(\mathbf{k})\\ &=&\frac{1}{3}\sqrt{\frac{3}{M}}\int\frac{d^{3}k}{(2\pi)^{3}}\left( 4+\frac{2m}{E}\right)\phi(\mathbf{k}) \end{array} $$
(4)

Now,

$$ \begin{array}{@{}rcl@{}} E=(k^{2}+m^{2})^{1/2}&=&\left[m^{2}\left( 1+\frac{k^{2}}{m^{2}}\right)\right]^{1/2}\\ &=& m\left( 1+\frac{k^{2}}{2m^{2}}\right)+\mathcal{O}\left( \frac{k^{4}}{m^{4}}\right) \end{array} $$
(5)

Substituting (5) in (4), we get

$$ \begin{array}{@{}rcl@{}} f_{V}&=&\frac{1}{3}\sqrt{\frac{3}{M}}\int\frac{d^{3}k}{(2\pi)^{3}}\left( 4+\frac{2m}{m(1+\frac{k^{2}}{2m^{2}})}\right)\phi(\mathbf{k})\\ &=&\frac{1}{3}\sqrt{\frac{3}{M}}\int\frac{d^{3}k}{(2\pi)^{3}}\left( 4+\frac{4m^{2}}{k^{2}+2m^{2}}\right)\phi(\mathbf{k}) \end{array} $$
(6)

Now using the three dimensional Fourier transform of ϕ(k) the second integral in (6) can be written as:

$$ \begin{array}{@{}rcl@{}} \int\frac{d^{3}k}{(2\pi)^{3}}\left( \frac{4m^{2}}{k^{2}+2m^{2}}\right)\phi(\mathbf{k})&=&4m^{2}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{k^{2}+2m^{2}}\int d^{3}r ~ e^{-i\mathbf{k}.\vec{r}}\tilde{\phi}(\vec{r})\\ &=&4m^{2} \int d^{3}r ~\tilde{\phi}(\vec{r})\int\frac{d^{3}k}{(2\pi)^{3}}\frac{e^{-i\mathbf{k}.\vec{r}}}{k^{2}+2m^{2}}\\ &=&4m^{2} \int d^{3}r ~\tilde{\phi}(\vec{r})\int\frac{d^{3}k}{(2\pi)^{3}}\frac{e^{-i\mathbf{k}.\vec{r}}}{k^{2}+(\sqrt{2}m)^{2}}\\ &=&4m^{2} \int d^{3}r ~\tilde{\phi}(\vec{r})\cdot \frac{1}{4\pi r} e^{-\sqrt{2}mr}\\ &=&\frac{m^{2}}{\pi}\int d^{3}r ~\frac{\tilde{\phi}(\vec{r})}{r} e^{-\sqrt{2}mr} \end{array} $$
(7)

For S-wave states, there is no angular dependence:

$$ \tilde{\phi}(\vec{r})=\frac{R(r)}{\sqrt{4\pi}} $$
(8)

Therefore, (7) becomes:

$$ \begin{array}{@{}rcl@{}} \int\frac{d^{3}k}{(2\pi)^{3}}\left( \frac{4m^{2}}{k^{2}+2m^{2}}\right)\phi(\mathbf{k})&=&\frac{m^{2}}{\pi}\frac{1}{\sqrt{4\pi}}\int d^{3}r~ \frac{R(r)}{r} e^{-\sqrt{2}mr}\\ &=&\frac{m^{2}}{\pi}\frac{1}{\sqrt{4\pi}} 4\pi \int dr~ \frac{R(r)}{r} r^{2} e^{-\sqrt{2}mr}\\ &=& \sqrt{\frac{4}{\pi}} m^{2} \int dr~R(r) r e^{-\sqrt{2}mr} \end{array} $$
(9)

Similarly the integral

$$ \int\frac{d^{3}k}{(2\pi)^{3}}\phi(\mathbf{k})=\tilde{\phi}(\vec{r}=0)=\frac{R(0)}{\sqrt{4\pi}} $$
(10)

Substituting (9) and (10) in (6), we get the final expression for vector decay constant of S-wave quarkonia as:

$$ f_{V}=\frac{1}{3}\sqrt{\frac{3}{M}}\sqrt{\frac{4}{\pi}}\left[R(0)+m^{2}\int dr~ R(r) r e^{-\sqrt{2}mr}\right] $$
(11)

Similarly setting \(m_{Q}=m_{\bar {Q}}=m\), \(E_{k}=E_{\bar {k}}=E\) and k2 = (E + m)(Em), (2) becomes

$$ \begin{array}{@{}rcl@{}} f_{P}&=&\sqrt{\frac{3}{M}}\int\frac{d^{3}k}{(2\pi)^{3}}\left( 1+\frac{m}{E}\right)\left( 1-\frac{(E+m)(E-m)}{(E+m)(E+m)}\right) \phi(\mathbf{k})\\ &=&\sqrt{\frac{3}{M}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{2m}{E}\phi(\mathbf{k}) \end{array} $$
(12)

Substituting (5) in (12), we get:

$$ \begin{array}{@{}rcl@{}} f_{P}&=&\sqrt{\frac{3}{M}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{4m^{2}}{k^{2}+2m^{2}}\phi(\mathbf{k}) \end{array} $$

Fourier tranforming ϕ(k) and performing the integrals as done earlier, we obtain the final expression for pseudoscalar decay constant of S-wave quarkonia as:

$$ f_{P}=\sqrt{\frac{3}{M}}\sqrt{\frac{4}{\pi}} m^{2}\int dr~ R(r) r e^{-\sqrt{2}mr} $$
(13)

Taking into account the QCD corrections [59,60,61], the decay constants of vector and pseudoscalar states equal to

$$ \begin{array}{@{}rcl@{}} \bar{f}_{V/P}&=&f_{V/P}\left( 1-\delta_{V/P}\frac{\alpha_{S}}{\pi}\right), \end{array} $$
(14)

where δV = 8/3 and δP = 2.

3 Theoretical Model

In order to obtain the bound state wavefunctions of heavy mesons, we have used a screened potential model suggested by Li and Chao [30, 31]. The quark-antiquark potential in this model is given by

$$ V(r)=-\frac{4}{3}\frac{\alpha_{C}}{r}+\lambda\left( \frac{1-e^{-\mu r}}{\mu}\right)+V_{0} $$
(15)

The hyperfine interaction term which provide the spin singlet-triplet splitting is

$$ H_{hf}=\frac{32 \pi \alpha_{C}}{9 {m_{Q}^{2}}} \delta_{\sigma} (r)~ \vec{S}_{Q}\cdot\vec{S}_{\bar{Q}}~, $$
(16)

where \(\delta _{\sigma } (r)=(\sigma /\sqrt {\pi })^{3} e^{-\sigma ^{2}r^{2}}\). The hyperfine interaction is treated perturbatively. The parameters used in (15) and (16) are same as in References [30, 31] and are listed in Table 1 for completeness.

Table 1 Model Parameters

Using the potential (15), we solve the corresponding non-relativistic Schrodinger equation numerically [62] and obtain the bound state wavefunctions for heavy mesons. Substituting the wavefunctions in (14) we obtain respectively the vector and pseudoscalar decay constants.

To check the validity of our numerical computation method, we calculate the mass spectra and the radial wavefunction at the origin (|R(0)|2) and compare the results with that of References [30, 31]. The results are presented in Tables 2 and 3 along with the spin-averaged masses (MSA) and hyperfine contributions. \(\langle H_{hf}^{V} \rangle \) and \(\langle H_{hf}^{P} \rangle \) are respectively the hyperfine contributions for spin-triplet (vector) and spin-singlet (pseudoscalar) states. The |R(0)|2 of References [30, 31] are extracted from the corresponding leptonic decay widths. The mass spectra and |R(0)|2 obtained in the present analysis are in good agreement with the corresponding results of References [30, 31], confirming the validity of our numerical computation method.

Table 2 Charmonium spectrum (in GeV) and radial wavefunctions at the origin (in GeV3)
Table 3 Bottomonium spectrum (in GeV) and radial wavefunctions at the origin (in GeV3)

4 Results and Discussions

The computed values of pseudoscalar and vector decay constants are presented in Tables 4 and 5. fNR and \(\bar {f}^{NR}\) are respectively the decay constants calculated using the non-relativistic Van Royen-Weisskopf formula (3) without and with the QCD correction factor. f and \(\bar {f}\) are respectively the decay constants calculated from our present analysis ((11) and (13)) without and with the QCD correction factor (14).

Table 4 Decay constants (in MeV) of charmonium states in various schemes
Table 5 Decay constants (in MeV) of bottomonium states in various schemes

In Tables 678 and 9, pseudoscalar (\(\bar {f}_{P}\)) and vector (\(\bar {f}_{V}\)) decay constants for charmonium and bottomonium states obtained from our analysis are compared with the available experimental results [3, 56] and predictions from lattice [19, 52, 63,64,65,66,67,68], nonrelativistic potential models [27, 69], QCD sum rules (QCDSR) [19, 57, 58], light front quark model (LFQM) [50] and formalisms based on Bethe-Salpeter and Dyson-Schwinger equations (BS/DSE) [44,45,46,47,48,49].

Table 6 Pseudoscalar decay constants (\(\bar {f}_{P}\)) of charmonium states (in MeV)
Table 7 Vector decay constants (\(\bar {f}_{V}\)) of charmonium states (in MeV)
Table 8 Pseudoscalar decay constants (\(\bar {f}_{P}\)) of bottomonium states (in MeV)
Table 9 Vector decay constants (\(\bar {f}_{V}\)) of bottomonium states (in MeV)

Vector decay constants are experimentally extracted from the leptonic decay widths through the relation

$$ {\varGamma} (V\rightarrow e^{+}e^{-})=\frac{4\pi}{3}\alpha^{2}{e_{Q}^{2}}\frac{{f_{V}^{2}}}{M_{V}} $$

Substituting the PDG averages for \({\varGamma } (V\rightarrow e^{+}e^{-})\) in the above expression, one can estimate the vector decay constants. As far as the pseudoscalar decay constants are concerned, only the \(f_{n_{c}}\) is experimentally measured through the \(B\rightarrow n_{c} K\) decay [56]. \(f_{n_{b}}\) is not yet experimentally determined.

As seen from Tables 69, the predictions of decay constants from various models show a wide range of variation. For most of the states, the experimental results are yet to be available. Results from our analysis are found to be in good agreement with the available experimental results. For both vector and pseudoscalar states, the decay constants keep decreasing with the radial excitation number which is also in accordance with existing experimental/lattice results and other theoretical studies.

The ratios of pseudoscalar to vector decay constant (\(\bar {f}_{P}/\bar {f}_{V}\))) for the 1S states are shown in Table 10. For both \(c\overline c\) and \( b\overline b\) mesons, we obtain \(\bar {f}_{P}/\bar {f}_{V} < 1\), which is consistent with the available experimental data (\(f_{\eta _{c}}/f_{J/\psi }=0.81\pm 0.18\) [3, 56]) and other theoretical predictions. This result is in sharp contrast with predictions of some of the nonrelativistic potential models [22,23,24,25, 32, 35] which predict fP/fV ≈ 1. This indicate that for quarkonia the relativistic effects are important in calculating the decay constants.

Table 10 Ratios of pseudoscalar to vector decay constants (\(\bar {f}_{P}/\bar {f}_{V}\))

In summary, we have derived expressions for pseudoscalar and vector decay constants for S-wave quarkonia including the k2/c2 corrections. Our predictions are consistent with the experiment and predictions from other theoretical studies.