R-Matrix Calculations for Few-Quark Bound States

The R--matrix method is implemented to study the heavy charm and bottom diquark, triquark, tetraquark and pentaquarks in configuration space, as the bound states of quark--antiquark, diquark--quark, diquark--antidiquark and diquark--antitriquark systems, respectively. The mass spectrum and the size of these systems are calculated for different partial wave channels. The calculated masses are compared with recent theoretical results obtained by other methods in momentum and configuration spaces and also by available experimental data.


Introduction
The original idea of the R-matrix theory was introduced by Kapur and Peierls [1], to remove unsatisfactory reliance on perturbation theory in nuclear reactions. It was a few years later that Wigner simplified the idea to the formulation of R-matrix theory in which all expressions are energy dependent [2][3][4]. Although the R-matrix theory was originally developed for treatment of nuclear resonances [4,5], it can also be used to describe all types of reaction phenomena and can be considered as an elegant method to solve the Schrödinger equation. These extensions especially became possible after the work of Bloch [6], by introducing a singular operator between internal and external regions.
The solution of the Schrödinger equation for the bound states of few-quarks in momentum space is numerically difficult to handle, because the confining part of the potential leads to singularity at small momenta. To overcome this problem we have successfully used a regularized form of the quark-antiquark [7] and recently diquark-antidiquark (DD) [8,9] interactions to solve the Lippmann-Schwinger (LS) equation in momentum space and calculate the mass spectrum of heavy quarkonia and tetraquarks. To this aim one can keep the divergent part of the potential fixed after exceeding a certain distance, called the regularization cutoff. This procedure creates an artificial barrier and the influence of tunneling barrier is manifested by significant changes in the energy eigenvalues at small distances.
Recently we have shown that the homogeneous LS equation can be formulated in configuration space to study the heavy tetraquarks as a bound state of DD system [10]. The variational methods are also used to study the few-quark bound states [11][12][13][14][15][16]. In order to obtain the variational energy, one must minimize the lowest eigenvalue with respect to the variational parameters after diagonalizing the Hamiltonian matrix. The variational energy arXiv:1609.03445v2 [hep-ph] 30 Sep 2016 can be obtained by differentiating the lowest eigenvalue with respect to the variational parameters. The successful application of the R-matrix theory to describe the resonance and scattering states resulting from the interaction of particles or systems of particles, which can be nucleons, nuclei, electrons, atoms and molecules [17,18], motivated us to implement it to the few-quark bound states. In this paper we have shown that R-matrix theory is an effective and efficient method to study the bound states of few-quarks by solving the Schrödinger equation in different partial wave channels. The successful implementation of the R-matrix method to heavy few-quark systems, paves the path to accurately predict the masses of few-quark bound states composed of light quarks.

R-matrix method for few-quark bound states in configuration space
In our study for diquark, triquark, tetraquark and pentaquark systems, we have used a twobody picture by considering them as the bound states of quark-antiquark (qq), diquarkquark (Dq), diquark-antidiquark (DD) and diquark-antitriquark (DT ) systems, respectively. The nonrelativistic bound state of any of these two-body systems with the pair relative distance r in a partial wave representation can be described by the Schrödinger equation. For simplicity of the notation we use AB for representation of these systems, where A and B stand for any of the subsystems. Since the AB interaction V (r) is a central force, thus the wave function of few-quark bound states consists of some type of radial function times a spherical harmonic function, i.e. Ψ nlm (r) = R nl (r)Y lm (θ , φ ). The radial function can be obtained by solution of the differential equation where E = m AB −m A −m B is AB binding energy (m AB is the mass of AB system composed of two subsystems A and B with masses m A and m B ). µ = m A m B m A + m B is the reduced mass and l is the orbital angular momentum of relative motion of AB system. In R-matrix calculations, the radial wave function is considered in two internal (r < r c ) and external (r > r c ) regions as R int and R ext , correspondingly. The parameter a is large enough to be sure that the mass spectrum is independent of it. The internal wave function R int is defined as combination of the basis functions u i (r) where In our study, the basis states parameters ρ i are Since in the external region the interaction between diquark and antidiquark is fixed, the external wave function can be considered as modified spherical Bessel function of the second kind R ext = C l k l (κr), where C l is a constant parameter and κ = 2µ (−E +V 0 ). Continuity of the internal and external wave functions and derivatives implies that The Hamiltonian H l is not Hermitian in the internal region. To avoid this the Bloch operator is defined as The dimensionless parameter b is an arbitrary real constant. The delta function indicates that the Bloch operator is a surface operator acting only on r = r c . The operator H l + L(b) is Hermitian, and therefore has a discrete spectrum in the finite region. Using Eq. (6), the Schrödinger equation in the internal region can be approximated by It means the logarithmic derivative of the wave function is continuous at r = r c . By multiplying Eq. (8) with u i and integrating in the internal region, we obtain the following equation to determine the unknown coefficients c i where Solving Eq. (9) for c i and substituting them into Eq.
where R(E, b) is the R-matrix given by The wave function in the internal region is then given by , the right hand side of Eq. (9) will be zero and consequently leads to the following Schrödinger-Bloch equation Since b depends on κ or the AB binding energy E that we want to calculate, we have written b ≡ b(E). The equation (14) can be written schematically as eigenvalue equation where the matrix elements of A and B matrices can be obtained as Since the A matrix is energy dependent, the solution of the eigenvalue equation (15) can be started by an initial guess for the energy E and the search in the binding energy can be stopped when λ − E E ≤ 10 −10 . In order to solve the Eq. (14), we have discretized the continuous variable r with Gauss-Legendre points using a hyperbolic-linear mapping [8]. To 15] GeV −1 using 75, 75 and 50 nodes in each subinterval, respectively. It indicates that the parameter r c , which divides the configuration space into internal and external regions, is chosen to be r c = 15 GeV −1 and we have numerically verified that the calculated masses of tetraquarks are independent of the regularization cutoff a. By having the AB binding energy and eigenvector c from the solution of eigenvalue Eq. (15), we can calculate the AB internal wave function by equations (10), (12) and (13). Of course, one can calculate the internal wave function directly using Eq. (2). Using AB wave function, we can evaluate the expectation value r for AB pair distance as where the AB radial wave function is normalized to 1, i.e. ∞ 0 dr r 2 R 2 (r) = 1.

Results and Discussion
For numerical solution of the integral equation (1) for qq, Dq and DD we have used the spin-independent interaction with the linear confining and the Coulomb-like one-gluon exchange potential F A and F B are the form factors of the subsystems A and B, correspondingly, and have the following functional form The parameters of this model are fixed from the analysis of heavy quarkonia masses and radiative decays [19][20][21].

Heavy quarkonia
For this first test of application of R-matrix method, we have solved the integral equation (1) to calculate the mass spectra of heavy quarkonia, mesons consisting heavy quark and antiquark. We have used the linear confining plus coulomb potential of Eq. (18) Table 1, our numerical results for masses of charmonium and bottomonium, obtained by R-matrix method are in excellent agreement with solution of Lippmann-Schwinger integral equation in momentum [7] and configuration [22] spaces and also with the experimental data [23].

Heavy baryons
In the next step we have calculated the masses of the ground state heavy baryons consisting of two light (u; d; s) and one heavy (c; b) quarks in the heavy-quark-light-diquark approximation. The used diquark mass and form factor parameters are given in Table 2. As we have shown in Table 3, our numerical results for different heavy baryons, calculated by nonrelativistic R-matrix method, are in good agreement with relativistic and spin-dependent results of EFG [24] and also with MLW results [25] of lattice nonrelativistic QCD. The relative difference between our and Ebert, Faustov, and Galkin (EFG) group results is less than 1.2 (4.8) % for bottom (charm) baryons. Clearly the difference comes from the relativistic effects and also spin terms of the potential that we have ignored in our calculations.

Heavy tetraquarks
In our calculations for heavy tetraquarks we have used the masses of diquark (antidiquark) and form factor parameters of Ref. [27] which are given in Table 2.
Our numerical results for the masses of charm (cqcq and cscs) and bottom (bqbq and bsbs) tetraquarks for s−, p− and d−wave channels with total spin S = 0 are listed in Tables  4 and 5. The tetraquark masses are calculated for scalar SS and axial-vector AĀ diquarkantidiquark contents. We have also calculated the expectation value of the relative distance between DD pair which can provide an estimate of the size of the tetraquarks. We have compared our results for tetraquark masses with recent results obtained in momentum and configuration spaces by solution of the nonrelativistic Lippmann-Schwinger integral equation [8,10] and also with those of previous relativistic studies by EFG reported in Refs. [27][28][29].
It indicates that the calculated masses for the ground state of charm tetraquarks (i.e. cqcq and cscs) with corresponding results from LS (in configuration space) [10] shows that they are in excellent agreement. Our results are also in good agreement with those of LS (in momentum space) and EFG with a relative percentage difference estimated to be at most 2.5 % and 4 %, respectively. While our R-matrix calculations and LS results are both done in a nonrelativistic spin-independent scheme, there is some difference between the results obtained in momentum and configuration spaces. Clearly, the difference between our Rmatrix results in configuration space and EFG results in momentum space is larger, and comes from the relativistic effects and also spin contribution in the DD interaction which appears in spin-orbit, spin-spin and tensor spin-space terms [28]. As we have shown in Ref. [8] the relativistic effect leads to a small reduction in the mass of heavy tetraquarks and decreases the masses of charm and bottom tetraquarks by less than 2 % and 0.2 %, respectively, whereas the spin contribution may lead to small decrese or increase in the masses of tetraquarks.

fm
In Fig. 1 we have shown few examples of the probability density |Ψ nlm (r)| 2 of bqbq tetraquark in SS state for s, p and d channels. As we can see the tetraquark probability functions for higher states have been expanded to larger distances which leads to larger expectation value for the relative distance between DD pair. As we can see in Tables 4 and  5, the expectation value of the relative distance between DD pair is almost the same for scalar SS and axial-vector AĀ diquark-antidiquark contents. It is larger for higher states and its size changes roughly with a factor of 3 from 1s to 3p state. In Table 6, we have compared our results for the masses of charm and bottom tetraquarks with the possible experimental candidates. They are in good agreement with a relative difference below 3.4 %.

Pentaquarks
By successful application of the R-matrix method for diquark, triquarks and tetraquarks we have also implemented it to study the pentaquarks as bound states of diquark-antitriquark systems (see Fig. 2).
In this study we have used two models of DT interaction. In the following sections we present the models and our numerical results.

One-pion exchange potential
One-pion exchange potential (OPEP) acting between a nucleon and a heavy meson (D or B) given as [40] V π (r) = where and the labels are: nucleon isospin I N , heavy-meson isospin I H (total isospin I = I N + I H ), tensor force S 12 , tensor potential V T (r), nucleon spin S N , light quark in heavy meson spin S l (sum of nucleon spin and light quark spin K = S N + S l ) and central potential V c (r). The parameters of OPEP potential are given in Table 7.  Our results for the masses of the pentaquarks composed of a nucleon and a meson (B and D mesons) for l = 0, S = 3 2 and J P = 3 2 + channel are given in Table 8. As we can see our results with R-matrix method are in excellent agreement with Cohen et al. results [40].

Cornell Potential
For the second test of Pentaquark calculations, the nonrelativistic linear-plus-Coulomb Cornell potential is used. It has the following form [41] V (r) = V c (r) +V spin−spin +V spin−orbit +V tensor (24) where and spin-spin, spin-orbit and tensor operators can be calculated as , J = L + 1 In Table 9, the parameter of Cornell potential used in pentaquark calculations are given.
Our results for the masses of charmoniumlike pentaquark P + c in s− and p−wave channels are given in Table 10. In our calculations we have ignored tensor force and we have considered central, spin-spin and spin-orbit terms of diquark-antitriquark interaction of Eq. (24). Beside the small difference between our results for the masses of s− and p−wave charmoniumlike pentaquarks and Lebed's results [42], which is about 5% and comes from neglected tensor force, the diquark-antitriquark separation r for s−wave channel with value of 0.306 fm is almost half of the obtained separation by Lebed with value of 0.64 fm. For the p−wave channel, neglecting the tensor force leads to the relative difference of 14% in calculated separations.
In conclusion, we have implemented the R-matrix method to calculate the mass spectra of heavy quarkonia, baryons, tetraquarks and pentaquarks in the two-body picture. Our numerical results for the masses of heavy charm and bottom few-quarks even by neglecting the relativistic effects, are in good agreement with other theoretical predictions and also with available experimental data.