Heavy flavor dissociation in the framework of multi-body Dirac equations

TL;DR

This paper investigates heavy flavor particles at finite temperature using a relativistic potential model and an improved three-body Dirac equation solution, revealing sequential dissociation patterns in quark-gluon plasma.

Contribution

It introduces an improved method for solving the three-body Dirac equation and applies lattice QCD-based potentials to study heavy flavor dissociation at finite temperature.

Findings

Mesons survive at higher temperatures than baryons.

Heavy flavor binding energies decrease with increasing temperature.

Sequential dissociation of mesons and baryons in quark-gluon plasma.

Abstract

We study heavy flavor properties at finite temperature in the framework of a relativistic potential model. With an improved method to solve the three-body Dirac equation, we determine a universal set of model parameters for both mesons and baryons by fitting heavy flavor masses in vacuum. Taking heavy quark potential from lattice QCD simulations in hot medium, we systematically calculate heavy flavor binding energies and averaged sizes as functions of temperature. The meson and baryons are separately sequentially dissociated in the quark-gluon plasma, and the mesons can survive at higher temperature due to the stronger potential between quark-antiquark pairs than that between quark-quark pairs.

Tables

Table 1. Table 1: The universal set of parameters of the potential model.

$m_u=m_d=0.135$ GeV
$m_s=0.263$ GeV
$m_c=1.400$ GeV
$m_b=4.773$ GeV
${\alpha }_{qq}={\alpha }_{q\overline{q}}/2.22=0.20$
${\sigma }_{qq}={\sigma }_{q\overline{q}}/2.04=0.09$ GeV2

Table 2. Table 2: The experimentally measured  Tanabashi:2018oca and model calculated heavy flavor meson masses M E subscript 𝑀 𝐸 M_{E} and M T subscript 𝑀 𝑇 M_{T} . D R = | ( M T − M E ) / M E | subscript 𝐷 𝑅 subscript 𝑀 𝑇 subscript 𝑀 𝐸 subscript 𝑀 𝐸 D_{R}=|(M_{T}-M_{E})/M_{E}| is the relative difference.

Meson	$J^P$	$M_m^E$	$M_m^T$	$D_R$
		(GeV)	(GeV)	(%)
$D^0$	$0^{-}$	1.865	1.940	$4.0$
$D^{\ast 0}$	$1^{-}$	2.007	2.066	$3.0$
$D^{+}$	$0^{-}$	1.870	1.940	$3.8$
$D^{\ast +}$	$1^{-}$	2.010	2.066	$2.8$
$D_s$	$0^{-}$	1.968	2.028	$3.1$
$D_s^{\ast }$	$1^{-}$	2.112	2.157	$2.1$
${\eta }_c$	$0^{-}$	2.984	2.990	$0.2$
${\eta }_c\left(2S\right)$	$0^{-}$	3.637	3.609	$0.8$
$h_{c1}$	$1^{+}$	3.525	3.506	$0.5$
$J/\psi$	$1^{-}$	3.097	3.123	$0.8$
$\psi \left(2S\right)$	$1^{-}$	3.686	3.701	$0.4$
${\chi }_{c0}$	$0^{+}$	3.415	3.442	$0.8$
${\chi }_{c1}$	$1^{+}$	3.511	3.504	$0.2$
${\chi }_{c2}$	$2^{+}$	3.556	3.519	$1.0$
$B^{-}$	$0^{-}$	5.279	5.326	$0.5$
$B^{-\ast }$	$1^{-}$	5.325	5.371	$0.9$
$B^0$	$0^{-}$	5.280	5.326	$0.9$
$B^{0\ast }$	$1^{-}$	5.325	5.371	$0.9$
$B_s$	$0^{-}$	5.367	5.408	$0.8$
$B_s^{\ast }$	$1^{-}$	5.415	5.458	$0.8$
${\eta }_b$	$0^{-}$	9.399	9.378	$0.2$
${\eta }_b\left(2S\right)$	$0^{-}$	9.999	9.964	$0.3$
$h_{b1}$	$1^{+}$	9.899	9.918	$0.2$
${\rm Y}\left(1S\right)$	$1^{-}$	9.460	9.507	$0.5$
${\rm Y}\left(2S\right)$	$1^{-}$	10.023	10.025	$0.0$
${\chi }_{b0}$	$0^{+}$	9.859	9.878	$0.2$
${\chi }_{b1}$	$1^{+}$	9.893	9.912	$0.2$
${\chi }_{b2}$	$2^{+}$	9.912	9.929	$0.2$

Table 3. Table 3: The experimentally measured  Tanabashi:2018oca and model calculated heavy flavor baryon masses M E subscript 𝑀 𝐸 M_{E} and M T subscript 𝑀 𝑇 M_{T} . D R subscript 𝐷 𝑅 D_{R} is the relative difference between the data and our calculation.

Baryon	$J^P$	$M_b^E$	$M_b^T$	$D_R$
		(GeV)	(GeV)	(%)
${\Lambda }_c^{+}$	${\left(1/2\right)}^{+}$	2.286	2.440	$6.8$
${\Sigma }_c^{++}$	${\left(1/2\right)}^{+}$	2.454	2.413	$1.6$
${\Sigma }_c^{+}$	${\left(1/2\right)}^{+}$	2.453	2.413	$1.5$
${\Sigma }_c^0$	${\left(1/2\right)}^{+}$	2.454	2.413	$1.6$
${\Xi }_c^{+}$	${\left(1/2\right)}^{+}$	2.468	2.557	$3.6$
${\Xi }_c^0$	${\left(1/2\right)}^{+}$	2.471	2.557	$3.5$
${\Xi }_c^{{}_{}{}^{\prime }+}$	${\left(1/2\right)}^{+}$	2.577	2.566	$0.4$
${\Xi }_c^{{}_{}{}^{\prime }0}$	${\left(1/2\right)}^{+}$	2.579	2.566	$0.5$
${\Omega }_c^0$	${\left(1/2\right)}^{+}$	2.695	2.681	$0.5$
${\Xi }_{cc}^{++}$	${\left(1/2\right)}^{+}$	3.621	3.632	$0.3$
${\Xi }_{cc}^{+}$	${\left(1/2\right)}^{+}$	3.619	3.632	$0.4$
${\Omega }_{cc}^{+}$	${\left(1/2\right)}^{+}$		3.745
${\Sigma }_c^{\ast ++}$	${\left(3/2\right)}^{+}$	2.518	2.429	$3.6$
${\Sigma }_c^{\ast +}$	${\left(3/2\right)}^{+}$	2.518	2.429	$3.6$
${\Sigma }_c^{\ast 0}$	${\left(3/2\right)}^{+}$	2.518	2.429	$3.6$
${\Xi }_c^{\ast +}$	${\left(3/2\right)}^{+}$	2.646	2.567	$3.0$
${\Xi }_c^{\ast 0}$	${\left(3/2\right)}^{+}$	2.646	2.567	$3.0$
${\Omega }_c^{\ast 0}$	${\left(3/2\right)}^{+}$	2.766	2.689	$2.8$
${\Xi }_{cc}^{\ast ++}$	${\left(3/2\right)}^{+}$		3.644
${\Xi }_{cc}^{\ast +}$	${\left(3/2\right)}^{+}$		3.644
${\Omega }_{cc}^{\ast +}$	${\left(3/2\right)}^{+}$		3.754
${\Omega }_{ccc}^{++}$	${\left(3/2\right)}^{+}$		4.784
${\Lambda }_b^0$	${\left(1/2\right)}^{+}$	5.620	5.793	$3.1$
${\Sigma }_b^{+}$	${\left(1/2\right)}^{+}$	5.811	5.769	$0.7$
${\Sigma }_b^0$	${\left(1/2\right)}^{+}$		5.769
${\Sigma }_b^{-}$	${\left(1/2\right)}^{+}$	5.816	5.769	$0.8$
${\Xi }_b^0$	${\left(1/2\right)}^{+}$	5.792	5.913	$2.1$
${\Xi }_b^{-}$	${\left(1/2\right)}^{+}$	5.795	5.913	$2.0$
${\Xi }_b^{{}_{}{}^{\prime }0}$	${\left(1/2\right)}^{+}$	5.792	5.903	$1.9$
${\Xi }_b^{{}_{}{}^{\prime }-}$	${\left(1/2\right)}^{+}$	5.795	5.903	$1.9$
${\Omega }_b^{-}$	${\left(1/2\right)}^{+}$	6.046	6.021	$0.4$
${\Xi }_{bb}^{-}$	${\left(1/2\right)}^{+}$		10.210
${\Xi }_{bb}^0$	${\left(1/2\right)}^{+}$		10.210
${\Omega }_{bb}^{-}$	${\left(1/2\right)}^{+}$		10.319
${\Sigma }_b^{\ast +}$	${\left(3/2\right)}^{+}$	5.832	5.781	$0.9$
${\Sigma }_b^{\ast 0}$	${\left(3/2\right)}^{+}$		5.7 81
${\Sigma }_b^{\ast -}$	${\left(3/2\right)}^{+}$	5.835	5.781	$0.9$
${\Xi }_b^{\ast 0}$	${\left(3/2\right)}^{+}$		5.915
${\Xi }_b^{\ast -}$	${\left(3/2\right)}^{+}$		5.915
${\Omega }_b^{\ast -}$	${\left(3/2\right)}^{+}$		6.033
${\Xi }_{bb}^{\ast +}$	${\left(3/2\right)}^{+}$		10.221
${\Xi }_{bb}^{\ast 0}$	${\left(3/2\right)}^{+}$		10.221
${\Omega }_{bb}^{\ast -}$	${\left(3/2\right)}^{+}$		10.331
${\Omega }_{bbb}^{-}$	${\left(3/2\right)}^{+}$		14.499

Table 4. Table 4: The experimentally measured and model calculated ground and excited states of Ω c subscript Ω 𝑐 \Omega_{c} .

	Experiment		Model
Baryon	$J^P$	$M_E$	$J^P$	$M_T$	$D_R$
		(MeV)		(MeV)	(%)
${\Omega }_c^0$	${\left(1/2\right)}^{+}$	2.695	${\left(1/2\right)}^{+}$(1S)	2.681	$0.5$
${\Omega }_c^{\ast }{\left(2770\right)}^0$	${\left(3/2\right)}^{+}$	2.766	${\left(3/2\right)}^{+}$(1S)	2.689	$2.8$
${\Omega }_c{\left(3000\right)}^0$		3.000	${\left(1/2\right)}^{-}$(1P)	2.990	$0.3$
${\Omega }_c{\left(3050\right)}^0$		3.050	${\left(3/2\right)}^{-}$(1P)	3.052	$0.1$
${\Omega }_c{\left(3065\right)}^0$		3.065	${\left(1/2\right)}^{-}$(1P)	3.074	$0.3$
${\Omega }_c{\left(3090\right)}^0$		3.090	${\left(3/2\right)}^{-}$(1P)	3.085	$0.2$
${\Omega }_c{\left(3120\right)}^0$		3.119	${\left(5/2\right)}^{-}$(1P)	3.252	$4.3$

Table 5. Table 5: The scaled dissociation temperatures T d / T c subscript 𝑇 𝑑 subscript 𝑇 𝑐 T_{d}/T_{c} for charmed mesons and baryons in the two limits of V = F 𝑉 𝐹 V=F and V = U 𝑉 𝑈 V=U .

	$J/\psi$	${\chi }_c$	${\psi }^{\prime }$	$D_s$	$D_s^{\ast }$	$D^0$	$D^{\ast 0}$
$V=F$	1.42	-	-	1.14	1.10	1.10	1.08
$V=U$	3.09	1.30	1.24	2.50	1.98	2.35	1.80
	${\Omega }_{ccc}$	${\Omega }_{cc}$	${\Xi }_{cc}$	${\Omega }_c$	${\Xi }_c$	${\Lambda }_c$
$V=F$	1.15	1.06	1.05	1.03	1.02	1.02
$V=U$	2.18	1.63	1.54	1.41	1.39	1.37

Equations

$$\left[\sum _{i=1}^3\frac{{\mathbf{p}}_i^2}{2{ϵ}_i}+\sum _{i<j}^3\frac{{ϵ}_i+{ϵ}_j}{2{ϵ}_i{ϵ}_j}{\mathcal{V}}_{ij}\right]\Psi =E\Psi ,$$

$$M=3{ϵ}_i+\sum _{j\ne i}\frac{m_i^2-m_j^2}{{ϵ}_i+{ϵ}_j}.$$

$${\mathcal{V}}_{ij}$$

$$\mathbf{R}$$

$$\mathbf{\rho }$$

$$\mathbf{\lambda }$$

$$|n_{\rho }{l}_{\rho }{m}_{\rho }{n}_{\lambda }{l}_{\lambda }{m}_{\lambda }⟩=|n_{\rho }{l}_{\rho }{m}_{\rho }⟩|n_{\lambda }{l}_{\lambda }{m}_{\lambda }⟩$$

$$\Phi (|{\mathbf{r}}_{ij}|)=\frac{{ϵ}_i^2{ϵ}_j^2{\omega }^2}{({ϵ}_i+{ϵ}_j)({ϵ}_1+{ϵ}_2+{ϵ}_3)}{r}_{ij}^2.$$

$$|{\Psi }_{FSC}⟩=|F⟩\times |S⟩\times |n_{\rho }{l}_{\rho }{m}_{\rho }{n}_{\lambda }{l}_{\lambda }{m}_{\lambda }⟩.$$

$$|\Psi ⟩=\sum _{FSC}{C}_{FSC}|{\Psi }_{FSC}⟩.$$

$$H_{FSC,F^{\prime }{S}^{\prime }{C}^{\prime }}=⟨{\Psi }_{FSC}|H|{\Psi }_{F^{\prime }{S}^{\prime }{C}^{\prime }}⟩.$$

$$\sum _{F^{\prime }{S}^{\prime }{C}^{\prime }}{H}_{FSC,F^{\prime }{S}^{\prime }{C}^{\prime }}{C}_{F^{\prime }{S}^{\prime }{C}^{\prime }}=E{C}_{FSC}$$

$$H_{FSC,F^{\prime }{S}^{\prime }{C}^{\prime }}$$

$${\Psi }_{n_{\rho }{l}_{\rho }{m}_{\rho }}(\mathbf{\rho }){\Psi }_{n_{\lambda }{l}_{\lambda }{m}_{\lambda }}(\mathbf{\lambda })$$

$$\frac{\partial E}{\partial \alpha }=0,\frac{{\partial }^{2}E}{\partial {\alpha }^2}>0,$$

$$A_{q\bar{q}}(r,T)$$

$$S_{q\bar{q}}(r,T)$$

$$ϵ(T_d)$$

$$⟨r(T_d)⟩$$

$$ϵ(T)=M(\infty ,T)-M(T).$$

$$M(\infty ,T)=V_{q\bar{q}}(\infty ,T)+\sqrt{V_{q\bar{q}}^2(\infty ,T)+(m_1+m_2{)}^2}.$$

$$\frac{1}{6}\sum _{ij}\frac{{ϵ}_j^2-m_j^2}{{ϵ}_i}=\sum _{i<j}\frac{{ϵ}_i+{ϵ}_j}{2{ϵ}_i{ϵ}_j}{\mathcal{V}}_{ij}(\infty ,T)$$

$${\mathcal{V}}_{ij}(\infty ,T)=\frac{2m_i{m}_j}{{ϵ}_i+{ϵ}_j}{V}_{qq}(\infty ,T)+V_{qq}^2(\infty ,T)$$

$${ϵ}_i=\frac{M(\infty ,T)}{3}+\frac{1}{3}\sum _{j\ne i}\frac{m_i^2-m_j^2}{{ϵ}_i+{ϵ}_j}.$$

$$-⟨n_{\rho }^{\prime }{l}_{\rho }^{\prime }{m}_{\rho }^{\prime }{n}_{\lambda }^{\prime }{l}_{\lambda }^{\prime }{m}_{\lambda }^{\prime }|{\nabla }_a^2|n_{\rho }{l}_{\rho }{m}_{\rho }{n}_{\lambda }{l}_{\lambda }{m}_{\lambda }⟩$$

$$⟨n_{\rho }^{\prime }{l}_{\rho }^{\prime }{m}_{\rho }^{\prime }{n}_{\lambda }^{\prime }{l}_{\lambda }^{\prime }{m}_{\lambda }^{\prime }|\Phi (|{\mathbf{r}}_{12}|)|n_{\rho }{l}_{\rho }{m}_{\rho }{n}_{\lambda }{l}_{\lambda }{m}_{\lambda }⟩={\delta }_{l_{\rho }}^{l_{\rho }^{\prime }}{\delta }_{m_{\rho }}^{m_{\rho }^{\prime }}{\delta }_{n_{\lambda }}^{n_{\lambda }^{\prime }}{\delta }_{l_{\lambda }}^{l_{\lambda }^{\prime }}{\delta }_{m_{\lambda }}^{m_{\lambda }^{\prime }}\int {\rho }^2\mathrm{d}\rho {\psi }_{n_{\rho }{l}_{\rho }}(\rho ){\psi }_{n_{\rho }^{\prime }{l}_{\rho }^{\prime }}(\rho )\Phi \left(\sqrt{\frac{\overline{m}({ϵ}_1+{ϵ}_2)}{{ϵ}_1{ϵ}_2}}\rho \right).$$

$${|F⟩}_{S/A}$$

$${|F⟩}_{M1S/A}$$

$${|F⟩}_{M2S/A}$$

$${|S⟩}_{S1/2}$$

$${|S⟩}_{S3/4}$$

$${|S⟩}_{MS1/2}$$

$${|S⟩}_{MA1/2}$$

$${\hat{r}}_z|l,m⟩$$

$${\hat{r}}^{\pm }|l,m⟩$$

Full text

I Introduction

Heavy quark masses are much larger than the typical temperature of the fireball formed in high energy nuclear collisions at Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) and the typical energy scale of the Quantum Chromodynamics (QCD), $m_{Q}\sim(1-5)$ GeV $\gg T\sim(0.3-0.5)$ GeV and $m_{Q}\gg\Lambda_{QCD}$, their thermal production in the fireball can then be safely neglected and they are entirely originated from the initial hard processes which can be solidly calculated in the perturbation theory. On the way out of the fireball, the heavy quarks pass through the quark-gluon plasma (QGP) created in the collisions and strongly interact with it. Since heavy flavor hadrons inherit the information of the QGP carried by the heavy quarks, they have long been considered as a sensitive probe of the QGP formation Matsui:1986dk .

For closed heavy flavors, one can neglect the parton creation-annihilation fluctuations and describe their static properties with non-relativistic potential models isgur . The study is directly extended from vacuum to hot medium, by taking heavy quark potential from lattice QCD at finite temperature Petreczky:2010yn ; Kaczmarek:1900zz . Solving the two-body Schrödinger equation for a $Q\bar{Q}$ system, people obtained the sequential dissociation temperatures for charmonium and bottominum states Satz:2005hx ; Burnier:2015tda . The potential model is also extended to studying open heavy flavor states like $D$ mesons and $\Xi_{cc}$ and $\Omega_{ccc}$ baryons, by using the two-body Dirac equation Crater:2002fq ; Crater:2008rt ; Guo:2012hx ; Shi:2013rga or three-body Schrödinger equation SilvestreBrac:1996bg ; He:2014tga ; Zhao:2016ccp ; Zhao:2017gpq .

For heavy flavor baryons with a typical radius $0.5$ fm, the uncertainty principle leads to a characteristic momentum around $0.4$ GeV, corresponding to a relative velocity $v\sim c$ between any two constituents. Therefore, the relativistic effect should seriously be taken into account, especially for singly charmed baryons. On the other hand, the spin interaction, which results in the fine structure of hadrons and becomes more important in multi-quark states, is a relativistic effect and should self-consistently be included through the covariant Dirac equation. In this paper, we study the three-body Dirac equation for heavy flavor baryons in hot medium and determine the dissociation temperatures for all the heavy flavor hadrons in a self-consistent way.

The covariant wave equation proposed by Sazdjian Sazdjian:1988be ; Sazdjian:1986aw provides a way to find relativistic bound-state wave functions. The interaction potential and the wave function of the bound-state are related in a definite way to the kernel and the wave function of the Bethe-Salpeter equation Salpeter:1951sz ; Sazdjian:1986aw . Crater and Alstine developed the two-body Dirac equation Crater:2002fq ; Crater:1983ew ; Crater:1987hm ; Crater:2008rt for mesons which covariantly converts the Dirac equation with 16 degrees of freedom to a Klein–Gordon-like equation with highly disentangled degrees of freedom. In this framework, one can clearly see the relativistic corrections to the non-relativistic potential model through the Darwin term and many spin interaction terms. Later on, the three-body Dirac equation for baryons was developed by Whitney and Crater Whitney:2011aa , which combines Sazdjian’s $N$-body relativistic potential model with the two-body Dirac equation.

To systematically describe heavy flavor mesons and baryons in a relativistic potential model at finite temperature, we first improve the method to solve the three-body Dirac equation. We will expand the baryon wave functions in terms of the two-body spherical harmonic oscillator states, which guarantees the completeness and orthogonality of the expansion, increases the accuracy of the calculation, and can be used to study both the ground and excited baryon states. We will also take a universal set of model parameters for both mesons and baryons. This makes the calculation more self-consistent and predictable. With the improved method, we then solve the two- and three-body Dirac equations in hot medium, and determine the heavy flavor dissociation temperatures by calculating their binding energies and averaged radii.

II Three-body Dirac equation

The two-body Dirac equation which is used to describe heavy flavor mesons is discussed in detail in Refs. Crater:2002fq ; Crater:2008rt in vacuum and Refs. Guo:2012hx ; Shi:2013rga at finite temperature. We will focus in this section on the three-body Dirac equation for heavy flavor baryons. In the framework firstly developed by Whitney et al Whitney:2011aa based on Sazdjian’s $N$-body relativistic potential model Sazdjian:1988be , the baryon wave function $\Psi({\bf r}_{1},{\bf r}_{2},{\bf r}_{3})$ is controlled by the Schrödinger-like equation,

[TABLE]

where ${\bf r}_{i}$ and ${\bf p}_{i}$ are the quark coordinate and momentum, $E=\sum_{ij}(\epsilon_{j}^{2}-m_{j}^{2})/(6\epsilon_{i})$ is the energy eigenvalue related to the effective quark mass $\epsilon_{i}$ and vacuum quark mass $m_{i}$, and the baryon mass $M$ is determined by the coupled equations,

[TABLE]

The interaction between (anti-)quarks is described by the relativistic potential ${\cal V}_{ij}({\bf r}_{ij})$Crater:2008rt ,

[TABLE]

with the mass and energy parameters $m_{ij}=m_{i}m_{j}/(\epsilon_{i}+\epsilon_{j})$ and $\epsilon_{ij}=((\epsilon_{i}+\epsilon_{j})^{2}-m_{i}^{2}-m_{j}^{2})/(2(\epsilon_{i}+\epsilon_{j}))$, relative coordinate ${\bf r}_{ij}={\bf r}_{i}-{\bf r}_{j}$, unit vector $\hat{\bf r}_{ij}={\bf r}_{ij}/|{\bf r}_{ij}|$, and orbital and spin angular momenta ${\bf L}_{ij}$ and ${\boldsymbol{\sigma}}_{i}$. The non-relativistic central potential $V_{q\bar{q}}(r)$ between quark and antiquark can be separated into two parts, $V_{q\bar{q}}(r)=A_{q\bar{q}}(r)+S_{q\bar{q}}(r)$ with $A_{q\bar{q}}(r)=-\alpha_{q\bar{q}}/r$ and $S_{q\bar{q}}(r)=\sigma_{q\bar{q}}r$, where $A_{q\bar{q}}$ and $S_{q\bar{q}}$ control, respectively, the behavior of the potential at short and long distances. In vacuum one usually takes the Cornell potential Eichten:1978tg , including a Coulomb part which dominates the wave function around $r=0$ and a linear part which leads to the quark confinement. For $V_{qq}$ between two quarks one can still take the separation, $V_{qq}(r)=A_{qq}(r)+S_{qq}(r)$ with $A_{qq}(r)=-\alpha_{qq}/r$ and $S_{qq}(r)=\sigma_{qq}r$. The two couplings $\alpha_{qq}$ and $\sigma_{qq}$ in quark-quark potential are in general very different from $\alpha_{q\bar{q}}$ and $\sigma_{q\bar{q}}$ in quark-antiquark potential. However, considering only one gluon exchange leads to $\alpha_{qq}=\alpha_{q\bar{q}}/2$ Schwartz:2013pla , and the recent lattice simulation Alexandrou:2002sn shows $\sigma_{qq}\approx\sigma_{q\bar{q}}/2$. All the relativistic corrections in the potential (3), including the Darwin term $\Phi_{D}$, spin-spin coupling $\Phi_{SS}$, spin-orbital couplings $\Phi_{SO},\Phi_{SOT},\Phi_{SOD}$ and $\Phi_{SOX}$, and tensor coupling $\Phi_{T}$, are explicitly shown in Ref.Crater:2008rt . Note that, any $\Phi$ is a function of the distance $|{\bf r}_{ij}|$, the dependence of the potential ${\cal V}_{ij}$ on the azimuthal angles is reflected in the coefficients in front of $\Phi$.

Solving the baryon mass $M$ and wave function $\Psi$ from the Schrödinger-like equation (1) is not a simple eigenstate problem, as the eigenvalue $E$ or $M$ appears also on the left hand side of the equation. We solve $M$ or the binding energy $\mathcal{E}=M-\widetilde{M}\ (\widetilde{M}=m_{1}+m_{2}+m_{3})$ by employing iteration method. For a given value $M^{(n)}$ at the $n$-th step, we first obtain the corresponding effective quark mass $\epsilon_{i}^{(n)}$ from the coupled equations (2), and then solve the Schrödinger-like equation (1) as an eigenstate problem with eigenvalue $E^{(n)}$. With sufficiently large number of iteration steps, one could obtain the baryon mass up to arbitrary accuracy.

To solve the Schrödinger-like equation for multi-body bound state problem, a usually used way is to expand the wave function in terms of known functions. While in some special cases, for instance the ground state, one can use variational method and expand the state in terms of Gaussian wave packets Hiyama:2018ivm , a general and systematic study including both ground and excited states should be carried out in a complete and orthogonal Hilbert space. We employ a numerical framework similar to Ref. Capstick:1986bm to treat the three-body bound state problem. We solve the baryon mass and wave function in the Hilbert space constructed by spherical harmonic oscillator states which are by definition complete and orthogonal. Considering the mass difference among the three quarks in a general baryon state, we take different constituent masses in constructing the spherical harmonic oscillator states. This is different from Ref. Capstick:1986bm where all the three constituents have the same mass. In such a framework we can study both the ground and excited states.

Like a two-body problem, we factorize the three-body motion into a center-of-mass motion and a relative motion. To this end, we introduce the coordinates

[TABLE]

instead of ${\bf r}_{1}$, ${\bf r}_{2}$ and ${\bf r}_{3}$, where $\overline{m}$ is an arbitrary parameter which automatically disappears in the end and does not affect the result 1998FBS….25..199K . Accordingly, we use the center-of-mass momentum ${\bf P}$ and relative momenta ${\bf p}$ and ${\bf q}$ instead of ${\bf p}_{1}$, ${\bf p}_{2}$ and ${\bf p}_{3}$.

We now focus on the relative motion of the three-body problem which controls the inner structure of the bound state. We expand the relative wave function in terms of two-body spherical harmonic oscillator states. For a single spherical harmonic oscillator, its Schrödinger equation with potential $\mu\omega^{2}r^{2}/2$ can be straightforwardly solved with solutions $\Psi_{nlm}({\bf r})=\psi_{nl}(r)Y_{l}^{m}(\theta,\phi)$ and $E_{nl}=(2n+l+3/2)\hbar\omega$, where $n,l,m$ are the principal, orbital and magnetic quantum numbers, and $Y_{l}^{m}$ and $\psi_{nl}$ are the angular and radial parts of the wave function. Note that the radial wave function $\psi_{nl}$ does not depend on the mass $\mu$ and frequency $\omega$ separately, it depends on the combined quantity $\alpha=\sqrt{\mu\omega/\hbar}$.

The two-body spherical harmonic oscillator states are defined as a direct product of two single spherical harmonic oscillator states,

[TABLE]

with the same scaling parameter $\alpha$ for the two oscillators.

The above defined two-body spherical harmonic oscillator states are the exact solutions of a three-body bound state problem with interaction potential

[TABLE]

With the two-body spherical harmonic oscillator states obtained, we define the Hilbert space as a direct product of the flavor, spin and coordinate spaces,

[TABLE]

The possible flavor and spin states are shown in Appendix A. In such a Hilbert space we expand the baryon wave function

[TABLE]

Correspondingly, the Hamiltonian operator $H$ of the system becomes a matrix in the Hilbert space,

[TABLE]

Taking into account the complete and orthogonal conditions for the states $\left|\Psi_{FSC}\right>$, the eigenstate problem of a three-body system, $H\left|\Psi\right>=E\left|\Psi\right>$, becomes a matrix equation for the coefficients $C_{FSC}$,

[TABLE]

with the Hamiltonian matrix elements

[TABLE]

The calculation for the kinetic energy $K=\sum_{i}p_{i}^{2}/(2\epsilon_{i})$, which is proportional to $\nabla_{\rho}^{2}+\nabla_{\lambda}^{2}$, is straightforward. The difficulty comes from the calculation for the potential part $W=\sum_{i<j}(\epsilon_{i}+\epsilon_{j}){\cal V}_{ij}({\bf r}_{ij})/(2\epsilon_{i}\epsilon_{j})$. For the potential ${\cal V}_{12}({\bf r}_{12})$ between the first and second constituents, since ${\bf r}_{12}$ is proportional to $\boldsymbol{\rho}$ and independent of $\boldsymbol{\lambda}$, the six dimensional integration $\int\mathrm{d}^{3}\boldsymbol{\rho}\,\mathrm{d}^{3}\boldsymbol{\lambda}$ is reduced to a simple radial integration $\int\rho^{2}\mathrm{d}\rho$. For the other two potentials ${\cal V}_{23}({\bf r}_{23})$ and ${\cal V}_{13}({\bf r}_{13})$, since ${\bf r}_{23}$ and ${\bf r}_{13}$ depend on both $\boldsymbol{\rho}$ and $\boldsymbol{\lambda}$, the integration over $\boldsymbol{\rho}$ and $\boldsymbol{\lambda}$ can not be so simplified. We then make a rotation in the coordinate space from $(\boldsymbol{\rho},\boldsymbol{\lambda})$ to $(\widetilde{\boldsymbol{\rho}},\widetilde{\boldsymbol{\lambda}})$ to guarantee ${\bf r}_{23}\sim\widetilde{\boldsymbol{\rho}}$ or ${\bf r}_{13}\sim\widetilde{\boldsymbol{\rho}}$. This rotation in coordinate space is equivalent to a particle index rotation from the ordering $\{1,2,3\}$ to $\{2,3,1\}$ or $\{3,1,2\}$. The rotation can be explicitly written as

[TABLE]

where $D$ is the rotation matrix. After the rotation, the six dimensional integration in calculating the expectation values of ${\cal V}_{23}$ and ${\cal V}_{13}$ becomes again a simple radial integration $\int{\widetilde{\rho}}^{2}\mathrm{d}\widetilde{\rho}$. The Hamiltonian matrix elements are explicitly listed in Appendix A.

III Heavy flavors in vacuum

Before we solve the two- and three-body Dirac equations at finite temperature, we first fix the model parameters by fitting the heavy flavor masses in vacuum. Note that, the completeness of the two-body oscillator states $\left|n_{\rho}l_{\rho}m_{\rho}n_{\lambda}l_{\lambda}m_{\lambda}\right>$ is independent of the value of the scaling parameter $\alpha$, if we take infinite number of states. When the number of the states is finite in practical computation, the accuracy of the calculation depends on the choice of the $\alpha$ value.

For a baryon ground state, the $\alpha$ value corresponding to the most fast convergence can be fixed by applying variational method. We take $\alpha$ as a parameter and calculate the energy eigenvalue $E$ and wave function $\Psi$ as functions of $\alpha$. By minimizing the eigenvalue

[TABLE]

we obtain the $\alpha$ value corresponding to the most fast convergence.

The number of the oscillator states or the size of the Hilbert space is controlled by the total principal quantum number $N=2n_{\rho}+2n_{\lambda}+l_{\rho}+l_{\lambda}$. We truncate the Hilbert space by taking $N\leq N_{max}=0,2,4,6,8$ and $10$. Considering the degeneracy of the wave functions, the corresponding number of states is $1,28,210,924,3003$ and $8008$. From the parity symmetry, which decouples states with odd and even $N$, and momentum conservation, which decouples states with different total magnetic quantum number $m_{\rho}+m_{\lambda}$, the effective numbers of states are reduced to $1,8,34,108,259$ and $560$. We first consider three specific baryon ground states $\left|ssc\right>$, $\left|scc\right>$ and $\left|ccc\right>$ containing one, two and three charm quarks. We calculate their binding energy $\mathcal{E}=M-\widetilde{M}$ and check its stability when $N_{max}$ is large enough. The result is shown in Figure 1 for the three states and six values of $N_{\max}$. With increasing $N_{\max}$, the $\alpha$-dependence of $\mathcal{E}$ becomes more and more smooth and reaches almost a constant for $N_{\max}=8$ and $10$. For the state $\left|ssc\right>$, the minimum binding energy obtained from the variational method is $\mathcal{E}=0.642,\ 0.630,\ 0.616,\ 0.613,\ 0.610$ and $0.609$ GeV corresponding to $\alpha=0.72,\ 0.72,\ 0.80,\ 0.80,\ 0.84$ and $0.84$ GeV, when $N_{\max}$ increases from [math] to $10$. For the other two states $\left|scc\right>$ and $\left|ccc\right>$, the $\alpha$ dependence is similar, but the binding energy reaches stability much faster. In the following calculation, we will take $N_{\max}=10$ which already guarantees the stability of the numerical calculation.

To self-consistently describe both the meson and baryon states, we take a universal set of parameters, by fitting the heavy flavor meson and baryon masses. The parameters we used, including the vacuum quark masses and coupling strengths $\alpha$ and $\sigma$, are shown in Table 1. Note that, the relations between the couplings $\alpha_{qq}=0.45\,\alpha_{q\bar{q}}$ and $\sigma_{qq}=0.49\,\sigma_{q\bar{q}}$ are close to the results from one gluon exchange and lattice QCD.

The calculated heavy flavor meson mass $M_{T}$ and the comparison with the experimentally measured mass $M_{E}$ Tanabashi:2018oca are shown in Table 2. From the relative difference $D_{R}=|(M_{T}-M_{E})/M_{E}|$ between the data and model calculation, the two-body Dirac equation describes heavy flavor mesons reasonably well, especially for charmonia and all bottomed mesons with relative difference $D_{R}<1\%$. Since the potential description of the interaction between a light and a heavy quark is no longer a good enough approximation, the relative difference $D_{R}$ becomes larger for $D$ mesons, $D_{R}\sim(2-4)\%$.

The baryon masses and the comparison with the experimental data are shown in Table 3. For doubly charmed baryons $\Xi_{cc}^{+}$ and $\Xi_{cc}^{++}$ which have been observed experimentally, the agreement between the calculation and the data is good with relative difference $D_{R}<0.4\%$. Theoretically, we expect that the relativistic potential model should be more reliable in describing triply charmed and doubly and triply bottomed baryons. For singly charmed or bottomed baryons, however, the model calculation is not so good, since the two light quarks in such a baryon make the potential approximation no longer good enough. Especially for those baryons with two light quarks $u$ and $d$, the relative difference becomes much larger. For $\Lambda_{c}^{+}$ it reaches $\sim 7\%$. It is phenomenologically suggested that SilvestreBrac:1996bg , a three body force $C/(m_{1}m_{2}m_{3})$ controlled by the constituent quark masses may improve the calculation for heavy flavor baryons, especially for the baryons with two light quarks.

With the expansion method in a complete and orthogonal Hilbert space, we can calculate not only the baryon ground states but also the excited states. Table 4 shows the result for $\Omega_{c}[ssc]$. $\Omega_{c}^{0}$ and $\Omega_{c}^{*}(2770)^{0}$ are the ground states with $J^{P}=(1/2)^{+}$ and $(3/2)^{+}$, and the other five baryons are excited states. The masses of the excited states have been experimentally measured, but the total spins and parities are not yet fixed. From our model calculation, they are all $1P$ states with $J^{P}=(1/2)^{-},(3/2)^{-}$ and $(5/2)^{-}$. The mass difference among the excited states is mainly from the spin-orbit and spin-spin interactions. The expansion method in our calculation guarantees the accuracy of the calculation, the relative difference for the masses is $D_{R}<4.3\%$.

IV Heavy flavors in the hot medium

We now turn to study the properties of heavy flavor hadrons in the hot medium. Different from light hadrons which are all produced at the end of the quark-gluon plasma, heavy flavor hadrons can be produced in initial hard scattering and survive in the quark-gluon plasma, due to their larger binding energies compared to light hadrons. However, because of the medium screening effect, the interaction between quarks becomes weaker when temperature increases. Above certain temperature, the interaction potential well is too narrow to form some (or all) of the bound states, and these hadrons dissociate. The surviving temperature or the dissociation temperature of a heavy flavor hadron in the hot medium can be determined by solving the corresponding dynamical equation with heavy quark potential at finite temperature. This is the original idea of Matsui and Satz Matsui:1986dk who suggested $J/\psi$ suppression as a signal of the quark-gluon plasma formation in heavy ion collisions. By solving the non-relativistic Schrödinger equation for charmonium states, $\chi_{c}$, $\psi^{\prime}$ and $J/\psi$ are found to be sequentially dissociated with increasing temperature Satz:2005hx . The potential which controls the absolute value of the dissociation temperature depends on the detail of dissociation process in the hot medium Brambilla:2008cx . Yet, its exact form is undetermined on the theoretical side, but can be constrained by two limits. One is the slow dissociation limit, which allows enough time for the constituent quarks to exchange heat with the hot medium and reach thermal equilibrium, and the potential takes the form of the free energy $F$. The other is the rapid dissociation limit, assuming no heat exchange between constituent quarks and the medium, so the potential is the internal energy $U$. The free energy $F$ has been computed via lattice QCD simulation Petreczky:2010yn , while the internal energy $U$ can be obtained via the thermodynamic relation $U=F+TS$ where $S=-\partial F/\partial T$ is the entropy. In a general case in high energy nuclear collisions, the potential $V$ is in between these two limits, $F<V<U$ Satz:2005hx . More discussions on the choice of interaction potential can be found in Refs. Liu:2018syc ; Lee:2013dca . From the above thermodynamic relation one can tell the potential well is deeper in the internal energy limit, hence the dissociation temperature with $V=U$ is higher than that with $V=F$. Considering the Debye screening at finite temperature, the short and long range interaction potential $V_{q\bar{q}}=F_{q\bar{q}}=A_{q\bar{q}}+S_{q\bar{q}}$ can be expressed respectively as Karsch:1987pv ; Satz:2005hx ; Shi:2013rga

[TABLE]

where $\Gamma$ is the Gamma function, $K$ is the modified Bessel function of the second kind, and the temperature dependent screening mass $\mu(T)$ can be extracted by fitting the lattice simulated free energy Petreczky:2010yn . From the known free energy $F$, one can then obtain the other limit of the potential, $V=U$.

From the definition of the dissociation temperature $T_{d}$ for a heavy quark bound state in hot medium, it is determined by the vanishing binding energy or infinite averaged radius,

[TABLE]

Considering the plateau structure of the heavy quark potential at large distance, see the lattice simulated free energy Petreczky:2010yn , the binding energy $\epsilon$ is the difference between the hadron masses calculated with potential $V(r=\infty,T)$ and $V(r,T)$,

[TABLE]

For heavy flavor mesons, $M(T)$ is calculated from the two-body Dirac equation and $M(\infty,T)$ is simply expressed as

[TABLE]

For heavy flavor baryons, $M(T)$ is calculated from the three-body Dirac equation, and $M(\infty,T)$ is determined by the upper limit of the binding energy,

[TABLE]

with

[TABLE]

and

[TABLE]

The scaled meson and baryon binding energy and averaged radius in the two limits of $V=F$ and $V=U$ are shown in Figures 2 and  3 as functions of scaled temperature. From the recent lattice simulation Bazavov:2011nk , the critical temperature for the QCD deconfinement phase transition is $T_{c}\simeq 155$ MeV. With increasing temperature, the binding energy decreases and the averaged radius increases, which indicate a more and more loosely bounded state. Finally, at the dissociation temperature $T_{d}$ the binding energy reaches zero and the averaged radius becomes infinity, which means the vanishing of the bound state. Since heavy flavors are in more deeply bound states in the limit of $V=U$, they can survive at higher temperatures in this case. The surviving temperatures for charmed mesons and baryons are listed in Table 5. It is clear to see a sequential dissociation for both mesons and baryons, due to their different binding energies.

V Summary

We investigated heavy flavor hadrons in vacuum and hot medium in the framework of multi-body Dirac equations. With the improved method to solve the three-body Dirac equation, we first fixed a universe set of model parameters by fitting the heavy flavor meson and baryon masses in vacuum. Taking the free energy and internal energy simulated by lattice QCD as the two limits of the heavy quark potential, we then solved the eigenvalues and eigenstates of the two- and three-body Dirac equations at finite temperature, from which we systematically determined the binding energies, averaged sizes and in turn dissociation temperatures for all the heavy flavors. Mesons and baryons are separately sequentially dissociated in the quark-gluon plasma. In comparison with mesons, baryons are easier to be melted in hot medium due to the weaker interaction between quark-quark pairs than that between quark-antiquark pairs. Considering the limitation of the potential model, the Dirac equations should be systematically improved to better approximate QCD and then used to make better predictions for comparison with the experimental data, especially for mesons with one light quark and baryons with two light quarks.

Acknowledgement: We thank Mr. Tiecheng Guo for the collaboration in the beginning of the work and Profs. Charles Gale and Sangyong Jeon for helpful discussions. JZ and PZ are supported by the NSFC grant Nos. 11575093 and 11890712, and SS is grateful to the Natural Sciences and Engineering Research Council of Canada. Computations were partially made on the supercomputer Belúga, managed by Calcul Québec and Compute Canada.

Appendix A Hamiltonian matrix elements

Let’s first consider the elements of the baryon kinetic energy ($\sim({\bf\nabla}_{\rho}^{2}+{\bf\nabla}_{\lambda}^{2})$) and averaged radius ($\sim(\rho^{2}+\lambda^{2})$), since they can be directly calculated without rotation from $({\boldsymbol{\rho}},{\boldsymbol{\lambda}})$ to $({\widetilde{\boldsymbol{\rho}}},{\widetilde{\boldsymbol{\lambda}}})$. Taking into account the orthogonality of the spherical harmonic oscillator states, $\left<n_{\rho}^{\prime}l_{\rho}^{\prime}m_{\rho}^{\prime}n_{\lambda}^{\prime}l_{\lambda}^{\prime}m_{\lambda}^{\prime}|n_{\rho}l_{\rho}m_{\rho}n_{\lambda}l_{\lambda}m_{\lambda}\right>=\delta_{n_{\rho}}^{n_{\rho}^{\prime}}\,\delta_{l_{\rho}}^{l_{\rho}^{\prime}}\,\delta_{m_{\rho}}^{m_{\rho}^{\prime}}\,\delta_{n_{\lambda}}^{n_{\lambda}^{\prime}}\,\delta_{l_{\lambda}}^{l_{\lambda}^{\prime}}\,\delta_{m_{\lambda}}^{m_{\lambda}^{\prime}}$, we have

[TABLE]

for $a=\rho,\lambda$. Note that, the kinetic energy and averaged radius are diagonal in spin and flavor spaces.

For the potential ${\cal V}_{ij}({\bf r}_{ij})$ between quark $i$ and quark $j$, the spin-independent terms and spin-spin coupling term depend only on the distance between two quarks. For $\Phi(|{\bf r}_{12}|)$, we have simply

[TABLE]

After the rotation (II), we have similar matrix elements for $\Phi(|{\bf r}_{13}|)$ and $\Phi(|{\bf r}_{23}|)$ with the replacement of $({\boldsymbol{\rho}},{\boldsymbol{\lambda}})$ by $({\widetilde{\boldsymbol{\rho}}},{\widetilde{\boldsymbol{\lambda}}})$.

The orbital angular momentum dependent terms and tensor term in ${\cal V}_{ij}({\bf r}_{ij})$ depend on the azimuthal angles $\theta$ and $\phi$ of ${\bf r}_{ij}$. Note that, the couplings $\Phi_{SO},\Phi_{SOT},\Phi_{SOD},\Phi_{SOX}$ and $\Phi_{T}$ depend only on the distance $|\hat{\bf r}_{ij}|$, the angle dependence comes from their coefficients $H({\boldsymbol{\sigma}}_{i},{\boldsymbol{\sigma}}_{j},{\bf L}_{ij},\hat{\bf r}_{ij})$ shown in (3). To calculate the matrix elements of the coefficients in flavor and spin spaces, we first construct all possible symmetric, antisymmetric, mixed symmetric and mixed antisymmetric flavor states $\left|F\right>$,

[TABLE]

where the symmetry of $\left|F\right>_{M1S}$ and $\left|F\right>_{M2A}$ and antisymmetry of $\left|F\right>_{M2S}$ and $\left|F\right>_{M1A}$ are in the sense of the flavor exchange between the first two quarks $q_{1}\leftrightarrow q_{2}$, corresponding to isospin triplet and singlet for $q_{1},q_{2}=u,d$. We then express all possible three-body spin states $\left|S\right>$ as product of two-body spin state $\left|S,S_{z}\right>_{ij}$ and single spin state $\left|s_{z}\right>_{k}$,

[TABLE]

where one can clearly see that, the mixed antisymmetric states correspond to the spin singlet of the first two quarks, and the symmetric and mixed symmetric states correspond to all the triplets.

Employing the algebraic method used in quantum mechanics for the single spin operator ${\bf s}_{i}$ and orbital angular momentum ${\bf L}_{ij}$, we introduce the ladder operators for the relative coordinate $\hat{\bf r}_{ij}$, $\hat{r}_{z}$ and ${\hat{r}}^{\pm}=({\hat{r}}_{x}\pm i{\hat{r}}_{y})/\sqrt{2}=\sin\theta e^{\pm i\phi}/\sqrt{2}$ with the raising and lowering relations,

[TABLE]

With the above preparation we are now ready to compute the matrix elements of the coefficients $H(\boldsymbol{\sigma_{i}},\boldsymbol{\sigma_{j}},{\bf L}_{ij},{\hat{\bf r}}_{ij})$. For a fixed spin state $\left|S,S_{z}\right>\left|s_{z}\right>$ and a fixed spherical harmonic oscillator $\left|n_{\rho}l_{\rho}m_{\rho}n_{\lambda}l_{\lambda}m_{\lambda}\right>$, the matrix elements $H_{LS,L^{\prime}S^{\prime}}=\left<l^{\prime}_{\rho}m^{\prime}_{\rho}\right|\left<S^{\prime},S^{\prime}_{z}\right|H({\boldsymbol{\sigma}}_{i},{\boldsymbol{\sigma}}_{j},{\bf L}_{ij},{\hat{\bf r}}_{ij})\left|S,S_{z}\right>\left|l_{\rho}m_{\rho}\right>$ in angular momentum space can be analytically calculated. They are

[TABLE]

for spin-spin coupling with $H={\boldsymbol{\sigma}}_{i}\cdot{\boldsymbol{\sigma}}_{j}$,

[TABLE]

with $D_{l}^{m}=\sqrt{(l+m+1)(l-m)}$ for spin-orbital coupling with $H={\bf L}_{ij}\cdot({\boldsymbol{\sigma}}_{i}+{\boldsymbol{\sigma}}_{j})$,

[TABLE]

for spin-orbital difference and cross couplings with $H={\bf L}_{ij}\cdot({\boldsymbol{\sigma}}_{i}-{\boldsymbol{\sigma}}_{j})$ and $H=i{\bf L}_{ij}\cdot({\boldsymbol{\sigma}}_{i}\times{\boldsymbol{\sigma}}_{j})$,

[TABLE]

for tensor coupling with $H=3({\boldsymbol{\sigma}}_{i}\cdot{\hat{\bf r}}_{ij})({\boldsymbol{\sigma}}_{j}\cdot{\hat{\bf r}}_{ij})-{\boldsymbol{\sigma}}_{i}\cdot{\boldsymbol{\sigma}}_{j}$, and

[TABLE]

for spin-orbital tensor coupling with $H=({\boldsymbol{\sigma}}_{i}\cdot{\hat{\bf r}}_{ij})({\boldsymbol{\sigma}}_{j}\cdot{\hat{\bf r}}_{ij})({\bf L}_{ij}\cdot({\boldsymbol{\sigma}}_{i}+{\boldsymbol{\sigma}}_{j}))$.