Application of the Thomas Fermi Quark Model to Multiquark Mesons

TL;DR

This paper applies the Thomas-Fermi quark model to investigate the stability and properties of multiquark mesons with multiple quark-antiquark pairs, revealing potential stability for certain multi-$Z$ mesons.

Contribution

It introduces a novel application of the Thomas-Fermi model to multiquark mesons, deriving energy expressions and analyzing stability for heavy-light quark systems.

Findings

Multi-$Z$ mesons may be stable for certain quark numbers.

Energy trends suggest stability of multi-charmonium and multi-bottomonium systems.

Quark interactions are primarily between quarks and antiquarks, simplifying the model.

Abstract

The possibility of the existence of mesons with two or more quark-antiquark pairs is investigated with a new application of the Thomas-Fermi (TF) statistical quark model. Quark color couplings are treated in a mean field manner similar to a previous application to baryons, and short and concise expressions for energies are derived. We find that, on average, quarks only interact with antiquarks in such systems. The TF differential equation is constructed and systems with heavy-light quark content are examined. Three types of mesonic systems are defined. In the case of charm quarks, multi-charmonium, multi-$Z$ meson and multi-$D$ meson family types are examined. System analogs for bottom quarks are also constructed. Quantitative trends for system energies of mesonic quark matter are extracted as a function of the number of quark pairs. We find indications from energy plots that multi-$Z$ type mesons (and their bottom quark analogs) are actually stable for a range of quark number pairs. At this initial stage we have not yet included explicit spin interaction couplings between quarks, but we can take one level of degeneracy into account in our two-inequivalent TF function construction.

Tables

Table 1. Table 1: Color singlet counting and representation.

Number of pockets	Type of pocket	Representation
$x$	$r\overline{r}$	$R$
$y$	$b\overline{b}$	$B$
$z$	$g\overline{g}$	$G$

Table 2. Table 2: Color interactions for pockets of the same color, pockets of different color and within a pocket.

$\mathrm{𝐑𝐑}$	$𝐫$	$\overline{𝐫}$
$𝐫$	$rr$	$r\overline{r}$
$\overline{𝐫}$	$\overline{r}r$	$\overline{r}\overline{r}$

Table 3. Table 3: Color interaction counting and indices.

Index(i)	Interaction type	${\left(RR\right)}_i$ pocket	${\left(RB\right)}_i$ pocket	${\left(R\right)}_i$ pocket
1	CCR	1	0	0
2	CCA	0	1	0
3	CAR	0	2	0
4	CAA	2	0	1
5	AAR	1	0	0
6	AAA	0	1	0

Table 4. Table 4: The coupling strength ( C i subscript 𝐶 𝑖 C_{i} ), total outcome ( E i subscript 𝐸 𝑖 E_{i} ) and probabilities ( P i subscript 𝑃 𝑖 P_{i} ) for all six interactions.

${\left(C\right)}_i$	Couplings	${\left(E\right)}_i$	Outcomes	${\left(P\right)}_i$	Probabilities
$C_1$	$\frac{4}{3}{g}^2$	$E_1$	$3^{\eta -1}\frac{\eta \left(\eta -1\right)}{2}$	$P_1$	$\frac{\eta -1}{6\left(2\eta -1\right)}$
$C_2$	$-\frac{2}{3}{g}^2$	$E_2$	$3^{\eta -1}\eta \left(\eta -1\right)$	$P_2$	$\frac{\eta -1}{3\left(2\eta -1\right)}$
$C_3$	$\frac{2}{3}{g}^2$	$E_3$	$3^{\eta -1}2\eta \left(\eta -1\right)$	$P_3$	$\frac{2\left(\eta -1\right)}{3\left(2\eta -1\right)}$
$C_4$	$-\frac{4}{3}{g}^2$	$E_4$	$3^{\eta -1}\eta \left(\eta +2\right)$	$P_4$	$\frac{\eta +2}{3\left(2\eta -1\right)}$
$C_5$	$-\frac{4}{3}{g}^2$	$E_5$	$3^{\eta -1}\frac{\eta \left(\eta -1\right)}{2}$	$P_5$	$\frac{\eta -1}{6\left(2\eta -1\right)}$
$C_6$	$-\frac{2}{3}{g}^2$	$E_6$	$3^{\eta -1}\eta \left(\eta -1\right)$	$P_6$	$\frac{\eta -1}{3\left(2\eta -1\right)}$

Table 5. Table 5: The coupling constants and probabilities for certain types of quark and antiquark interactions in mesons.

New symbol	Old symbol	Probability values
$P_{ii}$	$\frac{P_1}{3}$	$\frac{\left(\eta -1\right)}{18\left(2\eta -1\right)}$
	$\frac{P_2}{3}$	$\frac{\eta -1}{9\left(2\eta -1\right)}$
	$\frac{P_3}{3}$	$\frac{2\left(\eta -1\right)}{9\left(2\eta -1\right)}$
${\overline{P}}_{ii}$	$\frac{P_4}{3}$	$\frac{\left(\eta +2\right)}{9\left(2\eta -1\right)}$
${\overline{\overline{P}}}_{ii}$	$\frac{P_5}{3}$	$\frac{\left(\eta -1\right)}{18\left(2\eta -1\right)}$
	$\frac{P_6}{3}$	$\frac{\left(\eta -1\right)}{9\left(2\eta -1\right)}$

Table 6. Table 6: Comparison between masses used for fitting and those obtained after fitting.

Fitted particle	Masses after fit(MeV)	Masses used for fit(MeV)
Case 1-charm	$3049$	$3069$
Case 2-charm	$4015$	$3990$
Case 3-charm	$1960$	$1973$
Case 1-bottom	$9469$	$9445$
Case 2-bottom	$10653$	$10641$
Case 3-bottom	$5239$	$5313$

Equations

$$V_F\equiv \frac{4}{3}\pi p_F^3\left(\vec{r}\right).$$

$$\Delta V_{ph}\equiv V_F\Delta V=\frac{4}{3}\pi p_F^3(\vec{r})\Delta V.$$

$$\Delta N\equiv \frac{g^I}{h^3}\Delta V_{ph}=\frac{4\pi g^I}{3h^3}{p}_F^3(\vec{r})\Delta V.$$

$$\Delta N\equiv n_q(\vec{r})\Delta V,$$

$$n_q(\vec{r})=\frac{4\pi g^I}{3h^3}{p}_F^3(\vec{r}).$$

$$F_{(\vec{r})}(p)dp=\{\begin{array}{l}\frac{4\pi p^{2}dp}{\frac{4}{3}\pi p_F^3(\vec{r})},\text{if}p\le p_F(\vec{r}),\\ 0,\text{otherwise}.\end{array}$$

$$t(\vec{r})={\int }^{p_F}\frac{p^2}{2m_I}{n}_q(\vec{r})F_{(\vec{r})}(p)dp=C_F{\left(n_q(\vec{r})\right)}^{5/3},$$

$$C_F=\frac{{\left(6{\pi }^2{\hslash }^3\right)}^{5/3}}{20{\pi }^2{\hslash }^3{m}_I{\left(g^I\right)}^{2/3}}.$$

$$T=\int t(\vec{r})d^{3}r.$$

$$U=V_N\int \int \frac{n_q(\vec{r})n_q(\vec{r}{}^{\prime })}{|\vec{r}-\vec{r}{}^{\prime }|}{d}^{3}r{d}^3{r}^{\prime },$$

$$E$$

$$=C_F\int {\left(n_q(\vec{r})\right)}^{5/3}{d}^{3}r+V_N\int \int \frac{n_q(\vec{r})n_q(\vec{r}{}^{\prime })}{|\vec{r}-\vec{r}{}^{\prime }|}{d}^{3}r{d}^3{r}^{\prime }.$$

$$E_i=\sum _{x=0}^{\eta }\sum _{y=0}^{\eta -x}\frac{\eta !}{x!y!z!}\times$$

$$\left\{\frac{x(x-1)+y(y-1)+z(z-1)}{2}{(RR)}_i+\left(xy+yz+zx\right){(RB)}_i+\eta {(R)}_i\right\}.$$

$$\sum _i$$

$$=\eta \left(2\eta -1\right),$$

$$\vec{Q}=\sum _{i=1}^{2\eta }{\vec{q}}_i.$$

$${\vec{Q}}^2$$

$$=2\eta \times \frac{4}{3}{g}^2+2\times \frac{1}{3^{\eta }}\sum _i{E}_i{C}_i$$

$$=0.$$

$$\sum _I{g}^I{N}^I=\eta ,$$

$$\sum _I{\bar{g}}^I{\bar{N}}^I=\eta .$$

$$\sum _I\int d^{3}r{n}_i^I(r)=\eta /3,$$

$$\sum _I\int d^{3}r{\bar{n}}_i^I(r)=\eta /3.$$

$$\sum _i\int d^{3}r{n}_i^I(r)=N^I{g}^I$$

$$\sum _i\int d^{3}r{\bar{n}}_i^I(r)={\bar{N}}^I{\bar{g}}^I$$

$${\hat{n}}_i^I\equiv \frac{3n_i^I}{N^I}$$

$${\hat{\bar{n}}}_i^I\equiv \frac{3{\bar{n}}_i^I}{{\bar{N}}^I}.$$

$$\left(\sum _I\frac{N^I\left(N^I-1\right)}{2}(g^I{)}^2+\sum _{I\ne J}\frac{N^I{N}^J}{2}{g}^I{g}^J+\sum _I{N}^I\frac{g^I(g^I-1)}{2}\right)$$

$$\left(\sum _I\frac{{\bar{N}}^I\left({\bar{N}}^I-1\right)}{2}({\bar{g}}^I{)}^2+\sum _{I\ne J}\frac{{\bar{N}}^I{\bar{N}}^J}{2}{\bar{g}}^I{\bar{g}}^J+\sum _I{\bar{N}}^I\frac{{\bar{g}}^I({\bar{g}}^I-1)}{2}\right)$$

$$\left(\sum _{I,J}{\bar{N}}^I{N}^J{\bar{g}}^I{g}^J\right)$$

$$\sum _I\frac{N^I\left(N^I-1\right)}{2}(g^I{)}^2+$$

$$=\sum _I\frac{N^I{g}^I}{2}\left(N^I{g}^I-1\right)+\sum _{I\ne J}\frac{N^I{N}^J}{2}{g}^I{g}^J$$

Full text

Preprint: BU-HEPP-19-04

††journal: Nucl. Phys. A

1 Introduction

Lattice QCD is a very important tool used by particle physicists to investigate the properties of baryons and mesons. Lattice techniques[1, 2, 3, 4, 5, 6] are presently being employed to understand and elucidate the pentaquark and tetraquark observations made by Belle[7, 8, 9, 10], BESIII[11, 12], LHCb[13, 14, 15, 16, 17], and other collaborations. However, as the quark content increases, it becomes computationally expensive and time-intensive to do the lattice calculations. Every state must be investigated separately, which means a great deal of analysis on Wick contractions and specialized computer coding. In addition, as one adds more quarks, the states will become larger and the lattice used must also increase in volume. There is therefore a need for reliable quark models that can give an overview of many states to help guide these expensive lattice calculations. The MIT bag model[18, 19, 20], the Nambu-Jona-Lasinio[21] model and explicit tetraquark potential models[22] are some of these. Another approach, the Thomas-Fermi (TF) statistical model has been amazingly successful in the explanation of atomic spectra and structure, as well as nuclear applications. Our group has adopted the TF model and applied it to collections of many quarks[23, 24, 25, 26]. One advantage our model has over bag models is the inclusion of nonperturbative Coulombic interactions. One would expect that the TF quark model would become increasingly accurate as the number of constituents is increased, as a statistical treatment is more justified. The main usefulness will be to detect systematic trends as the parameters of the model are varied. It could also be key to identifying families of bound states, rather than individual cases. The TF quark model has already been used to investigate multi-quark states of baryons[24]. In this paper we have extended the TF quark model to mesonic states in order to investigate the stability of families built from some existing mesons and observed new exotic states, concentrating on heavy-light quark combinations. Although our model is nonrelativistic, we will see that this assumption is actually numerically consistent as quark content is increased.

Our paper is organized as follows. In Sec. 2 we will define the energies of the model in terms of the TF density of states. In Sec. 3 we will examine the classical color couplings in mean field theory and systematically determine the probabilities of interactions for quark-quark, quark-antiquark and antiquark-antiquark interactions. We obtain the system energies in Sec. 4, the TF quark equations in Sec. 5, and re-characterize the model in terms of new dimensionless variables in Sec. 6. The application of the model to heavy quark-light quark mesonic matter is initiated in Sec. 7, where we define three types of multi quark-pair families involving charm quarks: charmonium (“Case 1”), $Z$-meson type (“Case 2”) and $D$-meson type (“Case 3”). These states as well as their bottom-quark analogs are constructed in Sec. 7, where we examine the energy slopes to determine family stability. Numerical results and discussions are presented in Secs. 8 and 9, and concluding thoughts and remarks are presented in Sec. 10.

2 The TF meson model

The TF statistical model is a semiclassical quantum mechanical theory developed for many-fermion systems. In this model the assumption is made that quarks are distributed uniformly in each volume element $\Delta V$, while at the same time the quark density $n_{q}(r)$ can vary from one small volume element to the next. For a small volume element $\Delta V$, and for the system of quarks in it’s ground state, we can fill out a spherical momentum space volume $V_{F}$ up to the Fermi momentum $p_{F}$, and thus

[TABLE]

The corresponding phase space volume is,

[TABLE]

Let’s say $g^{I}$ is the degeneracy of a quark flavor $I$. Then quarks in $\Delta V_{{p}{h}}$ are distributed uniformly with $g^{I}$ quarks per $h^{3}$ of this phase space volume, where $h$ is Planck’s constant. The number of quarks in $\Delta V_{{p}{h}}$ is

[TABLE]

The number of quarks in $\Delta V$ is

[TABLE]

where $n_{q}{(\vec{r})}$ is the quark density. Equating Eqs. (3) and (4), we obtain

[TABLE]

The fraction of quarks at position $\vec{r}$ that have momentum between $p$ and $p+dp$ in spherical momentum space is

[TABLE]

Using the classical expression for kinetic energy of a quark with mass $m_{I}$, the kinetic energy per unit volume for a system of quarks is

[TABLE]

where

[TABLE]

Eq. (7) shows that kinetic energy per volume is proportional to the $5/3$-rd power of quark density. To obtain the total kinetic energy ($T$), we will have to integrate this expression over all the spatial volume:

[TABLE]

On the other hand, the total potential energy ($U$) due to the coulombic interactions of system of quarks with one another is given by

[TABLE]

where $V_{N}$ is a factor depending on the type of interaction between quarks. It will be explained in the sections to follow. The total energy for a system of quarks is therefore

[TABLE]

In order to minimize energy while keeping the number of quarks constant, we can add a Lagrange multiplier term of the form $\lambda\left(\int n_{q}(\vec{r})d^{3}(r)-N_{q}\right)$ to the expression for total energy, where $N_{q}$ is the fixed total number of quarks in the given system.

A given quark actually has a definite color and flavor. In the sections to follow we will replace the number density of quarks, $n_{q}(\vec{r})$, by the quantity ${n_{i}^{I}(\vec{r})}$, where color is represented by index $i$ and flavor by index $I$. It will then be summed over flavor and color indices to account for the energies of all of the quarks in the given system.

3 Residual Coulombic coupling and Interaction probabilities

The color couplings of quarks and antiquarks in our model originate from the Coulombic interactions expected at the classical level [27]. In the following we define $\eta$ to be the number of quark/antiquark pairs in the meson, which is assumed to be a color singlet. In addition, $g$ represents the strong coupling constant.

The types of interactions between the particles can then be categorized as:

Color-Color Repulsion (CCR)

Interactions between quarks with same colors is repulsive with coupling constant $4/3g^{2}$. The interactions are red-red ($rr$), green-green ($gg$) and blue-blue ($bb$).

Color-Color Attraction (CCA)

Interactions between quarks with different colors is attractive with coupling constant $-2/3g^{2}$. The interactions are $rb$, $rg$, $bg$, $br$, $gr$, and $gb$.

Color-Anticolor Repulsion (CAR)

Interactions between quarks and antiquarks with different color/anticolors is repulsive with coupling constant $2/3g^{2}$. The interactions are $r\bar{b}$, $r\bar{g}$, $b\bar{g}$, $b\bar{r}$, $g\bar{r}$, and $g\bar{b}$.

Color-Anticolor Attraction (CAA)

Interactions between quarks and antiquarks with same color/anticolors is attractive with coupling constant $-4/3g^{2}$. The interactions are $r\bar{r}$, $b\bar{b}$ and $g\bar{g}$.

Anticolor-Anticolor Repulsion (AAR)

Interactions between antiquarks with same anticolors is repulsive with coupling constant $4/3g^{2}$. The interactions are $\bar{r}\bar{r}$, $\bar{b}\bar{b}$ and $\bar{g}\bar{g}$.

Anticolor-Anticolor Attraction (AAA)

Interactions between antiquarks with different anticolors is attractive with coupling constant $-2/3g^{2}$. The interactions are $\bar{r}\bar{b}$, $\bar{r}\bar{g}$, $\bar{b}\bar{g}$, $\bar{b}\bar{r}$, $\bar{g}\bar{r}$, and $\bar{g}\bar{b}$.

When $\eta$ pairs of quarks interact with each other, the six types of color interactions appear with different probabilities. So, we need to find a way to count and formulate color-interaction probabilities in order to apply the TF model for mesonic matter.

Each of the $\eta$ pairs must be a color singlet and each singlet can be achieved three different ways: red-antired, green-antigreen and blue-antiblue. We will term each quark-antiquark pair a “pocket”. Out of $\eta$ pockets let us say we have $x$ number of red-antired, $y$ number of blue-antiblue and $z=\eta-x-y$ number of green-antigreen combinations, which will be represented by $R$, $B$ and $G$, respectively. This is depicted in Table 1.

When two $R$ pockets interact for example, there are four color combinations as depicted in the left box of Table 2. Out of these, $r\bar{r}$ and $\bar{r}r$ belong to interaction type CAA, $rr$ to CCR and $\bar{r}\bar{r}$ to AAR. The middle box $RB$ depicts color combinations when different types of pocket interact. The rightmost box $R$ represents a color combination within this pocket, which means there is interaction type CAA only. These color combinations are counted and depicted in Table 3. Of course $BB$ and $GG$ have the same interaction numbers as in the left box. In addition, $BG$ and $GR$ have same interaction numbers as in the middle box.

In Table 3 the six interaction types are given an index. The last three columns give the number of interactions contained in the previous table.

There are $x(x-1)/2$ number of ways the $RR$ type of interaction takes place. Similarly, ${y(y-1)}/{2}$ and ${z(z-1)}/{2}$ are the number of times $BB$ and $GG$ happens. Also, $RB$, $BG$ and $RG$ type of interactions occur $xy$, $yz$ and $xz$ times respectively.

In the next step we varied $x$ from 0 to $\eta$ and $y$ from 0 to $\eta-x$ thereby giving equal footing to all the color combinations and counted all the possible color interactions, $E_{i}$. The corresponding expression is given in the equation below.

[TABLE]

Here $E_{i}$ refers to total occurrence of the $i^{th}$ color coupling. For example, $E_{1}$ refers to total number of times the AAA type of interaction occurs out of $3^{\eta}\times\eta(2\eta-1)$ possibilities.

To understand this equation, let us break it into two parts. First, if we take the terms inside the curly braces and add over the indices, we get $\eta(2\eta-1)$, i.e.,

[TABLE]

regardless of $x$, $y$ and $z$ values as long as the constraint $x+y+z=\eta$ is satisfied. This is just the number of ways $2\eta$ quarks can be paired. Second, if we replace the curly braces by unity, we get $3^{\eta}$, the total number of pocket combinations in a meson with $\eta$ pairs of quarks. Thus, $E_{i}$ gives the total occurrence of the $i^{th}$ color interaction out of $3^{\eta}\eta(2\eta-1)$ possibilities.

After solving Eq. (12) using Mathematica, we obtained the outcomes, $E_{i}$. They are listed in Table 4. In the same table, we divided each event by $3^{\eta}\eta(2\eta-1)$ and obtained probabilities for each type of interaction. In addition, the coupling strength, $C_{i}$, is defined for each interaction type.

If we add the product of coupling and probabilities from Table 4, we get $-{\frac{4}{3}g^{2}}/{(2\eta-1)}$, very similar to the baryon case[24]. The negative sign indicates that the system is attractive because of the collective residual color coupling alone, even in the absence of volume pressure. This gives rise to a type of matter that is bound, but does not correspond to confined mesonic matter, as discussed in Ref. [23]. We are interested in confined matter and will need to add a pressure term to the energy to enforce this; see Eq. (53) below.

As a check on our counting, consider the total color charge of the quarks:

[TABLE]

Squaring on both sides, we get

[TABLE]

This shows our model is an overall color singlet. Note that we divided by $3^{\eta}$ in the second term to average over all the possible configurations.

The TF quark model replaces the sum over particle number in particle interaction models with an integral over the density of state particle properties. We will weight interaction strengths, taken from the classical theory, by the probabilities of various interactions in the color sector, which we assume to be flavor independent. So, we need a connection between particle number and probability. We will assume as in other TF quark models that these interaction probabilities are proportional to the number of particle interaction terms. Also, treating color-combinations on an equal footing, each color now shares one-third of the total probability. Therefore, we divide the probabilities in Table 4 by three as shown in Table 5. Also, notice that we have generated a new way of representing of probabilities using double color indices $i$ and $j$ and barred and double barred symbols for $P$. No bar is for color-color probabilities, one bar is for color-anticolor probabilities, and two bars is for anticolor-anticolor probabilities. From Table 5, we can see that we correctly obtain $\sum_{i\leq j}^{3}{P_{ij}+\bar{P}_{ij}}+\bar{\bar{P}}_{ij}=1$ for the sum of all the probabilities.

4 System energies and equations

Our multi quark-pair system consists of an equal number of $\eta$ quarks and $\eta$ anti-quarks. Each quark or antiquark carries both color and flavor. So, let us introduce $N^{I}$ as the number of quarks with flavor index $I$, and $\bar{N}^{I}$ as number of anti-quarks with anti-flavor index $I$ such that

[TABLE]

and

[TABLE]

In terms of quark density these equations can be expressed as

[TABLE]

and

[TABLE]

In Eqs.(18) and (19) degeneracy factors are already included in quark densities $n^{I}_{i}(r)$ and ${\bar{n}}^{I}_{i}(r)$.

Similarly, the sum over the color index gives the total number for a given flavor. Thus

[TABLE]

and

[TABLE]

Also, for the convenience, we will introduce the single-particle normalized density

[TABLE]

and

[TABLE]

This form of quark density will be helpful in correctly normalizing the TF interaction energy when continuum sources are used.

In the earlier section, we calculated probability for six types of interactions between colors. We will now see how these color interactions are associated with flavor numbers.

For a quark-antiquark system with $2\eta$ total number of particles, the number of interactions possible is $\eta(2\eta-1)$. Out of these CCA and CCR types occur

[TABLE]

times, AAR and AAA occurs

[TABLE]

times, and CAR and CAA occurs

[TABLE]

times. For a consistency check, we are going to add all of these seven terms. We begin with the three CCA and CCR terms, for which we have

[TABLE]

The interaction between antiflavors, i.e, the AAR and AAA terms, also yields $\eta^{2}/2-\eta/2$. The last term involving one flavor and one anti-flavor gives $\eta^{2}$. Adding all of them yields $\eta(2\eta-1)$ as it should.

We will now use this information to write expressions for kinetic and potential energies. Building up the expression from Section 2, we can write

[TABLE]

for the kinetic energy from quarks and antiquarks. Here we have assumed that radius of the objects are finite. Similarly, assigning the probabilities and couplings from the interaction terms shown in Table 5, we can now define total potential energy to be

[TABLE]

In Eq. (26) $\hat{n}_{i}^{I}$ are number densities, which are normalized to one when integrated over space[23]. Note that the degeneracy factors are already contained in the expression for number densities. Therefore, we divided the third and sixth terms by the square of degeneracy factors.

From Table 5, we can see that the probability of interaction type CCA is twice the CCR type. Also, interaction probability AAA is twice the AAR type. This simple finding amazingly removes six out of seven terms from the interaction energy, leaving us with the last term from (26). Thus we conclude, for mesons on average, quarks only interact with antiquarks. This cancellation makes our analytical solutions easier to achieve.

The total energy $E$ can now be written as the sum of kinetic and potential energies such that

[TABLE]

Now, we switch back to normalization $n^{I}_{i}$ and $\bar{n}^{{I}}_{i}$ and assume equal Fermi color momenta, $n^{I}\equiv n^{I}_{1}\equiv n^{I}_{2}\equiv n^{I}_{3}$ for each I. The same assumption is made for antiparticles. This gives

[TABLE]

Eqs. (20) and (21) can now be averaged over colors as

[TABLE]

and

[TABLE]

We will use this pair of equations to set up the normalization conditions.

5 Thomas Fermi Quark Equations

We can now introduce Lagrange’s undetermined multipliers $\lambda^{I}$ and $\bar{\lambda}^{{I}}$ associated with constraints in Eqs.(29) and (30) and add them to the expression for energy. This gives

[TABLE]

The purpose of adding these terms involving Lagrange multipliers is to allow a minimization of the total energy while keeping particle number constant. The density variations $\delta\bar{n}^{I}(r)$ and $\delta n^{{I}}(r)$ give

[TABLE]

[TABLE]

Assuming spherical symmetry, the TF spatial functions $f^{I}(r)$ and $\bar{f}^{I}(r)$ are defined such that

[TABLE]

[TABLE]

where

[TABLE]

gives the scale, where $m^{1}$ is the mass of lightest quark, and

[TABLE]

is the strong coupling constant. Equation (32) after using Eqs.(34) - (37) can now be written as

[TABLE]

Here we have used the integral

[TABLE]

resulting from spherical symmetry. Equation (38) can be further simplified into

[TABLE]

At this stage we introduce the mass ratios

[TABLE]

for both the particle and antiparticle cases.

Let us introduce the dimensionless parameter $x$ such that $r=Rx$ where

[TABLE]

Eq. (40) can now be written as

[TABLE]

Similarly, starting from Eq. (33) we obtain the other TF integral equation:

[TABLE]

Taking first derivatives of Eqs. (43) and (44), we have

[TABLE]

and

[TABLE]

Similarly, second derivatives of Eqs. (43) and (44) yield

[TABLE]

and

[TABLE]

Equations (47) and (48) are the differential form of the TF quark equations in the case of mesons. The interchangeability of these equations shows the TF equations are invariant with respect to particle and antiparticle. As was mentioned before, it also shows that quarks interact only with antiquarks in mesonic matter; quark/quark and antiquark/antiquark interactions sum to zero in the TF model. When there is an explicit particle/antiparticle symmetry, we assume $f=\bar{f}$ to reduce the TF differential equations into one incredibly simple form:

[TABLE]

We use this equation to form system energies and equations for charmonium/bottomonium families and $Z$-meson family systems. However, for $D$ and $B$-meson families, the explicit asymmetry in the masses of the particle and antiparticle will require us take a different approach.

6 Energy equations in terms of dimensionless radius

The expressions for kinetic and potential energies now need to be given in terms of the dimensionless distance in order for us to be able to obtain analytical solutions. In addition, the volume energy term that produces an external pressure on the system needs to be appropriately introduced.

To obtain the kinetic and potential energy expressions, we begin from Eq. (28) and apply Eqs. (34), (35) and (42). This gives the kinetic energy as

[TABLE]

The potential energy

[TABLE]

becomes

[TABLE]

The volume energy ($E_{v}$) term gives the inward pressure from the vacuum. We assume that[18]

[TABLE]

where $B$ is the bag constant. Now that we have $T$, $U$ and $E_{v}$ terms, we can find the total energy of a desired multi quark-pair state. The total energy of such a state is simply given by

[TABLE]

where $m_{q}$ and $\bar{m}_{q}$ are the mass of the quark and anti-quark respectively.

In Section 8 we will fit model parameters for a given set of mesons, including $m_{q}$ and $\bar{m}_{q}$. However, to assess the stability of multi quark-pair mesons we omit the trivial mass part and will examine the energy as a function of quark content.

7 Application of model

In this section we will now apply our model and obtain energy expressions for different families of quarks. In all the cases that follow charm can be read as bottom also.

7.1 Charmonium family: Case 1

This type of family contains quarks and anti-quarks of equal mass. The system of quarks in this family can be represented as $Q\bar{Q},Q\bar{Q}Q\bar{Q},Q\bar{Q}Q\bar{Q}Q\bar{Q}$ and so on. We have investigated the case where $Q$ is a charm quark or bottom quark but not both in the same system. We call it the charmonium family. Charmonium and all the multi quark-pair families of charmonium consists of charm and anti-charm (or bottom and anti-bottom) only. Therefore, $I$ takes a single value and hence is dropped. $g_{0}$ refers to degeneracy and, for this application, can have the value of one or two to represent spin. This can give spin splittings, but we have not yet included explicit spin-spin interactions in the model, like we have done in our baryon model[24]. Eq. (49) becomes

[TABLE]

where

[TABLE]

We choose the normalization equation to be

[TABLE]

Expressing the normalization equation in terms of boundary conditions, we have

[TABLE]

With these modified boundary conditions, we can derive the expressions for kinetic and potential energies. For the kinetic energy we can start with Eq. (50) and assume a single flavor. In this case, both quarks and antiquarks contribute equally to the kinetic energy, which gives

[TABLE]

The integral may be done using the TF differential equation, and results in

[TABLE]

Thus, the kinetic energy depends only on the derivative and value of the TF function at the boundary. When there is a single flavor as in the case of the charmonium family, and symmetry of particle and antiparticle occurs, the potential energy Eq. (52) can be written as

[TABLE]

This can be further simplified to

[TABLE]

which like the expression for $T$ depends on the derivative and values of the TF function at the boundary. It is important to realize that the “a” value in this section differs from others where the lightest quark is not charm or bottom.

7.2 $Z$-meson family: Case 2

The constituents of $Z$-mesons are charm ($c$) (or bottom($b$)), anti-charm ($\bar{c}$) (or anti-bottom ($\bar{b}$)), light ($u$ or $d$) and anti-light quarks. The particles in this family can be represented as $Q\bar{Q}q\bar{q}$, $Q\bar{Q}q\bar{q}Q\bar{Q}q\bar{q}$, $Q\bar{Q}q\bar{q}Q\bar{Q}q\bar{q}Q\bar{Q}q\bar{q}$, and so on, where $Q$ represents a heavy quark and $q$ is a light quark. We treat the mass of the up and down quarks as the same. This means the $Z$-meson and all multi quark-pair families of $Z$-mesons have a total quark mass equal to the antiquark mass. As before, we set $f=\bar{f}$ in the TF equations and obtain Eq. (49) with $I=1$ or 2. Let $f^{1}(x)$ be the TF function of the light quark and $f^{2}(x)$ be the TF function of the heavier quark. For $N_{1}$ quarks with degeneracy factor $g_{1}$, and $N_{2}$ quarks with degeneracy $g_{2}$, we have

[TABLE]

In our application $g_{1}$ can have values 1, 2, 3, or 4 whereas $g_{2}$ can have value of 1 or 2 only. We assume a linear relation exists between $f^{1}(x)$ and $f^{2}(x)$ in the region $0<x<x_{2}$, and that $f^{2}(x)$ vanishes for a dimensionless distance greater than $x_{2}$, i.e.,

[TABLE]

Eq. (49) can now be written as two set of equations:

[TABLE]

[TABLE]

To make Eqs. (65) and (66) consistent, we need

[TABLE]

which is just the inverse ratio of the given masses from (41). The similar step in the case of baryons gives a much more complicated consistency condition[24]. The normalization conditions are

[TABLE]

and

[TABLE]

In region $0<x<x_{2}$,

[TABLE]

where

[TABLE]

In region $x_{2}<x<x_{1}$,

[TABLE]

where

[TABLE]

With the equations above we can express normalization conditions in the form of boundary conditions:

[TABLE]

[TABLE]

The energies are derived as for Case 1. For $g_{1}$ flavors with $N_{1}$ particles and $g_{2}$ flavors with $N_{2}$ particles we have

[TABLE]

Using the TF function differential equations, the consistency condition for $k$ and boundary conditions allows one to relate the integrals to TF function values and derivatives on the surfaces:

[TABLE]

The expression for potential energy in the same case becomes

[TABLE]

which can be placed into the form

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

One can show that the $K_{12}$ integral is equivalent to the $K_{21}$ integral.

At this point we would like to take some time to explain the reductions of these expressions. We will see some interesting cancellations in the analytical solutions. The expressions for $K_{1}$, $K_{2}$ and $K_{12}$ after integrations and reductions can be written as

[TABLE]

[TABLE]

and

[TABLE]

Similarly, the potential energy now simplifies to

[TABLE]

which on further reduction yields

[TABLE]

Here we observe that the coefficient of $\frac{N_{2}}{\eta}\frac{df^{2}(x)}{dx}|_{x_{2}}$ vanishes; all the above eleven terms actually cancel. Also all the other seemingly difficult integrals boil down to a simple equation, which gives

[TABLE]

So, the interesting thing we observe in both the kinetic and potential energy expressions is that there is no dependency on the derivative of TF function at the inner boundary, unlike the baryon case.

7.3 $D$-meson family: Case 3

The $D$-meson and multi quark-pair families of $D$-mesons consist of charm (or bottom) quarks (or antiquarks) and light antiquarks (or light quarks). All the quarks in this family can be represented as $Q\bar{q}$, $Q\bar{q}Q\bar{q}$, $Q\bar{q}Q\bar{q}Q\bar{q}$, and so on. We will study the systems where the heavier mass is a quark and the lighter mass is an antiquark in each member of the $D$-meson family. By the symmetry inherent in the TF equations, this also covers the antiparticle state where the particle and antiparticle are interchanged. There is asymmetry in total mass of quarks and antiquarks, unlike the other cases we have studied. Also, because quarks interact only with antiquarks, the mathematics is different than earlier cases. We have two different TF functions for the charm and anti-light quark, $f$ and $\bar{f}$, respectively. We assume there is a universal TF function $\bar{f}=k_{0}f$ in the region where the TF functions overlap. This means they are related linearly in the smaller region where the TF function of the charmed quark is nonzero. Outside of this, only the TF function of the light antiquark exists.

For the region $0<x<x_{2}$, we have the differential equation

[TABLE]

where

[TABLE]

and

[TABLE]

is the consistency condition. $\alpha$ is given by Eq. (41), and in this case $\bar{\alpha}=1$.

In region $x_{2}<x<x_{1}$ with two TF equations, one function is zero and the other becomes

[TABLE]

where

[TABLE]

and

[TABLE]

Proceeding as before we obtained expressions for kinetic and potential energies. The energies are

[TABLE]

[TABLE]

As a consistency check, if we assume light and heavy quarks to have equal mass, we have found that the potential and kinetic energies with $\eta=1$ in Case 3 are equal to that of Case 1. Also, the non-degenerate Case 2 with $\eta=2$ is the same as degenerate Case 1 with $\eta=2$. We confirmed this both in our analytical and numerical results.

8 Method and Remarks

The phenomenological parameters we need for our model are the strong coupling constant $\alpha_{s}$, the bag constant $B$, the charm and bottom quark masses, $m_{c}$ and $m_{b}$, as well as the light quark mass, $m_{1}$. Previously, we used baryon phenomenology to obtain these parameters[25]. In addition, we did a fitting using mesonic states that involves only the charm quarks[26]. Here, we have included bottom quarks as well and obtained the fitted parameters. Since we do not yet include spin interactions in our model, we need to weight spin-split states to “remove”this interaction for our model fits.

Thus, we first fit the model expressions for the masses of charmonium states, bottomonium states, spin-weighted $D$-meson states, spin-weighted $B$-meson states, a spin-weighted combination of spin 1 likely tetraquark states involving charm and bottom quarks. In this way, we obtained the spin-weighted mass for these six states and then fit the parameters above to obtain the minimized chi-square.

For the mass of charmonium and bottomonium we weighted masses of $1S$ states such that

[TABLE]

[TABLE]

which we will refer to as Case 1-charm and Case 1-bottom mass.

To obtain Case 2-charm, we weighted spin 1 likely tetraquark states called the $Z_{c}(3900)$ and $X(4020)$. There are a number of charmonium-like exotic resonances which have been discovered in recent years[28]. The ground state of this set is the $\chi_{c1}(3872)$, with a mass of 3871.7 MeV, discovered by the Belle Collaboration in 2003[7]. However, this state appears to be a isospin singlet, although the charge states might not have been seen yet[22]. On the other hand, the $Z_{c}(3900)$ resonance, with mass 3886.6 MeV, spin 1 and $C=-1$, is a triplet as we would expect from the hidden charm nonrelativistic tetraquark model. In addition, the fact that the $Z_{c}(3900)$ is above the $D^{0}\bar{D}^{*+}$ threshold makes it more likely that it is a true tetraquark rather than a molecular state of these same two particles. In the same way, the $X(4020)$ with mass 4024.1 MeV is a spin 1, $C=-1$ state which in this case is just above $D^{*0}\bar{D}^{*+}$ threshold. We will adopt these two particles as our spin-split charmed tetraquark states. Note that there are hidden bottom analogs of the $Z_{c}(3900)$ and $X(4020)$ states in the $Z_{b}(10610)$ and $Z_{b}(10650)$ states. We will use these to define the Case 2 $Z$-mesons.

Note that a very simple model which can accommodate the $\chi_{c1}(3872)$, $Z_{c}(3900)$ and $X(4020)$ states has a spin interaction Hamiltonian given by

[TABLE]

where $s_{1}$ is the light quark spin, $s_{2}$ is the heavy quark spin and $J$ is the total spin. The three terms represent the light-light, heavy-heavy and light-heavy spin interactions, respectively. We would expect that $\kappa_{2}<\kappa_{1},\kappa_{3}$ based on the quark masses. By neglecting $\kappa_{2}$ and fitting the three states mentioned, one predicts that there are additional hidden charm spin 0 states at $3720$ and $3887$ MeV, as well as a spin 2 state at $4176$ MeV, all with $C=1$. There is a very broad spin 0 candidate state called the $\chi_{c0}$ at $3862$ MeV as well as a narrower spin 0/2 sate called the $X(3915)$ at 3918 MeV, which are both apparently isospin singlets, but the other predicted states do not presently have candidates[28].

The Case 2 calculations are now simply,

[TABLE]

and

[TABLE]

which we will refer to as Case 2-charm and Case 2-bottom mass respectively.

For the mass of $D$ mesons we weighted the spin 0 charge states such that

[TABLE]

and similarly for spin 1 $D^{*}$ mesons. We thus obtained the Case 3-charm mass

[TABLE]

Similarly for B mesons, we obtained the Case 3-bottom value as 5313 MeV.

Using Mathematica, we fitted parameters such that the mass chisquare was minimized using a grid search, obtaining a total $\sqrt{\chi^{2}}=86.1$ MeV spread over 6 masses. We solved the differential equations using an iterative implementation of NDSolve in Mathematica. We obtained $\alpha_{s}=0.346$, $B^{1/4}=107.6$ MeV, charm quark mass $m_{c}=1553$ MeV, and bottom quark mass $m_{b}=4862$ MeV. Our light quark mass, $m_{1}=306$ MeV, we take from our previous TF baryon fit[24]. Our baryon paper found $\alpha_{s}=0.371$ and $B^{1/4}=74.5$ MeV. Note that it would be premature to compare these sets of parameters as the baryon fit included spin interactions, which are quite significant, and the present fit does not. We will have additional comments on the inclusion of spin for TF mesons at the end of Section 9. Table 6 shows the difference between expected and obtained masses.

We will examine TF functions and energies for the three cases defined above. In nuclear physics, one examines the binding energy per nucleon in order to assess the stability of a given nucleus. We will do a similar investigation here. Thus, the important figure of merit in these evaluations is the total energy per quark, for if this increases as one adds more quarks, the family is unstable under decay to lower family members, whereas if it decreases, the family is stable. The static mass quark dependence of the the mesons does not play a role in these considerations and so will be left off.

9 Results and Discussions

First, let us first discuss the behavior of the TF density functions. The particle density function is proportional to $(f(x)/x)^{3/2}$ for all states. The charmonium and bottomonium function decreases smoothly in distance and has a discontinuity at the boundary as seen in Fig. 1 and Fig. 2, where the dimensionless $x$ variable is used. The discontinuity for the bottomonium density function actually is almost zero at the boundary.

The density function of $Z$-mesons has a long tail for the light quarks, while for charmed quarks the value is large and is concentrated near the origin, as seen in Fig. 3 and Fig. 4. This suggests an atomic-like structure with heavy charm, anti-charm quarks at the center while light quarks and antiquarks spread out like electrons. Figs. 5 and Fig. 6 are an enlargement of the density function of the light quark TF function for the $Z$-meson. It drops down abruptly until it reaches the boundary of the heavy quark TF function, then inflects and decreases. In the case of D and B-mesons, Fig. 7 and Fig. 8 respectively, the density function of light and heavy quarks are relatively closer. We can see that when charm is interchanged with bottom quarks the density function is more concentrated near the center, as one would expect.

We increased the quark content and compared density functions of a family of multi-mesons in all three cases. We observed similar density functions for a given multi-meson family regardless of the quark content. This is the embodiment of the Thomas-Fermi approach.

The Isgur-Wise symmetry methodology[29] in $D$ and $B$ mesons is a capable tool in particle physics and has been implemented in various studies[30, 31, 32, 33]. It arises when the heavy quark appears to the light quark degrees of freedom as a static color source[34]. There is a natural question of how close our simulations for multiquarks approach such a heavy-light symmetry scenario. As one can see by an examination of the TF functions in Figs. 4 and 8, there is never a case in which one can consider the radial extent of the heavy quark density function exceptionally compact, and thus kinematically nondynamical, compared to the light quark. Although Figs. 7 and 8 refer to $D$ and $B$ systems, remember that the TF approach actually embodies multiquarks, as discussed above. Also, note that the Case 2 and Case 3 heavy quark profiles are very different, reflecting the different interactions present among the quarks and antiquarks. The Case 2 heavy quarks interact with both heavy and light antiquarks, whereas the Case 3 heavy quarks interact only with the light antiquarks. Since the interactions are attractive this makes the Case 2 heavy quark profile more compact than Case 3. For light atomic elements one can neglect electron-electron interactions to lowest order, which is actually more like Case 3. Thus, even though the Case 2 heavy quark profile is more compact like an atomic nucleus, the interactions make the Case 2 mesons less atomic-like than Case 3. These TF functions and interactions suggest it is unlikely that a nondynamical heavy flavor symmetry will be manifest in heavy-light multiquark states for charm and bottom quarks.

Fig. 9 and Fig. 10 includes 11 possibilities for the physical radius. The physical radius is plotted versus quark number and compared with a generic baryon with three degenerate light flavors. From Eq. (42) we can see that the radius is proportional to the product of $\eta^{2/3}$ and the dimensionless radius $x$. As the quark content increases, the dimensionless radius, $x$, becomes smaller, whereas $\eta^{2/3}$ increases. In most cases, the result is an increase in radius with increasing quark content. We also observe that the curve of the radius plot for each case tends to flatten out for larger numbers of quarks. case1nf1xmax refers to multi quark-pair families of charmonium with no degeneracy. In this case all the charm quarks have the same spin and hence cannot occupy the same state. case1nf2xmax is instead the plot of the charmonium family with a degeneracy of two. In this case, spin up and down are assigned to a pair of charm quarks, and thus they do not occupy the same state. Thus, the physical radius of case1nf1xmax is larger than case1nf2xmax reflecting this fact. Note that the dotted lines refer to the inner boundary associated with the charmed quark in Cases 2 and 3. For Case 2 the difference between dotted and continuous lines is the largest. That means, like the $Z$-meson, all the higher quark family members have heavy charm-anticharm concentrated at the center while light and anti-light quarks are spread throughout. case2nf1x2 refers to the radius plot of the inner boundary of the multi quark-pair family of $Z$-mesons with degeneracy of one while case2nf2x2 refers to the same plot with degeneracy of two, and similarly for Case 3. In all cases we see that the plot of physical radius is larger for degeneracy one compared to degeneracy two. The $Z$ and $D$-meson family members are found to have equally large outer boundaries. For a given quark number, the outer radius of all types of mesons was found to be smaller than the generic baryon. Note that these radii are determined by the size of the TF functions; the electromagnetic radii have not yet been evaluated.

There are 3 types of energies in this model: kinetic, potential and volume. The kinetic energy per quark (in MeV) depends strongly on the meson family, as seen in Fig. 11 and Fig. 12. We see that the energy per quark is relatively small and tends to decrease slowly, which seems to provide some justification for this nonrelativistic model. Fig. 13 and Fig. 14 shows the corresponding graph for the potential energy. The energies rise quickly from rather large negative values to saturate near small negative values for all cases. Fig. 15 and Fig. 16 shows the graph for volume energies. The graphs are similar when we interchange bottom with charm quarks. The case2nf2 and case3nf2 results in these figures deserve some extra comments. These lines are determined by the most degenerate state, and thus the smallest energy per quark, available for the given quark content. If we denote the total number of quarks and antiquarks as $N$ ($=2\eta$), the $N=$ 8, 16 and 24 quark cases for case2nf2, which all have $g_{2}=2$, are associated with $g_{1}=2$, $g_{1}=4$ and $g_{1}=3$, respectively. Likewise, the $N=$ 4, 8, 12, 16 and 20 quark cases for case3nf2 are associated with $g_{1}=2$, $g_{1}=4$, $g_{1}=3$, $g_{1}=4$ and $g_{1}=2$, respectively.

Figs. 17 and 18 show the total energy per quark without the mass term, i.e., the sum of kinetic, potential and volume energy, plotted against the quark content. The generic baryon rises slowly and smoothly for increasing quark content, implying these are unstable; i.e., a higher quark content state can decay into lower members of the same family. The Case 1 mesons rise quickly, and then continue the rise more slowly; these are also unstable. The Case 2 and 3 nondegenerate mesons rise less quickly but still continue a slow increase for higher quark numbers.

The lowest mass Case 2 and 3 patterns have been highlighted in Figs. 19 and 20. For Case 2 we can observe negative slope from $N=4$ to $N=8$ and again from $N=8$ to $N=16$ when we account for degeneracy. Here we have plotted five points for Case 2 with both charm-light quarks and bottom-light quarks. The charm and bottom curves are very close to one another. After that there is an inflection at $N=24$. This hints at the possible existence of stable multi-mesons with $N=8$ and also with $N=16$, which we term octaquarks and hexadecaquarks. On the other hand, the Case 3 charm and bottom patterns are separated from one another, but are still very similar. There is an overall upward trend to this data, although there is an isolated downward dip from $N=12$ to $N=16$. We will comment on the $N=4$ results from this figure below.

There are lattice QCD results for Case 3 tetraquarks. Ref. 1 has two-point function results for quark flavor content $ud\bar{c}\bar{c}$ with spin 1 and $uu\bar{c}\bar{c}$ content for spin 0 for dynamical pion masses ranging from $257-688$ MeV. Their continuum extrapolations show binding of $-23.3\pm 11.4$ MeV for the spin 1 sector and antibinding (resonance energy) of $25.9\pm 10.9$ MeV above threshold in the spin 0 sector. On the other hand, the $ud\bar{b}\bar{b}$ with spin 1 has a reported binding of $-143.3\pm 33.9$ MeV and $uu\bar{b}\bar{b}$ has a small spin 0 binding of $-5.5\pm 17.7$ MeV. We can not compare directly with these spin sector results because we have not yet constructed the spin states as in our baryon study[24]. The best we can do at this stage is to point out that our Case 3 results from Fig. 20 double charm quarks show a very small downward jump in energy per quark from $N=2$ to $N=4$ of $-0.67$ MeV, corresponding to a total binding energy of $-2.68$ MeV, whereas our double bottom results from the same figure show a small upward jump of $1.98$ MeV, leading to total antibinding of $7.9$ MeV. Our lowest energy states are always the most degenerate, which would correspond to the spin 0 tetraquark state in the doubly heavy sector. Although it is a bit premature to compare, the reported Monte Carlo error bars in this reference can accommodate our results at the $\sim 2\sigma$ (double charm case) or $\sim 1\sigma$ (double bottom case) level. More precise studies will be necessary to make definitive conclusions concerning the binding state of these types of particles.

In our future work, we will need to include the spin interactions. In mesonic models, the color probabilities only allow quark-antiquark interactions, which greatly simplifies the spin calculations. A derivation of the color magnetic interaction similar to that in Ref. [24] gives

[TABLE]

The sums are over generalized flavor indices $\cal I,J$ (which include spin) as well as particle indices $i,j$. The $(S_{z}^{\cal I})_{i},(\bar{S}_{z}^{\cal J})_{j}$ are the spin components and the $\gamma_{\cal I},\bar{\gamma}_{\cal J}$ factors in 105 above are given by

[TABLE]

where $m_{\cal I},\bar{m}_{\cal J}$ is the quark or antiquark mass and g is the color gyromagnetic factor. Note that the energy shift in 105 is for each TF configuration; as in [24] the particle states are weighted sums of the various TF configurations. As pointed out above, including the spin interaction will affect the other model parameters and energy extrapolations. However, comparison of the previous baryon, Ref. [23], with newly fit meson parameters will allow us to evaluate systematic uncertainties within the TF quark model approach in future spectroscopic work.

10 Conclusions and Acknowledgements

We have initiated the study of multi quark-pair mesons using the TF quark model. After specifying the explicit interactions and summing on colors, we have formulated system interactions and energies in a mean field approximation. We have investigated three cases of mesonic states: charmonium family, $Z$-meson family and $D$-meson family, as well as their bottom quark analogs. We have not yet included explicit spin interactions in our model, but we can take one level of degeneracy into account in our two-TF function construction.

We have observed interesting patterns of single quark energies. Similar to our findings for baryons, the energy per quark is slowly rising for the nondegenerate Case 1 and 3 mesons, implying family instability. Our Case 3 finding for tetraquarks and other most degenerate states can not be considered definitive because of the lack of spin splitting terms in our interaction, but the overall trend argues against a family of such higher quark number states. Our Case 2 findings are the most interesting of our study. It is the only case where we see an actual sustained decrease in the energy of introduced quark pairs. This would indicate that stable octaquark and hexadecaquark versions of the charmed and bottom $Z$-meson exist. However, the energy trend is not smooth and the spin interactions have also not been taken into account. These findings are complicated dynamical results and could not have been predicted from first principles.

Our first order of business as we extend the model will be the inclusion of explicit spin interactions, to bring our meson model to the same level of development as the baryon model. These interactions can be determined from the nonrelativistic ground state wave functions of these states and the use of Eq. (105). A further extension of this model would be to examine mixed baryonic-mesonic states such as pentaquark families.

We thank the Baylor Quantum Optics Initiative and the University Research Committee of Baylor University for their partial support of this project. We also gratefully acknowledge discussions with N. Mathur as well as helpful considerations from G. Chandra Kaphle.