(\[eqn:eigenvalue\_equation\]) can be solved in a basis function approach [@Vary:2009gt; @Li:2015zda]. The basis function approach is particularly advantageous for the present model with the holographic confining potential, since, in the absence of the one-gluon exchange term, it can be diagonalized analytically. On the other hand, the confining interactions in momentum space are highly singular. The solutions can be expressed in terms of the analytic functions $\phi_{nm}$ and $\chi_l$. For the transverse direction, we have (see Fig. \[fig:basis\_functions\_a\]): $$\phi_{nm}(\vec q_\perp; b) = b^{-1} \sqrt{\frac{4\pi n!}{(n+|m|)!}} \bigg(\frac{q_\perp}{b}\bigg)^{|m|}
\exp\big(-q^2_\perp/(2b^2)\big) L_n^{|m|}(q^2_\perp/b^2) \exp\big(\imag m \theta_q),$$ where $\vec q_\perp \triangleq \vec k_\perp/\sqrt{x(1-x)}$, $q_\perp = |\vec q_\perp|$, $\theta_q = \arg \vec q_\perp$. $L_n^a(z)$ is the associated Laguerre polynomial. $b$ is the harmonic oscillator (HO) basis parameter in mass dimension. Following Ref. [@Li:2015zda], we choose $b\equiv \kappa$ to match the confining strength. For simplicity, we will often omit the label $b$ though it is implicit throughout. In the longitudinal direction, we have (see Fig. \[fig:basis\_functions\_b\]): $$\chi_l(x; \alpha,\beta) =
\sqrt{4\pi(2l+\alpha+\beta+1)}\sqrt{\frac{\Gamma(l+1)\Gamma(l+\alpha+\beta+1)}{\Gamma(l+\alpha+1)\Gamma(l+\beta+1)}}
x^{{\mathsmaller{\frac{\beta}{2}}}}(1-x)^{{\mathsmaller{\frac{\alpha}{2}}}} P^{(\alpha,\beta)}_l(2x-1).$$ Here $P_l^{(\alpha, \beta)}(z)$ is the Jacobi polynomial. $\alpha$ and $\beta$ are dimensionless basis parameters. In the model, they are $\alpha = 2m_{\bar q}(m_q+m_{\bar q})/\kappa^2$, $\beta = 2m_{q}(m_q+m_{\bar q})/\kappa^2$. Again, we will drop the explicit dependence on $\alpha$ 