# Charmonium Properties Using the Discrete Variable Representation (DVR) Method

Authors: Bhaghyesh A.
Journal: Advances in High Energy Physics
 (2021)

---

## Abstract

The Schrödinger equation is solved numerically for charmonium using the discrete variable representation (DVR) method. The Hamiltonian matrix is constructed and diagonalized to obtain the eigenvalues and eigenfunctions. Using these eigenvalues and eigenfunctions, spectra and various decay widths are calculated. The obtained results are in good agreement with other numerical methods and with experiments.


---

## Body

## 1. Introduction

The first quarkonium state was discovered independently at SLAC [1] and BNL [2], confirming the existence of heavy quark bound states. Since then, quarkonium has always been of great interest to particle physicists, being one of the extensively investigated system both theoretically and experimentally [3, 4]. New states are continuously being detected at various experiments. Recently, the LHCb collaboration [5] has detected a new state X3842 which is interpreted as a candidate for the unobserved ψ313D3 state. Both charmonium and bottomonium have rich spectrum of states below the open flavor threshold which have been experimentally observed, and various decay widths of these states have also been measured [6]. Studies on heavy quark systems are important because it gives information about quark interaction potential, confinement, QCD coupling constant, CKM matrix elements, and various other inputs to the standard model, some of which cannot be directly obtained from QCD.Theoretically, quarkonium systems have been studied by various formalisms based on phenomenological potential models [7–11], effective field theory [12], lattice gauge theory [13–16], Bethe Salpeter equation [17–20], etc. Among these, owing to its simplicity, formalism based on potential models is the widely chosen method to investigate quarkonium systems. In this method, both relativistic and quantum corrections can be easily incorporated. Potential models have been highly successful in predicting the spectra and decay widths [7, 9, 10]. In potential models, the usual method is to extract the properties of quarkonium by solving the Schrödinger equation using a chosen quark-antiquark potential. The widely used quark-antiquark potential in phenomenological models is the so-called Cornell potential [21–25], which includes a short-range Coulomb term and a linear confinement term. The form of this potential is also confirmed by lattice QCD calculations [26, 27].The Schrödinger equation for most of theqq¯ potentials (including the Cornell potential) cannot be solved analytically; hence, numerical solutions are called for. Some of the methods found in literature for solving the Schrödinger equation for qq¯ systems are numerical methods based on Runge-Kutte approximation [28, 29], Numerov matrix method [30–32], asymptotic iteration method [33–35], Fourier grid Hamiltonian method [36], variational method [37, 38], etc. Another method for numerically solving the Schrödinger equation is the discrete variable representation (DVR) method. This method has not been applied to quarkonium spectroscopy. Hence, in this article, we numerically solve the Schrödinger equation for cc¯ system using the discrete variable representation (DVR) scheme of Colbert and Miller [39]. DVR method was initially introduced by Harris et al. [40] and was extensively developed by Light and coworkers [41–46]. DVRs provide highly efficient and accurate solutions to quantum dynamical problems and have been widely used in atomic physics and quantum chemistry [47–56]. More details on DVR methods can be found in refs. [57, 58].This paper is organised as follows: a brief discussion on the potential model used to describe thecc¯ system and the DVR scheme used to solve the Schrödinger equation are presented in Section 2. The various decay properties calculated in the present analysis are given in Section 3. Results and discussions of the present work are given in Section 4.


## 2. Formalism

As a minimal model describing charmonium, we have used a nonrelativistic potential model, with the Hamiltonian(1)H=M+p22μ+Vr,wherep is the relative momentum, μ=mcmc¯/mc+mc¯ is the reduced mass of the cc¯ system, M=mc+mc¯, and Vr is the quark-antiquark potential. mc and mc¯ are the masses of individual quark and antiquark, respectively. For describing the quark-antiquark interaction, we use the standard Cornell potential plus a Gaussian-smeared contact hyperfine interaction [7]:
(2)Vr=−43αcr+br+V0+32παc9mc2σπ3e−σ2r2S⟶c·S⟶c¯.Parameters used in eq. (2) are given in Table 1 and are obtained by fitting the spectrum. Charmonium properties can be obtained by solving the Schrödinger equation corresponding to the Hamiltonian given in eq. (1) with potential given in eq. (2). In this work, to solve the Schrödinger equation, we have used the DVR scheme of Colbert and Miller [39]. In the DVR, the Hamiltonian is represented by a matrix on a uniform grid of points (ri=iΔr,i=1,2,3,..) in the coordinate space. Once the H-matrix is constructed, diagonalization gives us the bound state eigenvalues and the amplitudes of eigenfunctions on the grid point chosen.Table 1
Parameters used in the model.
ParameterValueαc0.54b0.136GeV2V00.149 GeVσ1.1 GeVmc1.4 GeVIn ref. [39], authors have shown that the kinetic energy matrix can be written as
(3)Tij=ℏ22mΔr2−1i−jπ23−12i2i=j,2i−j2−2i+j2i≠j,with ri=iΔr,i=1,2,⋯, where Δr is the grid spacing. The potential energy matrix is diagonal
(4)Vij=Vriδij.We have used eqs. (1)–(4) to construct the Hamiltonian matrix in the present model, which upon diagonalization returns the bound state eigenvalues and the amplitudes of eigenfunctions on the chosen grid points. In the present analysis, we have chosen a grid of length 10 fm with 1000 grid points. The eigenvalue problem for the matrix of the Hamiltonian (1) was solved using Mathematica. For a given eigenvalue, in order to obtain the eigenfunction in the entire range of coordinates, we have used the built-in interpolation function in Mathematica through the obtained eigenfunctions on grid points. This interpolation function was used as the representation of the reduced radial wavefunction for our further analysis. Obtained wavefunctions for 3S1 and 3PJ states are shown in Figure 1. In order to compute fine structure of the L≠0 states, we add the spin-orbit and tensor terms perturbatively [7]:
(5)Vfs=1mc22αcr3−b2rL⟶·S⟶+4αcr3T.Figure 1
The wavefunctions of charmonium3S1 states (a) and 3PJ states (b). Solid curve represents the n=1 state while dashed curve represents the n=2 state.
(a)(b)The computed mass spectra of charmonium are listed in Table2. Using the obtained wavefunction, we also compute the root mean square radii (<r2>) and the square of the radial wavefunction at the origin (R02) for these states, and our results are listed in Table 3.Table 2
Mass spectrum of charmonium states (in MeV).
StatePresentExp [6][7][10][59][60][61][62]J/ψ30973096.900±0.00630903097309430943090.03091.7ηc(1S)29862983.9±0.529822979299529892981.62992.4ψ(2S)36633686.097±0.02536723673364936813671.83671.4ηc(2S)36203637.6±1.236303623360636023630.33631.7ψ(3S)40554039±140724022403641294071.64075.5ηc(3S)402540433991400040584043.24048.1ψ(4S)43854421±444064273436245144406.14415.0ηc(4S)436043844250432844484383.74393.3ψ(5S)46784463465448634703.8ηc(5S)46574446462247994685.0ψ(6S)49474608492551854976.9ηc(6S)49294595489351244960.4χc2(1P)35463556.17±0.0735563554355634803549.03548.1χc1(1P)34983510.67±0.0535053510352334683505.43501.8χc0(1P)34193414.71±0.3034243433345734283424.53425.8hc(1P)35073525.38±0.1135163519353434703515.63510.5χc2(2P)395539723937395639553964.83970.0χc1(2P)390939253901392539383924.93925.8χc0(2P)383938523842386638973852.33856.7hc(2P)391839343908393639433933.63933.4ψ31D33379138063799380137553805.33800.6ψ21D32378638003798380537723800.43796.7ψ13D137713773.13±0.3537853787379937753785.03783.1ηc211D2378537993796380237653799.43795.1ψ323D3414641674103415141764165.54167.1ψ223D2413841584100415241884158.24160.2ψ23D141224191±541424089414541884141.54145.1ηc221D2413841584099415041824157.64159.1Table 3
Average radii (in fm) and square of the radial wave function at the origin (inGeV3).
State<r2>[62][61][10]R02[62][61]ηc(1S)0.3800.3750.36551.6491.54051.2294ηc(2S)0.8630.8390.83280.7310.75410.8717ηc(3S)1.2501.2101.20720.5730.60880.683ηc(4S)1.5841.5311.53060.5020.54300.5994ηc(5S)1.8851.82250.4610.5503J/ψ0.4340.4210.41430.410.9761.18611.97675ψ(2S)0.8970.8670.86270.910.8970.70920.7225ψ(3S)1.2741.2301.22871.381.2740.59140.6006ψ(4S)1.6031.5471.54781.871.6030.53400.5417ψ(5S)1.9021.83702.391.9020.505381P10.7000.6780.6738≈00≈01P30.7120.6890.71730.71≈00≈02P11.1081.0711.0697≈00≈02P31.1201.0821.11071.19≈00≈01D10.9310.8990.8984≈00≈01D30.9320.9010.91790.96≈00≈02D11.3041.2581.2595≈00≈02D31.3051.2611.19141.44≈00≈0


## 3. Decay Properties

For quarkonium, most of the decay properties are dependent on the wave function. Hence, to test the wavefunctions as obtained in the previous section, we calculate leptonic decay widths and radiative decay widths (M1 and E1) of some charmonium states.


### 3.1. Leptonic Decay Widths

The leptonic decay widths of the vector states are calculated using the Van Royen-Weisskopf formula [63, 64].
(6)Γee=4α2ec2MnSRnS021−16αs3π,where MnS is the mass for nS state, ec is the charm quark charge in unit of electron charge, α is the fine structure constant, αs≈αs2mc is the strong coupling constant, and Rns0 is the radial nS wave function at the origin. The terms in parenthesis are the QCD radiative correction factor. Obtained results are listed in Table 4.Table 4
Leptonic decay widths (in keV).
StatePresentExp [6][59][60][65][10]J/ψ4.9795.55±0.14±0.023.6232.9251.85326.60ψ(2S)2.1372.33±0.041.0851.5330.59832.40ψ(3S)1.4600.86±0.070.7481.0910.38121.42ψ(4S)1.1310.58±0.070.5990.8560.28470.97ψ(5S)0.9300.5080.7070.22860.70


### 3.2.M1 Radiative Transitions

Magnetic dipole (M1) radiative transitions obey the selection rule ΔL=0 and ΔS=±1. The M1 widths are evaluated using the formula [7].
(7)ΓM1n2S+1LJ⟶n′2S′+1LJ′′=432J′+12L+1δLL′δS,S′+1ec2αmc2<ψf∣ψi>2Eγ3EfMi,where Eγ is the emitted photon energy, <ψf∣ψi> is the overlap integral involving initial and final radial wavefunctions, Ef is the total energy of the final state, and Mi is the mass of the initial state. Calculated M1 widths are listed in Table 5.Table 5
M1 radiative partial widths (in keV).
TransitionPresentExp [6][7][60][59][66]1S31⟶1S102.661.58±0.372.92.7221.6472.392S31⟶2S100.170.21±0.150.211.1720.1350.192S31⟶1S105.021.00±0.154.67.50669.577.80


### 3.3.E1 Radiative Transitions

Electric dipole (E1) radiative transitions obey the selection rule ΔL=±1 and ΔS=0. The E1 widths are evaluated using the formula [7].
(8)ΓE1n2S+1LJ⟶n′2S′+1LJ′′=43CfiδSS′ec2α<ψfrψi>2Eγ3EfMi,where <ψfrψi> is the spatial matrix element involving the initial and final radial wavefunctions, and Cfi is the angular matrix element given by
(9)Cfi=maxL,L′2J′+1L′J′SJL12.E1 widths obtained from the present analysis are listed in Table 6.Table 6
E1 radiative partial widths (in keV).
TransitionPresentExp [6][7][60][59][66][10]2S31⟶1P3229.7827.99±0.963862.3127.0736432S31⟶1P3149.3128.67±1.055443.29210.3945622S31⟶1P3049.3828.78±0.986321.86311.9327741P32⟶1S31436.45374.30±19.73424157.225233.853274731P31⟶1S31319.08288.12±16.09314146.317189.862693541P30⟶1S31175.78151.2±9.99152112.030118.291411672D31⟶1P326.07<17.44.95.7226.455.45.82D31⟶1P31159.0567.73±6.7312593.775139.521151502D31⟶1P30425.59187.68±17.72403161.504343.87243486


## 3.1. Leptonic Decay Widths

The leptonic decay widths of the vector states are calculated using the Van Royen-Weisskopf formula [63, 64].
(6)Γee=4α2ec2MnSRnS021−16αs3π,where MnS is the mass for nS state, ec is the charm quark charge in unit of electron charge, α is the fine structure constant, αs≈αs2mc is the strong coupling constant, and Rns0 is the radial nS wave function at the origin. The terms in parenthesis are the QCD radiative correction factor. Obtained results are listed in Table 4.Table 4
Leptonic decay widths (in keV).
StatePresentExp [6][59][60][65][10]J/ψ4.9795.55±0.14±0.023.6232.9251.85326.60ψ(2S)2.1372.33±0.041.0851.5330.59832.40ψ(3S)1.4600.86±0.070.7481.0910.38121.42ψ(4S)1.1310.58±0.070.5990.8560.28470.97ψ(5S)0.9300.5080.7070.22860.70


## 3.2.M1 Radiative Transitions

Magnetic dipole (M1) radiative transitions obey the selection rule ΔL=0 and ΔS=±1. The M1 widths are evaluated using the formula [7].
(7)ΓM1n2S+1LJ⟶n′2S′+1LJ′′=432J′+12L+1δLL′δS,S′+1ec2αmc2<ψf∣ψi>2Eγ3EfMi,where Eγ is the emitted photon energy, <ψf∣ψi> is the overlap integral involving initial and final radial wavefunctions, Ef is the total energy of the final state, and Mi is the mass of the initial state. Calculated M1 widths are listed in Table 5.Table 5
M1 radiative partial widths (in keV).
TransitionPresentExp [6][7][60][59][66]1S31⟶1S102.661.58±0.372.92.7221.6472.392S31⟶2S100.170.21±0.150.211.1720.1350.192S31⟶1S105.021.00±0.154.67.50669.577.80


## 3.3.E1 Radiative Transitions

Electric dipole (E1) radiative transitions obey the selection rule ΔL=±1 and ΔS=0. The E1 widths are evaluated using the formula [7].
(8)ΓE1n2S+1LJ⟶n′2S′+1LJ′′=43CfiδSS′ec2α<ψfrψi>2Eγ3EfMi,where <ψfrψi> is the spatial matrix element involving the initial and final radial wavefunctions, and Cfi is the angular matrix element given by
(9)Cfi=maxL,L′2J′+1L′J′SJL12.E1 widths obtained from the present analysis are listed in Table 6.Table 6
E1 radiative partial widths (in keV).
TransitionPresentExp [6][7][60][59][66][10]2S31⟶1P3229.7827.99±0.963862.3127.0736432S31⟶1P3149.3128.67±1.055443.29210.3945622S31⟶1P3049.3828.78±0.986321.86311.9327741P32⟶1S31436.45374.30±19.73424157.225233.853274731P31⟶1S31319.08288.12±16.09314146.317189.862693541P30⟶1S31175.78151.2±9.99152112.030118.291411672D31⟶1P326.07<17.44.95.7226.455.45.82D31⟶1P31159.0567.73±6.7312593.775139.521151502D31⟶1P30425.59187.68±17.72403161.504343.87243486


## 4. Discussion and Summary

In the present work, we have numerically solved the Schrödinger equation for charmonium system using the DVR scheme of Colbert and Miller. The Hamiltonian matrix was constructed and diagonalised to obtain the masses and wavefunctions of charmonium states. In Table2, we compare the masses of radially and orbitally excited cc¯ states with experiment [6] and other theoretical predictions [7, 10, 59–62]. Authors in refs. [7, 59–62] have also used Cornell type potential to study the cc¯ system, where as in ref. [10], authors use a screened potential. From Table 2, we see that the predictions using the DVR method are in good agreement with experiment and other theoretical predictions. In Table 3, we have compared our predictions for the root mean square radii (<r2>) and the square of the radial wavefunction at the origin (R02) with other theoretical predictions [10, 61, 62]. The values of radial wavefunctions at the origin are important inputs for calculating quarkonium production cross-sections [22] and various decay amplitudes. We present our results for leptonic decays in Table 4 in comparison with experiment and other models. Our predictions for lower states are in good agreement with the experimental results. For higher excited states, our predictions are higher than the experimental results. We present results of E1 and M1 radiative transitions in Tables 5 and 6, respectively. Radiative transitions in quarkonia are important because they are one of the few mechanisms that produce transitions among qq¯ states with different quantum numbers. This decay mechanism also helps to produce excited P-wave states and F-wave states which are otherwise difficult to achieve [7]. M1 decays in particular allow to access spin-singlet states. From Tables 5 and 6, we see that there is a wide range of precautions for the radiative decay widths even though all these models [7, 59, 60, 66] employ a Cornell type potential. This may be due to the difference in wavefunctions of charmonium states used in these models. Our predictions for radiative decays are in accordance with experiment and other theoretical predictions. Inclusion of higher multipole contributions, coupled channel effects, relativistic corrections, etc., would give a better fit to the experimental results.In summary, in this article, we have successfully employed the DVR method to investigate the spectra and decays of charmonium. The obtained results of the present study are in good agreement with experimental data and with other theoretical models.



---

*Source: 9991152-2021-07-30.xml*