Quarkonium as relativistic bound state on the light front

TL;DR

This paper models charmonium and bottomonium as relativistic bound states using light-front holography, incorporating a longitudinal confinement and one-gluon exchange to improve spectroscopy and compute various properties.

Contribution

It introduces a novel light-front holographic model with longitudinal confinement and running coupling for heavy quarkonia, enhancing spectral accuracy and wave function analysis.

Findings

Improved mass spectra for charmonium and bottomonium

Accurate decay constants and radii calculations

Detailed light-front wave functions and distributions

Abstract

We study charmonium and bottomonium as relativistic bound states in a light-front quantized Hamiltonian formalism. The effective Hamiltonian is based on light-front holography. We use a recently proposed longitudinal confinement to complete the soft-wall holographic potential for the heavy flavors. The spin structure is generated from the one-gluon exchange interaction with a running coupling. The adoption of asymptotic freedom improves the spectroscopy compared with previous light-front results. Within this model, we compute the mass spectroscopy, decay constants and the r.m.s. radii. We also present a detailed study of the obtained light-front wave functions and use the wave functions to compute the light-cone distributions, specifically the distribution amplitudes and parton distribution functions. Overall, our model provides a reasonable description of the heavy quarkonia.

Tables

Table 1. Table 1: Summary of the model parameters (see text).

	$N_f$	${\alpha }_s\left(0\right)$	${\mu }_{\text{g}}$ (GeV)	$\kappa$ (GeV)	$m_q$ (GeV)	rms (MeV)	$\overline{{\delta }_{j}M}$ (MeV)	$N_{\text{exp}}$	$N_{\max }=L_{\max }$
$c\overline{c}$	4	0.6	0.02	0.966	1.603	31	17	8	32
$b\overline{b}$	5	0.6	0.02	1.389	4.902	38	8	14	32

Table 2. Table 2: Comparison of differences between fits and PDG experimental data between results of Ref.  Li:2015zda and those presented here. δ ​ M c ​ c ¯ 𝛿 subscript 𝑀 𝑐 ¯ 𝑐 \delta M_{c\bar{c}} is the rms mass deviation for charmonium from the PDG data. δ j ​ M ¯ c ​ c ¯ subscript ¯ subscript 𝛿 𝑗 𝑀 𝑐 ¯ 𝑐 \overline{\delta_{j}M}_{c\bar{c}} is the mean mass spread for charmonium. “fix- α s subscript 𝛼 𝑠 \alpha_{s} (refitted)” improves the bottomonium fits by ∼ similar-to \sim 10 MeV.

	${\overline{{\delta }_{j}M}}_{c\overline{c}}$	$\delta M_{c\overline{c}}$ (rms)	${\overline{{\delta }_{j}M}}_{b\overline{b}}$	$\delta M_{b\overline{b}}$ (rms)	$N_{\max }=L_{\max }$
fix-${\alpha }_s$ Li:2015zda	49 MeV	52 MeV	17 MeV	*58 MeV	24
fix-${\alpha }_s$ (refitted)	—	—	15 MeV	48 MeV	24
running-${\alpha }_s$	17 MeV	31 MeV	7 MeV	39 MeV	24
running-${\alpha }_s$	17 MeV	31 MeV	8 MeV	38 MeV	32

Table 3. Table 3: Model sensitivity with respect to the basis size N max = L max subscript 𝑁 subscript 𝐿 N_{\max}=L_{\max} . The model parameters fits and the r.m.s. deviations are well converged as N max = L max subscript 𝑁 subscript 𝐿 N_{\max}=L_{\max} increases.

	${\alpha }_s\left(0\right)$	$N_f$	${\mu }_{\text{g}}$ (GeV)	$\kappa$ (GeV)	$m_q$ (GeV)	rms (MeV)	$\overline{{\delta }_{j}M}$ (MeV)	$N_{\text{exp}}$	$N_{\max }=L_{\max }$
$c\overline{c}$	0.6	4	0.02	0.985	1.570	41	15	8 states	8
				0.979	1.587	32	21		16
				0.972	1.596	31	17		24
				0.966	1.603	31	17		32
$b\overline{b}$	0.6	5	0.02	1.387	4.894	48	6	14 states	8
				1.392	4.899	41	6		16
				1.390	4.901	39	7		24
				1.389	4.902	38	8		32

Table 4. Table 4: Comparison of heavy quarkonia moments from NRQCD Bodwin:2006dn , QCD sum rule Braguta:2006wr ; Braguta:2007fh ; Braguta:2007tq , light-front quark model Choi:2007ze and DSE Ding:2015rkn . The DSE results are obtained at 2 GeV. Results from other approaches are evaluated at quark mass scale μ ≃ m q similar-to-or-equals 𝜇 subscript 𝑚 𝑞 \mu\simeq m_{q} . The BLFQ results are given at N max = L max = 8 subscript 𝑁 subscript 𝐿 8 N_{\max}=L_{\max}=8 for charmonium and N max = L max = 32 subscript 𝑁 subscript 𝐿 32 N_{\max}=L_{\max}=32 for bottomonium, roughly corresponding to UV cutoffs μ = κ ​ N max ≈ 1.7 ​ m q 𝜇 𝜅 subscript 𝑁 1.7 subscript 𝑚 𝑞 \mu=\kappa\sqrt{N_{\max}}\approx 1.7m_{q} . For the convenience of comparison with other approaches, we also provide the extrapolated (ext.) or interpolated (int.) results at the effective quark mass scale (BLFQ*). The 3 ​ σ 3 𝜎 3\sigma ( ∼ 99.75 % similar-to absent percent 99.75 \sim 99.75\% C.L.) statistical errors (prediction intervals) are included to indicate the quality of the extrapolation or interpolation procedure. The pQCD asymptotics Lepage:1980fj ⟨ ξ n ⟩ asy = 3 / ( n + 1 ) ​ ( n + 2 ) subscript delimited-⟨⟩ superscript 𝜉 𝑛 asy 3 𝑛 1 𝑛 2 \langle\xi^{n}\rangle_{\mathrm{asy}}=3/(n+1)(n+2) and the AdS/QCD results of Brodsky and de Téramond (AdS/QCD, Brodsky:2014yha ) ⟨ ξ n ⟩ lfh = 2 ​ ( n − 1 ) !! / ( n + 2 ) !! subscript delimited-⟨⟩ superscript 𝜉 𝑛 lfh 2 double-factorial 𝑛 1 double-factorial 𝑛 2 \langle\xi^{n}\rangle_{\textsc{lfh}}=2(n-1)!!/(n+2)!! , and the IMA modified AdS/QCD results (IMA, Vega:2009zb ; Swarnkar:2015osa ) are also provided for comparison.

		NRQCD	QCDSR	LFQM	DSE	BLFQ*	BLFQ	AdS/QCD	IMA	pQCD
${\eta }_c$	$⟨{\xi }^2⟩$	0.075(11)	0.070(7)	0.0084${}_{-0.007}{}^{}{}_{}{}^{+0.004}$	0.10	0.096(13)	0.12	0.25	0.0058	0.20
	$⟨{\xi }^4⟩$	0.010(3)	0.012(2)	0.017${}_{-0.003}{}^{}{}_{}{}^{+0.001}$	0.032	0.019(2)	0.036	0.13	0.0084	0.086
	$⟨{\xi }^6⟩$	0.0017(7)	0.0032(9)	0.0047${}_{-0.0010}{}^{}{}_{}{}^{+0.0006}$	0.015	0.0036(27)	0.014	0.078	0.0018	0.047
	$⟨{\xi }^8⟩$				0.0059	$-$0.0005(46)	0.0068	0.055	0.00047	0.030
	$\mu$	$m_c$	$m_c$	$m_c$	2 GeV	$m_c$ (ext.)	$1.7m_c$			$\infty$
$J/\psi$	$⟨{\xi }^2⟩$	0.075(11)	0.070(7)	0.082${}_{-0.006}{}^{}{}_{}{}^{+0.004}$	0.039	0.096(20)	0.11	0.25	0.0058	0.20
	$⟨{\xi }^4⟩$	0.010(3)	0.012(2)	0.016${}_{-0.002}{}^{}{}_{}{}^{+0.002}$	0.0038	0.021(9)	0.030	0.13	0.0084	0.086
	$⟨{\xi }^6⟩$	0.0017(7)	0.0031(8)	0.0046${}_{-0.0010}{}^{}{}_{}{}^{+0.0005}$	7.3$\times {10}^{-4}$	0.0060(41)	0.011	0.078	0.0018	0.047
	$⟨{\xi }^8⟩$				3.3$\times {10}^{-4}$	0.0015(15)	0.0053	0.055	0.00047	0.030
	$\mu$	$m_c$	$m_c$	$m_c$	2 GeV	$m_c$ (ext.)	$1.7m_c$			$\infty$
${\eta }_c^{\prime }$	$⟨{\xi }^2⟩$	0.22(14)	0.18${}_{-0.07}{}^{}{}_{}{}^{+0.005}$			0.157(9)	0.179
	$⟨{\xi }^4⟩$	0.085(110)	0.051${}_{-0.031}{}^{}{}_{}{}^{+0.031}$			0.043(7)	0.059
	$⟨{\xi }^6⟩$	0.039(77)	0.017${}_{-0.014}{}^{}{}_{}{}^{+0.016}$			0.013(3)	0.025
	$⟨{\xi }^6⟩$					0.0036(5)	0.012
	$\mu$	$m_c$	$m_c$	$m_c$		$m_c$ (ext.)	$1.7m_c$
${\eta }_b$	$⟨{\xi }^2⟩$				0.070	0.052(2)	0.071	0.25		0.20
	$⟨{\xi }^4⟩$				0.015	0.0081(61)	0.015	0.13		0.086
	$⟨{\xi }^6⟩$				0.0042	0.0020(48)	0.0051	0.078		0.047
	$⟨{\xi }^8⟩$				0.0013	0.0006(31)	0.0021	0.055		0.030
	$\mu$	$m_b$	$m_b$	$m_b$	2 GeV	$m_b$ (int.)	$1.6m_b$			$\infty$
${\rm Y}$	$⟨{\xi }^2⟩$				0.014	0.047(17)	0.061	0.25		0.20
	$⟨{\xi }^4⟩$				4.3$\times {10}^{-4}$	0.0066(73)	0.012	0.13		0.086
	$⟨{\xi }^6⟩$				$4.4\times {10}^{-5}$	0.0014(63)	0.0036	0.078		0.047
	$⟨{\xi }^8⟩$				$3.7\times {10}^{-6}$	0.0004(30)	0.0014	0.055		0.030
	$\mu$	$m_b$	$m_b$	$m_b$	2 GeV	$m_b$ (int.)	$1.6m_b$			$\infty$
${\eta }_b^{\prime }$	$⟨{\xi }^2⟩$					0.082(13)	0.10
	$⟨{\xi }^4⟩$					0.013(15)	0.022
	$⟨{\xi }^6⟩$					0.003(10)	0.0068
	$⟨{\xi }^8⟩$					0.0007(44)	0.0027
	$\mu$	$m_b$	$m_b$	$m_b$	2 GeV	$m_b$ (int.)	$1.6m_b$

Table 5. Table 5: Spinor matrix elements u ¯ s 1 ′ ​ ( p 1 ′ ) ​ γ μ ​ u s 1 ​ ( p 1 ) ​ v ¯ s 2 ​ ( p 2 ) ​ γ μ ​ v s 2 ′ ​ ( p 2 ′ ) subscript ¯ 𝑢 superscript subscript 𝑠 1 ′ superscript subscript 𝑝 1 ′ subscript 𝛾 𝜇 subscript 𝑢 subscript 𝑠 1 subscript 𝑝 1 subscript ¯ 𝑣 subscript 𝑠 2 subscript 𝑝 2 superscript 𝛾 𝜇 subscript 𝑣 superscript subscript 𝑠 2 ′ superscript subscript 𝑝 2 ′ \bar{u}_{s_{1}^{\prime}}(p_{1}^{\prime})\gamma_{\mu}u_{s_{1}}(p_{1})\bar{v}_{s_{2}}(p_{2})\gamma^{\mu}v_{s_{2}^{\prime}}(p_{2}^{\prime}) . m q subscript 𝑚 𝑞 m_{q} ( m a subscript 𝑚 𝑎 m_{a} ) is the mass of the quark (antiquark). x = p 1 + / P + 𝑥 subscript superscript 𝑝 1 superscript 𝑃 x=p^{+}_{1}/P^{+} and x ′ = p 1 ′ ⁣ + / P + superscript 𝑥 ′ subscript superscript 𝑝 ′ 1 superscript 𝑃 x^{\prime}=p^{\prime+}_{1}/P^{+} are longitudinal momentum fractions of the quark, p → ⟂ = p → 1 ⟂ − x ​ P → ⟂ subscript → 𝑝 perpendicular-to subscript → 𝑝 perpendicular-to 1 absent 𝑥 subscript → 𝑃 perpendicular-to \vec{p}_{\perp}=\vec{p}_{1\perp}-x\vec{P}_{\perp} and p → ⟂ ′ = p → 1 ⟂ ′ − x ′ ​ P → ⟂ subscript superscript → 𝑝 ′ perpendicular-to subscript superscript → 𝑝 ′ perpendicular-to 1 absent superscript 𝑥 ′ subscript → 𝑃 perpendicular-to \vec{p}^{\prime}_{\perp}=\vec{p}^{\prime}_{1\perp}-x^{\prime}\vec{P}_{\perp} are relative transverse momenta. For convenience, we use the complex representation for the transverse vectors, viz, p ≜ p x + i ​ p y ≜ 𝑝 subscript 𝑝 𝑥 i subscript 𝑝 𝑦 p\triangleq p_{x}+{\mathrm{i}}p_{y} and p ∗ ≜ p x − i ​ p y ≜ superscript 𝑝 subscript 𝑝 𝑥 i subscript 𝑝 𝑦 p^{*}\triangleq p_{x}-{\mathrm{i}}p_{y} .

$s_1$	$s_2$	$s_1^{\prime }$	$s_2^{\prime }$	$\frac{{\overline{u}}_{s_1^{\prime }}\left(p_1^{\prime }\right){\gamma }_{\mu }{u}_{s_1}\left(p_1\right){\overline{v}}_{s_2}\left(p_2\right){\gamma }^{\mu }{v}_{s_2^{\prime }}\left(p_2^{\prime }\right)}{2\sqrt{x\left(1-x\right)x^{\prime }\left(1-x^{\prime }\right)}}$
$+$	$+$	$+$	$+$	$m_q^2\frac{1}{x{x}^{\prime }}+m_a^2\frac{1}{\left(1-x\right)\left(1-x^{\prime }\right)}+\frac{p{p}^{\prime \ast }}{x\left(1-x\right)x^{\prime }\left(1-x^{\prime }\right)}$
$-$	$-$	$-$	$-$	$m_q^2\frac{1}{x{x}^{\prime }}+m_a^2\frac{1}{\left(1-x\right)\left(1-x^{\prime }\right)}+\frac{p^{\ast }{p}^{\prime }}{x\left(1-x\right)x^{\prime }\left(1-x^{\prime }\right)}$
$+$	$-$	$+$	$-$	$m_q^2\frac{1}{x{x}^{\prime }}+m_a^2\frac{1}{\left(1-x\right)\left(1-x^{\prime }\right)}+\left(\frac{p^{\prime \ast }}{x^{\prime }}+\frac{p^{\ast }}{1-x}\right)\left(\frac{p}{x}+\frac{p^{\prime }}{1-x^{\prime }}\right)$
$-$	$+$	$-$	$+$	$m_q^2\frac{1}{x{x}^{\prime }}+m_a^2\frac{1}{\left(1-x\right)\left(1-x^{\prime }\right)}+\left(\frac{p^{\ast }}{x}+\frac{p^{\prime \ast }}{1-x^{\prime }}\right)\left(\frac{p^{\prime }}{x^{\prime }}+\frac{p}{1-x}\right)$
$+$	$+$	$+$	$-$	$m_a\frac{x^{\prime }}{\left(1-x\right)\left(1-x^{\prime }\right)}\left(\frac{p^{\prime }}{x^{\prime }}-\frac{p}{x}\right)$
$-$	$-$	$-$	$+$	$m_a\frac{x^{\prime }}{\left(1-x\right)\left(1-x^{\prime }\right)}\left(\frac{p^{\ast }}{x}-\frac{p^{\prime \ast }}{x^{\prime }}\right)$
$-$	$+$	$-$	$-$	$m_a\frac{x}{\left(1-x\right)\left(1-x^{\prime }\right)}\left(\frac{p^{\prime }}{x^{\prime }}-\frac{p}{x}\right)$
$+$	$-$	$+$	$+$	$m_a\frac{x}{\left(1-x\right)\left(1-x^{\prime }\right)}\left(\frac{p^{\ast }}{x}-\frac{p^{\prime \ast }}{x^{\prime }}\right)$
$+$	$+$	$-$	$+$	$m_q\frac{1-x^{\prime }}{x{x}^{\prime }}\left(\frac{p}{1-x}-\frac{p^{\prime }}{1-x^{\prime }}\right)$
$-$	$-$	$+$	$-$	$m_q\frac{1-x^{\prime }}{x{x}^{\prime }}\left(\frac{p^{\prime \ast }}{1-x^{\prime }}-\frac{p^{\ast }}{1-x}\right)$
$+$	$-$	$-$	$-$	$m_q\frac{1-x}{x{x}^{\prime }}\left(\frac{p}{1-x}-\frac{p^{\prime }}{1-x^{\prime }}\right)$
$-$	$+$	$+$	$+$	$m_q\frac{1-x}{x{x}^{\prime }}\left(\frac{p^{\prime \ast }}{1-x^{\prime }}-\frac{p^{\ast }}{1-x}\right)$
$+$	$-$	$-$	$+$	$-m_q{m}_a\frac{{\left(x-x^{\prime }\right)}^2}{x\left(1-x\right)x^{\prime }\left(1-x^{\prime }\right)}$
$-$	$+$	$+$	$-$
$+$	$+$	$-$	$-$	0
$-$	$-$	$+$	$+$

Equations

$$H_{\mathrm{eff}}=\frac{\vec{k}^{2}_{\perp}+m_{q}^{2}}{x}+\frac{\vec{k}^{2}_{\perp}+m_{\bar{q}}^{2}}{1-x}+\kappa^{4}\vec{\zeta}_{\perp}^{2}-\frac{\kappa^{4}}{(m_{q}+m_{\bar{q}})^{2}}\partial_{x}\big{(}x(1-x)\partial_{x}\big{)}\\ -\frac{C_{F}4\pi\alpha_{s}(Q^{2})}{Q^{2}}\bar{u}_{s^{\prime}}(k^{\prime})\gamma_{\mu}u_{s}(k)\bar{v}_{\bar{s}}(\bar{k})\gamma^{\mu}v_{\bar{s}^{\prime}}(\bar{k}^{\prime}).$$

$$Q^{2}=\frac{1}{2}\Big{(}\sqrt{\frac{x^{\prime}}{x}}\vec{k}_{\perp}-\sqrt{\frac{x}{x^{\prime}}}\vec{k}^{\prime}_{\perp}\Big{)}^{2}+\frac{1}{2}\Big{(}\sqrt{\frac{1-x^{\prime}}{1-x}}\vec{k}_{\perp}-\sqrt{\frac{1-x}{1-x^{\prime}}}\vec{k}^{\prime}_{\perp}\Big{)}^{2}+\frac{1}{2}(x-x^{\prime})^{2}\Big{(}\frac{m_{q}^{2}}{xx^{\prime}}+\frac{m_{\bar{q}}^{2}}{(1-x)(1-x^{\prime})}\Big{)}+\mu_{g}^{2}.$$

$$\alpha_{s}(Q^{2})=\frac{1}{\beta_{0}\ln\big{(}Q^{2}/\Lambda^{2}+\tau\big{)}}\triangleq\frac{\alpha_{s}(M_{\textsc{z}}^{2})}{1+\alpha_{s}(M_{\textsc{z}}^{2})\beta_{0}\ln(\mu^{2}_{\textsc{ir}}+Q^{2})/(\mu^{2}_{\textsc{ir}}+M^{2}_{\textsc{z}})},$$

$$H_{\text{eff}}|{\psi }_h(P,j,m_j)⟩=M_h^2|{\psi }_h(P,j,m_j)⟩.$$

$$|\psi_{h}(P,j,m_{j})\rangle=\sum_{s,\bar{s}}\int_{0}^{1}\frac{{\mathrm{d}}x}{2x(1-x)}\int\frac{{\mathrm{d}}^{2}k_{\perp}}{(2\pi)^{3}}\,\psi^{(m_{j})}_{s\bar{s}/h}(\vec{k}_{\perp},x)\\ \times\frac{1}{\sqrt{N_{c}}}\sum_{i=1}^{N_{c}}b^{\dagger}_{s{}i}\big{(}xP^{+},\vec{k}_{\perp}+x\vec{P}_{\perp}\big{)}d^{\dagger}_{\bar{s}{}i}\big{(}(1-x)P^{+},-\vec{k}_{\perp}+(1-x)\vec{P}_{\perp}\big{)}|0\rangle.$$

$$\big{\{}b_{si}(p^{+},\vec{p}_{\perp}),b_{s^{\prime}i^{\prime}}^{\dagger}(p^{\prime+},\vec{p}^{\prime}_{\perp})\big{\}}=\big{\{}d_{si}(p^{+},\vec{p}_{\perp}),d_{s^{\prime}i^{\prime}}^{\dagger}(p^{\prime+},\vec{p}^{\prime}_{\perp})\big{\}}=2p^{+}(2\pi)^{3}\delta^{3}(p-p^{\prime})\delta_{ss^{\prime}}\delta_{ii^{\prime}},$$

$$⟨{\psi }_h(P,j,m_j)|{\psi }_{h^{\prime }}(P^{\prime },j^{\prime },m_j^{\prime })⟩=2P^{+}(2\pi {)}^3{\delta }^3(P-P^{\prime }){\delta }_{j{j}^{\prime }}{\delta }_{m_j,m_j^{\prime }}{\delta }_{h{h}^{\prime }},$$

$$\sum _{s,\bar{s}}{\int }_0^1\frac{\mathrm{d}x}{2x(1-x)}\int \frac{{\mathrm{d}}^2{k}_{\perp }}{(2\pi {)}^3}{\psi }_{s\bar{s}/h^{\prime }}^{(m_j^{\prime })\ast }({\vec{k}}_{\perp },x){\psi }_{s\bar{s}/h}^{(m_j)}({\vec{k}}_{\perp },x)={\delta }_{h{h}^{\prime }}{\delta }_{m_j,m_j^{\prime }}.$$

$${\tilde{\psi }}_{s\bar{s}}({\vec{r}}_{\perp },x)\equiv \frac{1}{\sqrt{x(1-x)}}\int \frac{{\mathrm{d}}^2{k}_{\perp }}{(2\pi {)}^2}{e}^{\mathrm{i}{\vec{k}}_{\perp }\cdot {\vec{r}}_{\perp }}{\psi }_{s\bar{s}}({\vec{k}}_{\perp },x).$$

$$\sum _{s,\bar{s}}{\int }_0^1\frac{\mathrm{d}x}{4\pi }\int {\mathrm{d}}^2{r}_{\perp }{\tilde{\psi }}_{s\bar{s}/h^{\prime }}^{(m_j^{\prime })\ast }({\vec{r}}_{\perp },x){\tilde{\psi }}_{s\bar{s}/h}^{(m_j)}({\vec{r}}_{\perp },x)={\delta }_{h{h}^{\prime }}{\delta }_{m_j,m_j^{\prime }}.$$

$${\hat{m}}_P|{\psi }_h(P,j,m_j)⟩=(-i{)}^{2j}\mathsf{P}|{\psi }_h(\tilde{P},j,-m_j)⟩,\hat{C}|{\psi }_h(P,j,m_j)⟩=\mathsf{C}|{\psi }_{\bar{h}}(P,j,m_j)⟩.$$

$${\vec{\mathcal{J}}}^2|{\psi }_h(P,j,m_j)⟩=j(j+1)|{\psi }_h(P,j,m_j)⟩,{\mathcal{J}}_z|{\psi }_h(P,j,m_j)⟩=m_j|{\psi }_h(P,j,m_j)⟩.$$

$$\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{(}{\mathrm{i}}m\theta_{q}),$$

$${\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}^{\frac{\beta }{2}}(1-x{)}^{\frac{\alpha }{2}}{P}_l^{(\alpha ,\beta )}(2x-1).$$

$${\psi }_{s{s}^{\prime }/h}({\vec{k}}_{\perp },x)=\sum _{n,m,l}{\psi }_h(n,m,l,s,s^{\prime }){\varphi }_{nm}({\vec{k}}_{\perp }/\sqrt{x(1-x)}){\chi }_l(x).$$

$${\tilde{\psi }}_{s{s}^{\prime }/h}({\vec{r}}_{\perp },x)=\sqrt{x(1-x)}\sum _{n,m,l}{\psi }_h(n,m,l,s,s^{\prime }){\tilde{\varphi }}_{nm}(\sqrt{x(1-x)}{\vec{r}}_{\perp }){\chi }_l(x).$$

$$\widetilde{\phi}_{nm}(\vec{\rho}_{\perp};b^{-1})=b\sqrt{\frac{n!}{\pi(n+|m|)!}}(b\rho_{\perp})^{|m|}\exp\big{(}-b^{2}\rho^{2}_{\perp}/2\big{)}L_{n}^{|m|}(b^{2}\rho_{\perp}^{2})\exp\big{[}{\mathrm{i}}m\theta_{\rho}+{\mathrm{i}}\pi(n+|m|/2)\big{]}.$$

$$2n+|m|+1\le N_{\max },0\le l\le L_{\max }.$$

$$(-i{)}^{2j}\mathsf{P}=$$

$$\mathsf{C}=$$

$$M_{n,m,l}^2=(m_q+m_{\bar{q}}{)}^2+2{\kappa }^2(2n+|m|+l+1)+\frac{{\kappa }^4}{(m_q+m_{\bar{q}}{)}^2}l(l+1).$$

$${\psi }_{nml}({\vec{k}}_{\perp },x)={\varphi }_{nm}({\vec{k}}_{\perp }/\sqrt{x(1-x)}){\chi }_l(x).$$

$$\psi_{\mathrm{gs}}(\vec{k}_{\perp},x)=N\exp\big{[}-\vec{k}^{2}_{\perp}/(2\kappa^{2}x(1-x))\big{]}\big{(}x(1-x)\big{)}^{2m^{2}_{q}/\kappa^{2}}.$$

$$N\exp\big{[}-\vec{k}^{2}_{\perp}/(2\kappa^{2}x(1-x))\big{]}\;\to\;N^{\prime}\exp\big{[}-(\vec{k}^{2}_{\perp}+m_{q}^{2})/(2\kappa^{2}x(1-x))\big{]}.$$

$$|\ell -s|\le j\le \ell +s,\mathsf{P}=(-1{)}^{\ell +1},\mathsf{C}=(-1{)}^{\ell +s},$$

$$\overline{M}\equiv \sqrt{\frac{M_{-j}^2+M_{1-j}^2+\cdots +M_j^2}{2j+1}},$$

$$\overline{\delta_{j}M}\equiv\sqrt{\frac{1}{N_{h}}\sum_{h}^{j\neq 0}(\delta_{j}M_{h})^{2}}.\qquad\Big{(}N_{h}\equiv\sum_{h}^{j\neq 0}1\Big{)}$$

$$⟨0|\overline{\psi }(0){\gamma }^{+}{\gamma }_5\psi (0)|P(p)⟩=$$

$$⟨0|\overline{\psi }(0){\gamma }^{+}\psi (0)|V(p,\lambda )⟩=$$

$$\frac{f_{P,V}}{2\sqrt{2N_c}}={\int }_0^1\frac{\mathrm{d}x}{2\sqrt{x(1-x)}}\int \frac{{\mathrm{d}}^2{k}_{\perp }}{(2\pi {)}^3}{\psi }_{\uparrow \downarrow \mp \downarrow \uparrow }^{(\lambda =0)}(x,{\vec{k}}_{\perp }).$$

$$\langle r^{2}_{\mathrm{c}}\rangle=-6\frac{\partial}{\partial Q^{2}}F_{\mathrm{ch}}(Q^{2})\Big{|}_{Q\to 0},\quad\langle r^{2}_{\mathrm{m}}\rangle=-6\frac{\partial}{\partial Q^{2}}F_{\mathrm{gr}}(Q^{2})\Big{|}_{Q\to 0}.$$

Full text

I Introduction

Non-perturbative calculations of quantum chromodynamics (QCD) provide insights into the fundamental structure of hadrons which constitute the majority of the visible matter in the Universe. Lattice gauge theory has produced high precision results for hadron spectroscopy and many other observables. It is expected that Lattice QCD will eventually provide a valid description of the experimental data arising from both the theoretical progress and the growth of computational capacity. On the other hand, QCD at high energy is most conveniently expressed through the light-front variables Lepage:1980fj . While the so-called “hard processes” may be evaluated through perturbation theory (pQCD), non-perturbative information from QCD is also needed and is encoded within the so-called “light-cone distributions”. The light-cone distributions are intrinsically Minkowskian, and cannot be easily extracted from a Euclidean formulation of quantum field theories. It is anticipated that the light-front111In this article, we use the words “light-front” and “light-cone” interchangeably. Hamiltonian formalism provides a complementary alternative to lattice gauge theory Bakker:2013cea , with convenient access to light-cone distributions and other observables.

In principle, the hadron mass spectrum and light-front wave functions (LFWFs) can be obtained from diagonalizing the light-front quantized QCD (LFQCD) Hamiltonian operator Brodsky:1997de . Ab initio light-front Hamiltonian approaches, such as Discretized Light-Cone Quantization (DLCQ, Pauli:1985ps ) and Basis Light-Front Quantization (BLFQ, Vary:2009gt ), have made important strides in tackling various test problems, and show promise of advancing towards more realistic field theories, including QCD Vary:2016emi . As a complementary method to these ab initio approaches, light-front holography constructs an effective Hamiltonian based on insights from string theory, and has been shown, notwithstanding criticisms (e.g., Ballon-Bayona:2014oma ; Glozman:2009bt ), to be a valuable approximation to QCD Brodsky:2014yha . The efforts to improve light-front holography can be roughly cast into two categories: one is on the holographic QCD side (see Brodsky:2014yha and the references therein); the other is on the light-front Hamiltonian side (see Hiller:2016itl for a recent review).

The present work falls into the second category. We generalize the light-front holographic QCD of Brodsky and de Téramond to incorporate quark masses and quarkonium spin structure by extending the “soft-wall” light-front Hamiltonian. Our model introduces a phenomenological effective Hamiltonian. Key elements include a confining potential in the longitudinal direction and an effective one-gluon exchange interaction derived from light-front QCD Li:2015zda ; Li:2016wwu . It was long pointed out by Lepage and Brodsky Lepage:1980fj that the dominant ultraviolet (UV) physics can be analyzed through one-gluon exchange. Here, we combine the one-gluon exchange physics at short distance and the holographic QCD at long distance. The present work improves our previous calculation Li:2015zda by including the evolution of the strong coupling as a function of invariant 4-momentum transfer. Incorporating the running coupling not only implements important QCD physics, but also improves the UV asymptotics of the kernel. In particular, a previous non-covariant UV counterterm is now removed and the hyperfine structure is readily improved as we present in this work.

The motivation of the present work is multi-fold. As stated, we supplement the light-front holographic QCD interaction with one-gluon exchange, rather than patching the holographic wave functions with, e.g., spin structures (see, e.g., Ref. Chen:2016dlk and the references therein). The spectroscopy and the wave functions are obtained as a natural output. More importantly, we solve the problem using the basis function method Vary:2009gt . Effectively, we are applying BLFQ to a phenomenological interaction that emulates features of QCD. Indeed, this work is a direct extension of the BLFQ approach to positronium in QED Wiecki:2014ola . Finally, we acknowledge the similarities between our work and the relativistic bound-state models in QCD (e.g., Refs. Sommerer:1994bk ; Spence:1995bm ; Spence:1999db ; Maris:1999nt ; Leitao:2016bqq ), especially the light-front QCD bound-state models Pauli:1996ne ; Brisudova:1995hv ; Brisudova:1996vw ; Glazek:2003ky ; Glazek:2006cu ; Choi:2015ywa .

We organize this paper as follows. In Sect. II, we introduce the theoretical model, including the longitudinal confinement and a running strong coupling. The formulation and the methods are detailed in Sect. III. Sect. IV summarizes and analyzes the numerical results, including the spectroscopy, decay constants and radii. Sect. V presents LFWFs and light-cone distributions computed from them. We summarize the paper in Sect. VI.

II Holographic Confinement plus One-Gluon Exchange

We extend light-front holography by introducing realistic QCD interactions such as the one-gluon exchange interaction with running coupling Li:2015zda . In addition we include finite quark masses, important for heavy flavors, as well as a longitudinal confining potential to complement the transverse holographic confining potential. Spin structure and excited states (radial and angular) naturally emerge from the one-gluon exchange and its non-perturbative interplay with the confining potential Li:2016wwu . The effective Hamiltonian $H_{\text{eff}}\equiv P^{+}P^{-}_{\text{eff}}-\vec{P}^{2}_{\perp}$ reads,

[TABLE]

where $\vec{\zeta}_{\perp}\equiv\sqrt{x(1-x)}\vec{r}_{\perp}$ is Brodsky and de Téramond’s holographic variable Brodsky:2014yha , $\partial_{x}f(x,\vec{\zeta}_{\perp})=\partial f(x,\vec{\zeta}_{\perp})/\partial x|_{\vec{\zeta}}$, $C_{F}=(N_{c}^{2}-1)/(2N_{c})=4/3$ is the color factor for the color singlet state. $\kappa$ is the strength of the confinement, and $m_{q}$ ($m_{\bar{q}}$) is the mass of the quark (anti-quark). $Q^{2}=-(1/2)(k^{\prime}-k)^{2}-(1/2)(\bar{k}^{\prime}-\bar{k})^{2}$ is the average 4-momentum squared carried by the exchanged gluon. In terms of kinematical variables,

[TABLE]

II.1 Longitudinal Confinement

In Eq. (1), the term $\kappa^{4}\vec{\zeta}^{2}_{\perp}\equiv\kappa^{4}x(1-x)\vec{r}^{2}_{\perp}$ is the “soft-wall” confinement from light-front holography, which is introduced in the massless case. For heavy quarkonium, the quark masses and the longitudinal dynamics cannot be ignored222Without a longitudinal confinement, the longitudinal excitations will not be separated by mass gaps. In light-front holography (no quark mass nor one-gluon exchange), these excitations are degenerate and the system is two-dimensional in nature. and we introduce a longitudinal confining interaction to complete the transverse holographic confinement. The form of the longitudinal confinement is designed to produce a power-law behavior for the distribution amplitudes $\phi(x)\sim x^{a}(1-x)^{b}$ at the endpoints (cf. Gutsche:2012ez ; Gutsche:2013zia ; Chabysheva:2012fe ).

We fix the strength of the longitudinal confinement by matching to the transverse holographic confinement in the non-relativistic limit. Therefore, rotational symmetry is retained in the heavy-quark limit. Another advantage of this choice for the longitudinal confinement is that it produces, without the one-gluon exchange, analytic solutions. Therefore, it affords computational convenience within the basis function method (see Sect. III.2). In the massless limit, our wave function (without the one-gluon exchange) reduces to the soft-wall wave function of Brodsky and de Téramond333See Eq. (8) for our normalization convention. Brodsky:2014yha . It has been suggested that in the massless limit one can choose the longitudinal confining strength to be independent of the quark mass to reproduce the Gell-Mann-Oakes-Renner relation Gutsche:2012ez . Our proposal shares some similarities with other proposals in the literature Gutsche:2014oua ; Trawinski:2014msa ; Chabysheva:2012fe .

II.2 Running Coupling

As mentioned, we employ a running coupling based on the 1-loop pQCD. The running coupling is a function of the 4-momentum transfer squared $Q^{2}=-q^{2}>0$ (see also Fig. 1), viz

[TABLE]

where $\beta_{0}=(33-2N_{f})/(12\pi)$, with $N_{f}$ the number of quark flavors, $N_{f}=4$ for charmonium and $N_{f}=5$ for bottomonium. A constant $\tau$ is introduced to avoid the pQCD IR catastrophe. Similar ansätze are widely adopted in the literature (e.g. Maris:1999nt ). $\Lambda$ and constant $\tau$ are obtained by fixing the strong coupling at the Z-boson mass $\alpha_{s}(M_{\textsc{z}}^{2})=0.1183$ and at $Q=0$. In practice, we choose $\alpha_{s}(0)=0.6$, corresponding to $\mu_{\textsc{ir}}=0.55\,\mathrm{GeV}$ for $N_{f}=4$. We find, however, the spectra are not sensitive to the choice of $\alpha_{s}(0)$ within the range of $0.4\leq\alpha_{s}(0)\leq 0.8$.

Introducing the evolution of the strong coupling implements asymptotic freedom for the one-gluon exchange through a natural dependence on the covariant 4-momentum transfer $Q^{2}$. The use of the running coupling also serves to improve the UV asymptotics of the one-gluon exchange kernel. In our previous work Li:2015zda , we used a fixed coupling. The effective one-gluon exchange kernel, as derived from the leading-order effective Hamiltonian approach, produces a divergent results, as is well known in the literature (e.g., Refs. Krautgartner:1991xz ; Gubankova:1997mq ; ManginBrinet:2003nm ). This divergence is the result of the high momentum contribution from the spin non-flip part of the Hamiltonian matrix elements. It can be easily seen from the power counting in transverse momenta. In Ref. Li:2015zda , we adopted a UV counterterm proposed by Krautgärtner, Pauli and Wölz (KPW) Krautgartner:1991xz (cf. Refs. Glazek:1992bs ; Trittmann:1997xz ; Wiecki:2014ola ; Lamm:2016djr ). However, the KPW counterterm is non-covariant, and introduces a major source of violation of the rotational symmetry that is manifested in the spectrum. With asymptotic freedom, the UV divergence associated with the one-gluon exchange kernel is absent. Therefore, the non-covariant KPW counterterm is not needed and we omit it in the present work. As we will see below, the rotational symmetry is improved compared to the results of Ref. Li:2015zda .

III Hamiltonian formalism

III.1 Eigenvalue Equation

The mass spectrum and the wave functions are obtained from diagonalizing the effective light-front Hamiltonian operator (1):

[TABLE]

where $P=(P^{-},P^{+},\vec{P}_{\perp})$ is the 4-momentum of the particle; $j$ and $m_{j}$ are the particle’s total angular momentum and the magnetic projection, respectively.

The Fock space representation of quarkonium reads:

[TABLE]

The coefficients of the expansion, $\psi^{(m_{j})}_{s\bar{s}/h}(\vec{k}_{\perp},x)$ are the valence sector LFWFs with $s$ ($\bar{s}$) representing the spin of the quark (antiquark). The quark and anti-quark creation operators $b^{\dagger}$ and $d^{\dagger}$ satisfy the canonical anti-commutation relations,

[TABLE]

where $\delta^{3}(p-p^{\prime})\equiv\delta(p^{+}-p^{\prime+})\delta^{2}(\vec{p}_{\perp}-\vec{p}^{\prime}_{\perp})$. We have kept only the $q\bar{q}$ sector while, in principle, the $q\bar{q}g$ sector can be included by, e.g., a perturbative treatment Leitner:2010nx . The hadron state vector can be orthonormalized according to the one-particle state [cf. Eq. (6)]:

[TABLE]

Then, the orthonormalization of the LFWFs reads,

[TABLE]

Note that different hadron states with the same quantum numbers, such as $J/\psi$ and $\psi^{\prime}$, are also orthogonal. It is also useful to introduce LFWFs in the transverse coordinate space:

[TABLE]

with orthonormalization,

[TABLE]

Parity $\mathcal{P}$ is a dynamical symmetry on the light front, as it swaps light-front coordinate $x^{-}$ and light-front time $x^{+}$. The mirror parity $m_{P}\equiv\mathcal{R}_{x}(\pi)\mathcal{P}$, which only flips one of the transverse spatial coordinates ($x^{1}$), survives as a kinematical symmetry in light-front dynamics. The eigenvalue equations related to the mirror parity $\hat{m}_{P}$ and the charge conjugation $\hat{C}$ are Krautgartner:1991xz ; Trittmann:1997xz ; Brodsky:2006ez ; Li:2015zda :

[TABLE]

Here $\mathsf{P}$ and $\mathsf{C}$ are the parity and charge conjugation quantum numbers, respectively; and $P=(P^{-},P^{+},P^{1},P^{2})$ is the total 4-momentum of the particle, $\widetilde{P}=(P^{-},P^{+},-P^{1},P^{2})$. $\bar{h}$ represents the antiparticle of hadron $h$.

Particles are further classified by the eigenvalues of the intrinsic angular momenta $\{\vec{\mathcal{J}}^{2},\mathcal{J}_{z}\}$, viz

[TABLE]

On the light front, $\vec{\mathcal{J}}^{2}$ is dynamical and, in principle, it should be diagonalized simultaneously with the light-front Hamiltonian operator $P^{-}$ to obtain the total angular momentum $j$ Brodsky:1997de . Accordingly, in a truncated and regularized model space, $\vec{\mathcal{J}}^{2}$ may not commute with $P^{-}$, and the rotational symmetry is only approximate (see Fig. 4). To extract $j$, we compute the mass eigenvalues from all $m_{j}$ sectors. We count the multiplicity of the nearly-degenerate mass eigenstates with the further help of the mirror parity, charge conjugation and other relevant quantities444For example, the decay constants.. For this scheme to succeed, the degeneracies have to be observed in the results with sufficient accuracy to resolve ambiguities.

III.2 Basis Representation

The eigenvalue equation (4) 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. 2a):

[TABLE]

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. 2b):

[TABLE]

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$ or $\beta$ from now on.

In the presence of the one-gluon exchange term, we use these analytic functions as a basis to expand the LFWFs in,

[TABLE]

Here the coefficients $\psi_{h}(n,m,l,s,s^{\prime})$ are obtained from diagonalization. The basis is constructed to conserve the magnetic projection of the total angular momentum: $m_{j}=m+s+s^{\prime}$.

Performing a 2D Fourier transformation gives the LFWFs in coordinate space. The Fourier transformation of a HO function is a HO function with a relative phase, which simplifies the expression greatly.

[TABLE]

Here $\widetilde{\phi}_{nm}$ is the 2D HO in coordinate space:

[TABLE]

In practical calculations, the basis is truncated and wave functions are obtained in the basis expansion. Following Refs. Vary:2009gt ; Wiecki:2014ola ; Li:2015zda , we truncate the transverse and the longitudinal bases separately by their energies:

[TABLE]

As such, the $N_{\max}$-truncation provides a natural pair of UV and IR cutoffs: $\Lambda_{\textsc{uv}}\simeq b\sqrt{N_{\max}}$, $\lambda_{\textsc{ir}}\simeq b/\sqrt{N_{\max}}$, where $b=\kappa$ is the oscillator basis energy scale parameter. $L_{\max}$ represents the resolution of the basis in the longitudinal direction. Namely, the basis cannot resolve physics at: $\Delta x\lesssim L_{\max}^{-1}$ Wiecki:2014ola . The complete basis is reached by taking $N_{\max}\to\infty,L_{\max}\to\infty$.

The eigenvalues of the parity and charge conjugation operators can be extracted from the basis representation of the LFWFs as Li:2015zda ,

[TABLE]

III.3 Generalizing Light-Front Holography

Before proceeding to the full diagonalization, it is worth looking at the results without the one-gluon exchange, where the solutions are analytical. The mass eigenvalues are:

[TABLE]

Here $l$ is the longitudinal quantum number, not the orbital angular momentum. The corresponding wave functions are:

[TABLE]

States may be identified according to their mass spectrum with the help of parity $\mathsf{P}$ and charge conjugation $\mathsf{C}$. The quarkonium ground state (1S) is identified with $n=0,m=0,l=0$:

[TABLE]

In the literature, a commonly-used way to incorporate quark masses in the AdS/QCD wave function is through the invariant mass ansatz (IMA) Brodsky:2008pg , viz,

[TABLE]

Figure 3 compares the purely longitudinal part of our ground-state wave function with that of the IMA wave function. Our longitudinal wave function becomes almost identical to the IMA wave function in both the chiral limit and the heavy quark limit, except near the endpoints. This reflects the fact that rotational symmetry is restored in the non-relativistic limit with our choice of longitudinal basis functions.

The first excited state (1P) is identified with $n=0,m=\pm 1,l=0$ or $n=0,m=0,l=1$, noting that for heavy quarkonium, the term $\kappa^{4}/(m_{q}+m_{\bar{q}})^{2}l(l+1)$ is small comparing to the remaining terms. There are four 1P states: $\chi_{0}$ ($0^{++}$), $\chi_{1}$ ($1^{++}$), $\chi_{2}$ ($2^{++}$) and $h$ ($1^{+-}$). Let us focus on $h$ and restrict the discussion to $m_{j}=0$. From Eq. (20), we conclude: $-1=\mathsf{C}=(-1)^{m+l+1}(-1)^{s+1}$, where $s$ is the total spin, viz $s=0$ for singlet and $s=1$ for triplet. Apparently, for both sets of quantum numbers ($m=\pm 1,l=0$ or $m=0,l=1$), $s=0$. From Eq. (19), $-1=(-1)^{j}\mathsf{P}=(-1)^{m}(-1)^{s+1}$, implying $m=0$. Therefore, the correct quantum numbers for $h$ meson ($1^{+-}$) are $n=0,m=0,l=1$ with a singlet spin configuration, which is consistent with the non-relativistic quantum number assignment $1{\,}^{1}\!P_{0}$. Note that the orbital motion is excited through the longitudinal direction but not the transverse direction. This cannot be obtained from IMA555In the literature, the longitudinal excitations are typically obtained from modeling the spin structure via the spinor wave function $\bar{u}\Gamma v$. However, the longitudinal profile of the spinor wave function is qualitatively different from the holographic wave function. .

IV Numerical Results

We apply the model to heavy quarkonia (charmonium and bottomonium), where the quark masses are large and the radiative corrections are negligible. Therefore these are ideal systems to test our model. The model parameters are summarized in Table 1.

As mentioned, we fixed $\alpha_{s}(0)=0.6$. For fixed $N_{\max}$ and $L_{\max}$, we use experimental data to fit the confining strength $\kappa$ and the effective quark mass $m_{q}$ ($m_{c}$ and $m_{b}$) using the mass eigenvalues in the $m_{j}=0$ sector. We employ the experimental values, compiled by the Particle Data Group (PDG) Agashe:2014kda , below the open charm or open bottom threshold. We also introduced a small mass parameter $\mu_{g}=0.02\,\mathrm{GeV}$ to regularize the integrable Coulomb singularity in the energy denominator and to avoid numerical instability666Our numerical method is designed such that no singularity is encountered in the actual calculation. Nevertheless, we introduced this parameter, smaller than all other energy scales, to further tame the integrable singularity. . As has been shown in previous work with fixed coupling, the mass eigenvalues are converged with respect to $\mu_{g}\to 0$ within the numerical precision.

The effective Hamiltonian (1) is diagonalized for various $m_{j}$ sectors. Fig. 4 shows a representative spectrum as a function of $m_{j}$. The spectrum is symmetric with respect to $\pm m_{j}$, a consequence of the mirror parity symmetry (11). The discrete quantum numbers $m_{P}=(-i)^{2j}\mathsf{P}$ and $\mathsf{C}$ are computed to help identify states as mentioned. Total spin $\langle\vec{s}^{2}\rangle=s(s+1)$ as an approximate quantum number is also exploited. States with the same $j$ but different $m_{j}$’s are not exactly degenerate owing to the violation of the rotational symmetry. As is seen in Fig. 4, the approximate degeneracies are easily visible, at least for low-lying states. So the multiplicities, together with $m_{P}$, $\mathsf{C}$, $s$ and the constraints:

[TABLE]

can be employed to deduce the full set of quantum numbers $n\,{}^{2s+1}\!\ell_{j}$ or $j^{\mathsf{P}\mathsf{C}}$, where $\ell$ is the total orbital angular momentum, $n$ the radial quantum number. We also cross-check the state identification with the decay constants and the wave functions themselves (see Sect. V).

IV.1 Spectroscopy

The reconstructed spectra at $N_{\max}=L_{\max}=32$ are presented in Fig. 5. In these figures, we use boxes to indicate the spreads of the mass eigenvalues from different $m_{j}$. The mean values, marked by dashed bars, are defined as:

[TABLE]

where $M_{m_{j}}$ is the mass eigenvalue associated with the magnetic projection $m_{j}$. This definition is motivated by the covariant light-front analysis of relativistic bound states in Refs. Carbonell:1998rj ; ManginBrinet:2003nm . On the other hand, the mass spreads $\delta_{j}M\equiv\max M_{m_{j}}-\min M_{m_{j}}$ measure the violation of the rotational symmetry. We also introduce the mean spread:

[TABLE]

For charmonium (bottomonium) states evaluated by PDG below the threshold, the mean mass spread is 17 MeV (8 MeV), improving our previous results Li:2015zda by a factor of $\sim$3 ($\sim$2). More comparison between the results of this work and those of Ref. Li:2015zda is collected in Table 2.

Our light-front Hamiltonian approach yields states with high angular and radial excitations, which are not easily accessible in some other methods. No exotic quantum numbers emerge from our calculation, as is expected from the two-body truncation. In bottomonium, predictions are made for various states below the $B\overline{B}$ threshold, as also predicted in other approaches (e.g., Crater:2010fc ; McNeile:2012qf ; Hilger:2014nma ). The quality of the spectra can be measured by the root mean squared (r.m.s.) deviation from the experimentally measured values. For charmonium (bottomonium), the r.m.s. mass deviation is 31 MeV (38 MeV), improving the fixed coupling results Li:2015zda by as much as $\sim$40% ($\sim$20%). See Table 2 for further comparisons. Our spectroscopy is competitive with those obtained from other relativistic models Crater:2010fc ; McNeile:2012qf ; Hilger:2014nma ; Leitao:2016bqq . Not only are the mass spectra improved, the spread of the mass eigenvalues $\overline{\delta_{j}M}$ due to the violation of rotational symmetry, is also significantly reduced as mentioned. A related issue is the quenching of the hyperfine splitting found within the fixed coupling results. With the running coupling, this issue is resolved and the hyperfine splittings are consistent with the experimental values, as shown in Fig. 5 (cf. Fig. 6). Therefore, the violation of the rotational symmetry is significantly reduced.

Figure 6 shows the trends of the charmonium mass eigenvalues as functions of $N_{\max}^{-1}$ (with $N_{\max}=L_{\max}$). The left panel presents the convergence trends of the ground-state masses ($\eta_{c}$ and $J/\psi$). The right panel presents the convergence trends for the hyperfine splittings between 1S ($M_{J/\psi}-M_{\eta_{c}}$) and 2S ($M_{\psi^{\prime}}-M_{\eta^{\prime}_{c}}$) states. Two sets of parameters are used: the fix-parameter results use model parameters from the $N_{\max}=L_{\max}=32$ fit; the refit-parameter calculation refits the model parameters for each $N_{\max}=L_{\max}$. Smooth extrapolations are made using three types of functions: $a+b/N_{\max}+c/N_{\max}^{2}$ (solid), $a+b\exp(-cN_{\max})$ (dashed), $a+b\exp(-c\sqrt{N_{\max}})$ (dot-dashed). Both hyperfine splittings, 1S and 2S, show reasonable convergence in the complete basis limit ($N_{\max}^{-1}=L_{\max}^{-1}=0$).

Table 3 compares the spectroscopy obtained from different $N_{\max}=L_{\max}$ fits. While results from different $N_{\max}=L_{\max}$ are well converged, the r.m.s. deviation decreases as $N_{\max}=L_{\max}$ increases. In the present work, we adopt $N_{\max}=L_{\max}=32$ for our presented results, unless otherwise specified.

IV.2 Decay Constants

Decay constants are defined as the local vacuum-to-hadron matrix elements:

[TABLE]

Here only the “good” currents (the “+” component) are used. The corresponding LFWF representation reads Lepage:1980fj ,

[TABLE]

For this calculation, we choose $N_{\max}=8$ for charmonium and $N_{\max}=32$ for bottomonium, roughly corresponding to $\Lambda_{\textsc{uv}}\triangleq\kappa\sqrt{N_{\max}}\approx 1.7m_{q}$, where $\Lambda_{\textsc{uv}}$ is the UV regulator, and $m_{q}$ is the heavy quark mass. This choice is motived by the competition between the needs for both a better basis resolution and a lower UV scale since our model does not incorporate radiative corrections. We also provide an indicator for sensitivity by altering the basis truncation parameter $N_{\max}$. The resulting charmonium and bottomonium decay constants are shown in Fig. 7, which also collects PDG values converted from dilepton or diphoton decay widths Agashe:2014kda , Lattice Davies:2010ip ; McNeile:2012qf ; Donald:2012ga ; Colquhoun:2014ica and Dyson-Schwinger/Bethe-Salpeter equations (DSE/BSE or DSE, Blank:2011ha ; see also Ding:2015rkn ) results for comparison. Our results fall into the ballpark of the PDG values as well as those from other approaches wherever available.

IV.3 Radii

Classically and in non-relativistic quantum mechanics, the root-mean-square charge (mass) radius is the expectation value of the displacement operator that characterizes the charge (mass) distribution of the system. In quantum field theory, no such local position operator is allowed and, instead, the form factors are defined as the slope of the charge (gravitational) form factor at zero momentum transfer:

[TABLE]

Remarkably, in LFWF representation Brodsky:2000ii , this definition exactly restores the charge (mass) distribution interpretation Li:2016wwu . For example, for (pseudo-)scalar mesons in the two-body approximation,

[TABLE]

Here $\widetilde{\psi}$ are LFWFs in transverse coordinate space. $\vec{\zeta}_{\perp}\triangleq\sqrt{x(1-x)}\vec{r}_{\perp}$ is Brodsky and de Téramond’s holographic variable Brodsky:2014yha , $\vec{b}_{\perp}\triangleq(1-x)\vec{r}_{\perp}$ is Burkardt’s impact parameter Burkardt:2000za . This relation is also valid when higher Fock sector contributions are included if the we define $\vec{\zeta}_{\perp}$ and $\vec{b}_{\perp}$ in the $n$-body Fock sector as,

[TABLE]

where $\vec{R}_{\perp}\triangleq\sum_{i}x_{i}\vec{r}_{i\perp}$ is the transverse center of the system, $e_{i}$ is the charge number of the $i$-th constituent, and $\sum_{i}e_{i}\equiv Q$.

Due to charge conjugation symmetry, the charge radii of quarkonium vanishes. Here we define a fictitious charge radii by considering only the charge of the quark. With this definition, the “charge” radii are the same as the mass radii in the non-relativistic limit, which suggests that their difference is a pure relativistic effect. Fig. 8 presents the r.m.s. charge and mass radii of scalar and pseudo-scalar mesons. In our results, the mass radii are in general smaller than the charge radii and the difference is reduced in the heavier system (bottomonium). Fixed $\alpha_{s}$ BLFQ results (BLFQ15, Li:2015zda ) as well as earlier results from quenched Lattice calculation Dudek:2006ej and DSE Maris:2006ea are included in Fig. 8 for comparison. Our results are systematically smaller. From the trend with respect to basis truncation $N_{\max}=L_{\max}$, UV physics and/or higher Fock sector contributions may be expected to produce significant corrections to our results for radii.

V Wave Functions, Amplitudes and Distributions

V.1 Light-Front Wave Functions

Wave functions offer first-hand insight into the system. They play a central role in evaluating hadronic observables and light-cone distributions, and are an indispensable tool for investigating exclusive processes in deep inelastic scattering Chen:2016dlk . Compared with the widely used phenomenological LFWFs in the literature, our wave functions generalize the AdS/QCD wave functions and provide unified access to ground and excited states. In particular, the spin structure is generated from the one-gluon exchange and its interplay with the confining interaction.

In this section, we present the valence sector wave functions. Heavy quarkonium is an ideal system to explore the qualitative features of the wave functions, as they can be compared with the familiar non-relativistic quantum mechanical wave functions. We show LFWFs with different polarizations and spin alignments: $\psi_{s\bar{s}}^{\lambda}(\vec{k}_{\perp},x)$. For each spin configuration, the orbital angular momentum projection $m_{\ell}=\lambda-s_{1}-s_{2}$ is definite ($\lambda\equiv m_{j}$). Hence, the angular dependence of the wave function factorizes: $\psi_{s\bar{s}}^{\lambda}(\vec{k}_{\perp},x)=\Psi_{s\bar{s}}^{\lambda}(k_{\perp},x)\exp({\mathrm{i}}m_{\ell}\theta)$, with $\theta\equiv\arg\vec{k}_{\perp},k_{\perp}\equiv|\vec{k}_{\perp}|$. To visualize the wave functions, we drop the phase $\exp({\mathrm{i}}m_{\ell}\theta)$, while retaining the relative sign $\exp({\mathrm{i}}m_{\ell}\pi)=(-1)^{m_{\ell}}$ for negative values of $k_{\perp}$. Namely, we plot:

[TABLE]

We also define: \psi_{\uparrow\downarrow\pm\downarrow\uparrow}^{\lambda}(\vec{k}_{\perp},x)\equiv\frac{1}{\sqrt{2}}\big{[}\psi_{\uparrow\downarrow}^{\lambda}(\vec{k}_{\perp},x)\pm\psi_{\downarrow\uparrow}^{\lambda}(\vec{k}_{\perp},x)\big{]}. The full set of results is collected in supplemental materials. Here we focus on some selected results.

Figure 9 shows the LFWFs of the charmed ground-state pseudo scalar $\eta_{c}(1S)$. There are two independent components: $\psi_{\uparrow\downarrow-\downarrow\uparrow}(\vec{k}_{\perp},x)$ and $\psi_{\downarrow\downarrow}(\vec{k}_{\perp},x)=\psi^{*}_{\uparrow\uparrow}(\vec{k}_{\perp},x)$. The number of independent components is not a priori the same in different relativistic approaches. One of the components is related to the non-relativistic wave functions, whereas the other one is of purely relativistic origin and becomes negligible in the non-relativistic limit. In covariant light-front dynamics, the extra component depends on the orientation of the quantization surface Carbonell:1998rj ; Leitner:2010nx . Its existence ensures the rotational symmetry, albeit not exactly in our model Carbonell:1998rj . The Lorentz structure of the pseudo scalar wave function can be written as Carbonell:1998rj ; Leitner:2010nx ,

[TABLE]

where $\gamma^{+}=\gamma^{0}+\gamma^{3}$. Let $n=(1,0,0,-1)$ be a null vector perpendicular to the quantization surface. $\gamma^{+}=n_{\mu}\gamma^{\mu}$, $P^{+}\equiv n_{\mu}P^{\mu}$, both depending on the orientation of the quantization surface.

For charmonium, the dominate component is the singlet $\psi_{\uparrow\downarrow-\downarrow\uparrow}$ and its wave function resembles an S-wave. In the non-relativistic limit, the longitudinal momentum fraction $x$ is reduced to: $x\to 1/2+{k_{z}}/(2m_{q})$. Hence, the $x$–$k_{\perp}$ plots in Fig. 9 (central panels) are reduced to the $k_{z}$–$k_{\perp}$ density plots of the non-relativistic wave function, i.e. a slice of the full 3D wave function, in the non-relativistic limit. To visualize the full 3D wave function, one may rotate the density plot along the vertical axis at $k_{\perp}=0$, applying a phase factor $\exp({\mathrm{i}}m_{\ell}\theta)$ as necessary777This is where the relative sign at negative $k_{\perp}$ is useful.. To facilitate the visualization in 3D, we also plot the real part of the wave functions in the transverse plane at $x=0.5$ in Fig. 9 (right panels).

Figure 10 shows the spin singlet components of the charmed and beautified pseudo scalars $\eta_{c}(nS)$ and $\eta_{b}(nS)$. Each of them is the dominant component in their respective systems. The 2S and 3S states show both longitudinal and transverse nodes, consistent with the non-relativistic wave functions. Therefore, the non-relativistic picture emerges in heavy quarkonium as expected. Note that the node structure spans a broad kinematical region [$x\sim(0.2\text{--}0.8)$] in charmonium, extending beyond the naïve non-relativistic scope: $|x-\frac{1}{2}|\ll 1$.

It is also interesting to compare the charmonium and bottomonium LFWFs, as shown in Fig. 11. Bottomonium is associated with a larger mass scale and is broader in the transverse momentum direction. On the other hand, bottomonium is more non-relativistic compared with charmonium, and hence in the longitudinal direction its wave functions are narrower. Recall that in the non-relativistic limit, the quarkonium distribution amplitude is a Dirac delta: $\phi(x)\propto\delta(x-\frac{1}{2})$.

Figure 12 compares selected spin configurations of the charmed vector mesons: $J/\psi$ with its “angular” excitation $\psi(1D)$. The dominant components of $J/\psi$ are $\psi^{\lambda=0}_{\uparrow\downarrow+\downarrow\uparrow}$ (Fig. 12a) and $\psi^{\lambda=1}_{\uparrow\uparrow}$ (see supplemental materials) — both are S-wave. The D-wave components (e.g. Figs. 12b & 12c) are small but non-vanishing in $J/\psi$ as a result of S-D mixing. Similar sub-dominant components due to relativity are often missing in phenomenological vector meson wave functions888Very often, the spin structure of the phenomenological vector meson wave function is borrowed from the photon wave function, which is obtained via light-cone perturbation theory., e.g., boosted Gaussian wave function Chen:2016dlk . The dominant components of $\psi(1D)$ are $\psi^{\lambda=0}_{\uparrow\downarrow+\downarrow\uparrow}$ (Fig. 12d), $\psi^{\lambda=0}_{\downarrow\downarrow}$ (Fig. 12e), and $\psi^{\lambda=1}_{\downarrow\downarrow}$ (Fig. 12f). It is evident that they resemble the non-relativistic D-waves $Y_{20}(\hat{k})$, $Y_{21}(\hat{k})$ and $Y_{22}(\hat{k})$, where $Y_{\ell m}(\hat{k})$ are the spherical harmonics. This becomes more evident when LFWFs in the transverse plane ($k_{x}$–$k_{y}$) are considered (see Fig. 13). Fig. 14 displays $\Upsilon(2D)$, a state consisting of both radial and angular excitations.

V.2 Distribution Amplitudes

LFWFs provide unique access to light cone distributions by integrating out the transverse momentum Brodsky:1997de . Among those, the distribution amplitudes (DAs) and the parton distribution functions (PDFs) control the exclusive and inclusive processes at large momentum transfer, respectively Lepage:1980fj .

DAs are defined from the light-like separated gauge invariant vacuum-to-meson matrix elements Lepage:1980fj ; Bodwin:2006dm . In light-front formalism, the leading-twist DAs within the light-cone gauge for pseudo-scalar and vector mesons999In the present work, we focus on the longitudinal DA for vector mesons. are Bodwin:2006dm ; Braguta:2006wr ; Braguta:2007fh :

[TABLE]

where $f_{P,V}$ are the decay constants (see Sect. IV.2). $M_{P,V}$ are the mass eigenvalues. $e^{\mu}_{\lambda}(p)$ is the polarization vector. The non-local matrix elements as well as the DAs depend on the scale $\mu$, the renormalization scale or UV cutoff. In these definitions, DAs are normalized to unity, viz:

[TABLE]

In LFWF representation, DAs can be written as Lepage:1980fj ,

[TABLE]

Here $\psi_{\uparrow\downarrow\pm\downarrow\uparrow}=(\psi_{\uparrow\downarrow}\pm\psi_{\downarrow\uparrow})/\sqrt{2}$ as defined above and the minus (plus) sign is associated with the pseudo-scalar (vector) state. The UV cutoff is taken as $k_{\perp}/\sqrt{x(1-x)}\lesssim\mu$ (see, e.g., Refs. Lepage:1980fj ; Krautgartner:1991xz ; Zhang:1993dd ). In the basis representation, the truncation parameter $N_{\max}$ provides a natural UV regulator $\mu\approx\kappa\sqrt{N_{\max}}$ and no hard cutoff is needed in the integration.

Figure 15 compares the ground-state vector meson ($J/\psi$ and $\Upsilon$) DAs with predictions from BLFQ and AdS/QCD with or without IMA Brodsky:2014yha ; Vega:2009zb ; Branz:2010ub ; Swarnkar:2015osa . Calculations using pure basis functions are also presented (AdS/QCD + LC), which turns out to be very close to AdS/QCD + IMA (cf. Fig. 3), but very different from the full diagonalization (“BLFQ”) results. In fact, the BLFQ results move towards the pQCD asymptotics as the scale increases. Obviously, the one-gluon exchange interaction plays an important role at short distance as is expected. DAs of S-wave heavy quarkonia are shown in Fig. 16. The difference between the pseudo-scalar mesons and the accompanying vector mesons are, again, driven by the one-gluon exchange interaction. The shape of the excited state DAs is consistent with what has been obtained from other methods, e.g. QCD sum rule Braguta:2007tq , wherever available. The basis functions are optimized for long-distance physics, i.e., confinement, and DAs are sensitive to short-distance physics. The mismatch as a finite-basis effect is clearly visible around the endpoints in these figures.

It is useful to compute the moments in order to quantitatively compare with other approaches. The $n$-th moment is defined as,

[TABLE]

Table 4 compares the first few moments of selected heavy quarkonia states obtained from various approaches. Results from other approaches, including non-relativistic QCD (NRQCD, Bodwin:2006dn ), QCD sum rule (QCDSR, Braguta:2006wr ; Braguta:2007fh ; Braguta:2007tq ), light-front quark model (LFQM, Choi:2007ze ) and Dyson-Schwinger/Bethe-Salpeter equations (DSE, Ding:2015rkn ), are shown for comparison. In all these approaches, moments are computed at the effective heavy quark mass scale $\mu\simeq m_{q}$, with the exception of DSE at $\mu=2\,\mathrm{GeV}$. We provide results at $\mu\approx 1.7m_{q}$, corresponding to $N_{\max}=L_{\max}=8$ for charmonium and $N_{\max}=L_{\max}=32$ for bottomonium. For the sake of convenience, we also provide moments at the effective heavy quark mass scale $\mu=m_{q}$ (“BLFQ*”) through simple extrapolation (for charmonium) or interpolation (for bottomonium). The $3\sigma$ (99.75% C.L.) extrapolation or interpolation errors (prediction intervals) are included. Our results are in reasonable agreement with various other approaches, though relativistic models, including ours, are systematically larger than those of NRQCD. Results from pQCD asymptotics and AdS/QCD of Brodsky and de Téramond (AdS/QCD, Brodsky:2014yha ) are not particularly applicable for heavy quarkonia at the heavy quark mass scale and are simply included for completeness. The second moment can be used to estimate the relative velocity of the partons: $\langle v^{2}\rangle\approx 3\langle\xi^{2}\rangle$, viz

[TABLE]

V.3 Parton Distributions

The quark Parton Distribution Function (PDF) $q(x;\mu)$ is the probability of finding a collinear quark carrying momentum fraction $x$ up to scale $\mu$. In the light-front formalism, it can be obtained by integrating out the transverse momentum of the squared wave function:

[TABLE]

Within the two-body approximation, the PDF and its first moment are normalized to unity [cf. Eq. (8)]:

[TABLE]

Figure 17 shows PDFs of (pseudo-)scalar quarkonia. They exhibit distinctive features compared with DAs. In particular, there is no dip at $x=1/2$ in excited-state PDFs, in contrast to DAs. There appear to be ripples on the downward slopes of PDFs for excited states as may be expected from contributions of longitudinally excited basis functions.

The generalization of PDFs, known as generalized parton distributions (GPDs), unifying PDFs and form factors, provide more insights into the system, and are directly related to experiments Burkardt:2000za ; Diehl:2003ny . Wigner distributions are more general quantities unifying GPDs and the transverse momentum distributions. In principle, all of them are accessible through LFWFs, at least in some kinematical regime (e.g. Swarnkar:2015osa ). For example, in the zero skewedness limit, the impact parameter GPD $q(x,\vec{b}_{\perp})$ of Burkardt Burkardt:2000za is related to the LFWFs simply by,

[TABLE]

VI Summary and Discussions

We present a light-front model for quarkonium that incorporates light-front holographic QCD and the one-gluon exchange interaction with a running coupling. We solve the model in the Hamiltonian approach with a basis function expansion. We obtain mass spectroscopy and the light-front wave functions. The spectroscopy agrees with the PDG data within 30–40 MeV of r.m.s. mass deviation for states below the open flavor threshold. The overall quality improves the previous work that employed a fixed strong coupling and a non-covariant counterterm. The wave functions reveal rich structures, especially for excited states. Through analysis and comparison, we find these structures are consistent with the standard non-relativistic picture. From these wave functions, we also compute the decay constants, r.m.s. radii, distribution amplitudes, and parton distributions. Our results appear to be in reasonable agreement with those from other approaches wherever available.

This work is an attempt to improve light-front holographic QCD approach by adding realistic QCD interactions. In particular, we show that while rotational symmetry is broken due to truncation, the extraction of angular momentum $j$ is feasible and reliable Trittmann:1997xz ; Brodsky:2006ez . The Hamiltonian formalism and the basis function approach enable us to access a wide range of states, including radial and angular excited states extending over all known excited states and beyond. The obtained light-front wave functions allow us to directly compute hadronic distributions such as distribution amplitudes as well as hadronic observables. It should be emphasized that these attractive features are not limited to the present effective model—they are the shared advantages within the light-front Hamiltonian formalism Bakker:2013cea .

We did not include self-energy in solving the heavy quarkonia. However, radiative corrections may become important in evaluating some observables as we employ more realistic field-theory dynamics. The calculation of the decay constants illustrates this particular challenge. As we move to the light sector, the consistent inclusion of self-energies and renormalization issues may become more acute if one wants to address additional phenomena within QCD such as chiral symmetry breaking. Nevertheless, we believe the present work may serve as a substantial step for developing an elaborate light-front model for hadrons as relativistic bound states.

While the advantages of the basis function expansion is obvious, it nevertheless requires more investigation. The IR and UV scales are tied to the basis truncation parameter $N_{\max}$ and $L_{\max}$. Compared to the wave-equation approach, the UV asymptotics is not easy to analyze. We typically rely on extrapolation of the basis parameters as developed in ab initio nuclear structure calculations Coon:2012ab . In BLFQ, the basis extrapolation requires further study. One investigation was conducted in the context of strong coupling light-front QED and the authors found robust basis extrapolations that are consistent with the wave-equation approach Wiecki:2014ola . The coupling ($\alpha=0.3$), the transverse basis as well as the one-photon exchange kernel used in Ref. Wiecki:2014ola are very similar to the present model.

Future developments should focus on the inclusion of higher Fock sectors and the non-perturbative renormalization (see Ref. Hiller:2016itl for a recent review). In the top-down approach, a systematic non-perturbative renormalization scheme should be developed and non-pertubative dynamics has to be addressed using efficient numerical methods. Notable examples include the full basis light-front quantization (BLFQ, Vary:2009gt ), the renormalization group procedure for effective particles (RGPEP, Glazek:2012qj ; Gomez-Rocha:2015esa ), the Fock sector dependent renormalization (FSDR, Karmanov:2008br ; Li:2015iaw ), and the light-front coupled cluster method (LFCC, Chabysheva:2011ed ). In the bottom-up approach, one is motivated to design appropriate kernels that incorporate important physics while preserving the symmetries. Notable physics goals for hadrons include the radiative corrections, asymptotic freedom and the dynamical chiral symmetry breaking. Incorporating the running coupling is the first step. In both approaches, the current model may serve as a first approximation. See also Refs. Ji:2012ux ; Chang:2013pq ; Ji:2013dva ; Karmanov:2005nv ; Leitao:2017 for some recent works bridging other approaches with the light-front approach.

The applicability of the current model is not restricted to heavy quarkonium. Extensions to other meson and baryon systems, in principle, are straightforward, although new issues have to be addressed in each of these systems.

Acknowledgements

We wish to thank X. Zhao, G. Chen, M.-j. Li, S. Leitão, S.J. Brodsky, G. de Téramond, E. Swanson, J.R. Spence, J.R. Hiller, S.S. Chabysheva, S.D. Głazek, A. Trawiński and M. Gomez-Rocha for valuable discussions. One of us (Y.L.) also wants to thank the hospitality of the High Energy Nuclear Theory Group at the Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou, China, where part of the work is being completed.

This work was supported in part by the Department of Energy under Grant Nos. DE-FG02-87ER40371 and DESC0008485 (SciDAC-3/NUCLEI). Computational resources were provided by the National Energy Research Supercomputer Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

Appendix A Light-Front Coordinates

We adopt natural units throughout the article: $\hbar=c=1$. We roughly follow the convention of Ref. Wiecki:2014ola . The light-front coordinates are defined as $x=(x^{-},x^{+},x^{1},x^{2})\equiv(x^{-},x^{+},\vec{x}_{\perp})$, where $x^{\pm}=x^{0}\pm x^{3}$. The inner product of two 4-vectors is defined as: $a\cdot b=\frac{1}{2}a^{-}b^{+}+\frac{1}{2}a^{+}b^{-}-\vec{a}_{\perp}\cdot\vec{b}_{\perp}$. It should be noted that the determinant of the metric tensor is $\det g=-(1/4)$.

The Lorentz invariant phase space measure is

[TABLE]

where $\vartheta(z)$ is the unit step function. The one-particle state is normalized as: $\langle p,j,m_{j}|p^{\prime},j^{\prime},m_{j}^{\prime}\rangle=2p^{+}\vartheta(p^{+})(2\pi)^{3}\delta^{3}(p-p^{\prime})\delta_{jj^{\prime}}\delta_{m_{j},m_{j}^{\prime}}$, where the Dirac delta is defined as $\delta^{3}(p)\equiv\delta(p^{+})\delta^{2}(\vec{p}_{\perp})$.

Appendix B Few-Body Kinematics

We define boost-invariant momenta from the single-particle momenta $\{p^{+}_{i},\vec{p}_{i\perp}\}$ as,

[TABLE]

$x_{i}$ are the longitudinal light-front momentum fractions; and $k_{i\perp}$ are the transverse relative momenta. They satisfy:

[TABLE]

The $n$-body phase space integration measure factorizes:

[TABLE]

The invariant mass squared of the $n$-body Fock state is:

[TABLE]

Appendix C Spinors

The $u$, $v$ spinors are defined as,

[TABLE]

where \chi_{+}=(1,0,0,0)^{\mathchoice{\raisebox{0.0pt}{\displaystyle\intercal}}{\raisebox{0.0pt}{\textstyle\intercal}}{\raisebox{0.0pt}{\scriptstyle\intercal}}{\raisebox{0.0pt}{\scriptscriptstyle\intercal}}},\chi_{-}=(0,1,0,0)^{\mathchoice{\raisebox{0.0pt}{\displaystyle\intercal}}{\raisebox{0.0pt}{\textstyle\intercal}}{\raisebox{0.0pt}{\scriptstyle\intercal}}{\raisebox{0.0pt}{\scriptscriptstyle\intercal}}}; $\gamma^{\pm}=\gamma^{0}\pm\gamma^{3}$; $s=\pm$ is the light-front helicity. The $u$, $v$ spinors defined above are polarized in the $z$-direction (or longitudinal direction):

[TABLE]

$S_{z}\equiv\frac{{\mathrm{i}}}{2}\gamma^{1}\gamma^{2}$ and follow the standard orthonormality

[TABLE]

and completeness

[TABLE]

Here are some useful identities:

[TABLE]

The spinor matrix elements for the one-gluon exchange are collected in Table 5 (see also Table I of Ref. Wiecki:2014ola ).

Appendix D Polarization Vectors

gauge bosons

The polarization vector of a gauge boson in light-cone gauge $A^{+}=0$ is:

[TABLE]

where $\vec{\epsilon}_{\pm\perp}=\frac{1}{\sqrt{2}}(-1,\mp{\mathrm{i}})$. The polarization vector defined here satisfies:

• •

$k_{\mu}\varepsilon^{\mu}_{\lambda}(k)=0$;

• •

$\varepsilon^{\mu}_{\lambda}(k)\varepsilon^{*}_{\lambda^{\prime}\mu}(k)=-\delta_{\lambda,\lambda^{\prime}}$;

• •

helicity sum:

[TABLE]

Here $n=(1,0,0,-1)$ is a light-like 4-vector ($n_{\mu}n^{\mu}=0$) perpendicular to the light front.

vector bosons

The polarization vector for the a vector boson:

[TABLE]

where $m$ is the mass of the vector boson. The polarization vector defined here satisfies:

• •

$k_{\mu}e_{\lambda}^{\mu}(k)=0$;

• •

$e^{\mu}_{\lambda}(k)e^{*}_{\lambda^{\prime}\mu}(k)=-\delta_{\lambda,\lambda^{\prime}}$;

• •

spin sum:

[TABLE]