Application of the Thomas Fermi Quark Model to Multiquark Mesons
Suman Baral Walter Wilcox
hep-ph hep-lat








mymainaddress,mysecondadd,mythirdadd]Suman Baral
[mymainaddress]Department of Physics, Baylor University, Waco, TX USA 76798-7316
[mysecondadd]Everest Institute of Science and Technology, 343 Ranibari marga, Kathmandu, Nepal
[mythirdadd]Neural Innovations LLC, 3225 Skinner Dr., Lorena, TX USA 76655
suman@everestsciencetech.com; suman@neuralinnovations.io



mymainaddress]Walter Wilcoxmycorrespondingauthor
walter_wilcox@baylor.edu
[mycorrespondingauthor]Corresponding author




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. 



quark model Thomas-Fermi tetraquarks quark matter
12.39.-x12.39.Mk12.40.Ee 21.65.Qr

Preprint: BU-HEPP-19-04








§ INTRODUCTION


Lattice QCD is a very important tool used by particle physicists to investigate the properties of baryons and mesons. Lattice techniques<cit.> are presently being employed to understand and elucidate the pentaquark and tetraquark observations made by Belle<cit.>, BESIII<cit.>, LHCb<cit.>, 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<cit.>, the Nambu-Jona-Lasinio<cit.> model and explicit tetraquark potential models<cit.> 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<cit.>. 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<cit.>. 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. <ref> we will define the energies of the model in terms of the TF density of states. In Sec. <ref> 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. <ref>, the TF quark equations in Sec. <ref>, and re-characterize the model in terms of new dimensionless variables in Sec. <ref>. The application of the model to heavy quark-light quark mesonic matter is initiated in Sec. <ref>, 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. <ref>, where we examine the energy slopes to determine family stability. Numerical results and discussions are presented in Secs. <ref> and <ref>, and concluding thoughts and remarks are presented in Sec. <ref>. 




§ 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 Δ 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 Δ 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

    V_F≡4/3π p_F^3( r⃗).

The corresponding phase space volume is,

    Δ V_ph≡ V_FΔ V = 4/3π p_F^3(r⃗)Δ V.

Let's say g^I is the degeneracy of a quark flavor I. Then quarks in Δ V_ph 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 Δ V_ph is

    Δ N ≡g^I/h^3Δ V_ph  = 4π g^I/3h^3 p_F^3(r⃗)Δ V.

The number of quarks in Δ V is

    Δ N ≡ n_q(r⃗)Δ V,

where
 n_q(r⃗) is the quark density.
Equating Eqs. (<ref>) and (<ref>), we obtain

    n_q(r⃗)=4π g^I/3h^3p_F^3(r⃗).

The fraction of quarks at position r⃗ that have momentum between p and p+dp in spherical momentum space is
 
    F_(r⃗)(p)dp =
    4π p^2 dp/4/3π p_F^3(r⃗),         if p ≤ p_F(r⃗), 
    
          0,                   otherwise.

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

    t(r⃗)  = ∫^p_Fp^2/2m_In_q(r⃗) F_(r⃗)(p)dp  = C_F ( n_q(r⃗))^5/3,

where

    C_F = ( 6π^2 ħ^3 )^5/3/20 π^2 ħ^3 m_I( g^I)^2/3.

Eq. (<ref>) 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:
 
    T = ∫  t(r⃗)  d^3r.

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

    U= V_N ∫∫ n_q(r⃗) n_q(r⃗ ') /|r⃗-r⃗ '| d^3r  d^3r',

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

    E    = T+U 
       
    = C_F ∫(n_q(r⃗) )^5/3 d^3r + V_N ∫∫ n_q(r⃗) n_q(r⃗ ') /|r⃗-r⃗ '| d^3r d^3r'.

In order to minimize energy while keeping the number of quarks constant, we can add a Lagrange multiplier term of the form λ( ∫ n_q(r⃗) d^3(r) - N_q) 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(r⃗), by the quantity n_i^I(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.



§ 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 <cit.>. In the following we define η 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 rb̅, rg̅, bg̅, br̅, gr̅, and gb̅.


Color-Anticolor Attraction (CAA)

Interactions between quarks and antiquarks with same color/anticolors is attractive with coupling constant -4/3g^2. The interactions are rr̅, bb̅ and gg̅.


Anticolor-Anticolor Repulsion (AAR)

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


Anticolor-Anticolor Attraction (AAA)

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


When η 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 η 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 η pockets let us say we have x number of red-antired, y number of blue-antiblue and z= η -x -y number of green-antigreen combinations, which will be represented by R, B and G, respectively. This is depicted in Table <ref>.




When two R pockets interact for example, there are four color combinations as depicted in the left box of Table <ref>. Out of these, rr̅ and r̅r belong to interaction type CAA, rr to CCR and r̅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 <ref>. 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 <ref> 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 η and y from 0 to η -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.

    E_i = ∑_x=0^η∑_y=0^η-xη!/x! y! z!×
       { x(x-1)+y(y-1)+z(z-1)/2(RR)_i + ( xy + yz + zx ) (RB)_i + η(R)_i } .

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^η×η (2η -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 η(2η-1), i.e.,

    ∑_i    { x(x-1)+y(y-1)+z(z-1)/2(RR)_i + ( xy + yz + zx ) (RB)_i + η(R)_i }
       = η( 2η -1 ),

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



After solving Eq. (<ref>) using Mathematica, we obtained the outcomes, E_i. They are listed in Table <ref>. In the same table, we divided each event by 3^ηη (2η -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 <ref>, we get -4/3g^2/(2η -1), very similar to the baryon case<cit.>. 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. <cit.>. We are interested in confined matter and will need to add a pressure term to the energy to enforce this; see Eq. (<ref>) below.

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

    Q⃗=∑_i=1^2ηq⃗_i .

Squaring on both sides, we get

    Q⃗^2      =∑_i=1^2ηq⃗_i^ 2 +2 ∑_i jq⃗_i·q⃗_j 
       
    =2η×4/3g^2 + 2×1/3^η∑_i E_i C_i 
       
    =0.

This shows our model is an overall color singlet. Note that we divided by 3^η 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 <ref> by three as shown in Table <ref>. 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 <ref>, we can see that we correctly obtain ∑_i≤ j^3P_ij+P̅_ij+P̅̅̅_i j =1 for the sum of all the probabilities. 





§ SYSTEM ENERGIES AND EQUATIONS


Our multi quark-pair system consists of an equal number of η quarks and η 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 N̅^I as number of anti-quarks with anti-flavor index I such that

    ∑_Ig^I N^I =η,

and

    ∑_Ig̅^IN̅^I =η.

In terms of quark density these equations can be expressed as

    ∑_I ∫ d^3r  n^I_i(r) =η /3,

and

    ∑_I ∫ d^3r n̅^I_i(r) =η /3.

In Eqs.(<ref>) and (<ref>) degeneracy factors are already included in quark densities n^I_i(r) and 
n̅^I_i(r). 

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

    ∑_i ∫ d^3r  n^I_i (r) = N^Ig^I

and

    ∑_i ∫ d^3r n̅^I_i (r) =N̅^Ig̅^I

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

    n̂^I_i ≡3n^I_i/N^I

    and

    n̂̅̂^I_i ≡3n̅^I_i/N̅^I.

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η total number of particles, the number of interactions possible is η(2η -1). Out of these CCA and CCR types occur 
    ( ∑_I N^I( N^I-1)/2(g^I)^2 + ∑_I≠ JN^I N^J/2g^Ig^J + ∑_IN^Ig^I(g^I-1)/2)
 times,
AAR and AAA occurs 
    ( ∑_I N̅^I( N̅^I-1)/2(g̅^I)^2 + ∑_I≠ JN̅^IN̅^J/2g̅^Ig̅^J + ∑_IN̅^Ig̅^I(g̅^I-1)/2)
 times,
and CAR and CAA occurs 
    ( ∑_I,JN̅^I N^Jg̅^Ig^J)
 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

    ∑_I N^I( N^I-1)/2(g^I)^2 +    ∑_I≠ JN^I N^J/2g^Ig^J + ∑_IN^Ig^I(g^I-1)/2
       
    = ∑_I N^Ig^I/2( N^Ig^I -1) + ∑_I≠ JN^I N^J/2g^Ig^J
       
    = ∑_I (N^Ig^I)^2/2  -∑_I (N^Ig^I)/2 +∑_I≠ JN^I N^J/2g^Ig^J
       
    =1/2(∑_I N^Ig^I)^2  -∑_I (N^Ig^I)/2
       
    =1/2(η^2-η).

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

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

    T=   ∑_i,I∫^r_maxd^3r (6π^2ħ^3N^In̂_i^I(r)    )^5/3/20π^2ħ^3m^I(g^I)^2/3 +∑_i,I∫^r_maxd^3r (6π^2ħ^3N̅^In̂̅̂_i^I(r)    )^5/3/20π^2ħ^3m̅^I(g̅^I)^2/3.

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 <ref>, we can now define total potential energy to be

    U= 4/3g^2×
       {∑_IN^I( N^I-1)/2∫∫d^3r d^3r'(  ∑_iP_iin̂_i^I(r)n̂_i^I(r') -1/2∑_i<jP_ijn̂_i^I(r)n̂_j^I(r'))/|r⃗-r⃗ '|.
    
       
    +∑_I≠ JN^I N^J/2∫∫d^3r d^3r'(  ∑_iP_iin̂_i^I(r)n̂_i^J(r') -1/2∑_i<jP_ijn̂_i^I(r)n̂_j^J(r'))/|r⃗-r⃗ '|
       
    +∑_IN^Ig^I( g^I-1)/2(g^I)^2∫∫d^3r d^3r'(  ∑_iP_iin̂_i^I(r)n̂_i^I(r') -1/2∑_i<jP_ijn̂_i^I(r)n̂_j^I(r'))/|r⃗-r⃗ '|
       
    +∑_IN̅^I( N̅^I-1)/2∫∫d^3r d^3r'(  ∑_iP̅̅̅_iin̂̅̂_i^I(r)n̂̅̂_i^I(r') -1/2∑_i<jP̅̅̅_ijn̂̅̂_i^I(r)n̂̅̂_j^I(r'))/|r⃗-r⃗ '|
       
    +∑_I≠JN̅^IN̅^J/2∫∫d^3r d^3r'(  ∑_iP̅̅̅_iin̂̅̂_i^I(r)n̂̅̂_i^J(r') -1/2∑_i<jP̅̅̅_ijn̂̅̂_i^I(r)n̂̅̂_j^J(r'))/|r⃗-r⃗ '|
       
    +∑_IN̅^Ig̅^I( g̅^I-1)/2(g̅^I)^2∫∫d^3r d^3r'(  ∑_iP̅̅̅_iin̂̅̂_i^I(r)n̂̅̂_i^I(r') -1/2∑_i<jP̅̅̅_ijn̂̅̂_i^I(r)n̂̅̂_j^I(r'))/|r⃗-r⃗ '|
       
    -.∑_I,JN̅^I N^J∫∫d^3r d^3r'(  ∑_iP̅_iin̂̅̂_i^I(r)n̂_i^J(r') -1/2∑_i<jP̅_ijn̂̅̂_i^I(r)n̂_j^J(r'))/|r⃗-r⃗ '|}.

In Eq. (<ref>) n̂_i^I are number densities, which are normalized to one when integrated over space<cit.>. 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 <ref>, 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 (<ref>). 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

    E=∑_i,I∫^r_maxd^3r(6π^2ħ^3N^In̂_i^I(r)  )^5/3/20π^2ħ^3m^I(g^I)^2/3 +∑_i,I∫^r_maxd^3r(6π^2ħ^3N̅^In̂̅̂_i^I(r)    )^5/3/20π^2ħ^3m̅^I(g̅^I)^2/3
       
    -4/3g^2∑_I,JN̅^I N^J∫∫d^3r d^3r'(  ∑_iP̅_iin̂̅̂_i^I(r)n̂_i^J(r') -1/2∑_i<jP̅_ijn̂̅̂_i^I(r)n̂_j^J(r'))/|r⃗-r⃗ '|.

Now, we switch back to normalization n^I_i and n̅^I_i and assume equal Fermi color momenta, n^I≡ n^I_1≡ n^I_2≡ n^I_3 for each I. The same assumption is made for antiparticles. This gives

    E=∑_I∫^r_max
    d^3r 3( 6π^2ħ^3N^In^I(r) )^5/3/20π^2ħ^3m^I(g^I)^2/3
    +   
       ∑_I∫^r_max
    d^3r3(6π^2ħ^3N̅^In̅^I(r)    )^5/3/ 20π^2ħ^3m̅^I(g̅^I)^2/3
    - 9×4/3 g^2/2η -1∑_I,J∫∫d^3r d^3r'/|r⃗-r⃗ '|n̅^I (r) n^J(r') .

Eqs. (<ref>) and (<ref>) can now be averaged over colors as 

    ∫ d^3r  n^I(r) = N^Ig^I/3,

and

    ∫ d^3r n̅^I(r) = N̅^Ig̅^I/3.

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



§ THOMAS FERMI QUARK EQUATIONS


We can now introduce Lagrange's undetermined multipliers λ^I and λ̅^I associated with constraints in Eqs.(<ref>) and (<ref>) and add them to the expression for energy. This gives

    E=   ∑_I∫^r_max
      d^3r 3( 6π^2ħ^3N^In^I(r) )^5/3/20π^2ħ^3m^I(g^I)^2/3
      +   ∑_I∫^r_max
      d^3r3(6π^2ħ^3N̅^In̅^I(r)    )^5/3/ 20π^2ħ^3m̅^I (g̅^I)^2/3
       
      - 9×4/3 g^2/2η -1∑_I,J∫∫d^3r d^3r'/|r⃗-r⃗ '|n̅^I    (r) n^J(r') 
     
       
     + ∑_I λ^I( 3∫ d^3r  n^I(r)-N^Ig^I)
     + ∑_Iλ̅^I( 3∫ d^3r n̅^I(r) -N̅^Ig̅^I).

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 δn̅^I (r) and δ n^I(r) give

    (6 π^2ħ^3)^5/3/π^2 ħ^3[1/4m̅^I(n̅^I(r)/g^I)^2/3] = -3 λ̅^I + 9 ×4/3 g^2/(2 η -1)∑_ I∫ d^3 r'/|r⃗-r⃗ '|n^I(r').

 
    (6 π^2ħ^3)^5/3/π^2 ħ^3 [1/4m^I(n^I(r)/g^I)^2/3] = -3 λ^I + 9 ×4/3 g^2/(2 η -1)∑_I∫ d^3 r'/|r⃗-r⃗ '|n̅^I(r'),

 
Assuming spherical symmetry, the TF spatial functions f^I(r) and f̅^I(r) are defined such that

    f^I(r)≡ra/( 8α_s/3)( 6π^2n^I(r)/g^I)^2/3,


    f̅^I(r)≡ra/(8α_s/3)( 6π^2n̅^I(r)/g̅^I)^2/3,

where

    a≡ħ/ m^1 c,

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

    α_s=g^2/ħ c,

is the strong coupling constant. Equation (<ref>) after using Eqs.(<ref>) - (<ref>) can now be written as

    3m^1/m̅^I4/3g^2 f̅^I(r)/r =     - 3 λ̅^I + 6·4/3g^2/(2η-1)π(2×4/3α_s/a)^3/2
       ×∑_I g^I[  1/r∫_0^r dr'r'^2(f^I(r')/r')^3/2   +∫_0^r_mdr'r'(f^I(r')/r')^3/2].

Here we have used the integral

    ∫^r_max d^3r'n^I(r')/|r⃗-r⃗ '| = 4π[ ∫_0^r dr' r'^2 n^I(r')/r+∫_r^rmax dr' r'^2 n^I(r')/r'],

resulting from spherical symmetry. Equation (<ref>) can be further simplified into

    α̅^I f̅^I (r) =-λ̅^I r/4/3g^2 +   2/(2η-1)π(2×4/3α_s/a)^3/2
       ×∑_I g^I[∫_0^r dr'r'^2(f^I(r')/r')^3/2  +r∫_0^r_mdr'r'(f^I(r')/r')^3/2].

At this stage we introduce the mass ratios

    α^I ≡m^1/m^I,   α̅^I ≡m^1/m̅^I,

for both the particle and antiparticle cases.

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

    R≡a/(8α_s/3)(3πη/2)^2/3.

Eq. (<ref>) can now be written as

    α̅^If̅^I (x)=   -λ^IRx/4/3g^2+
       3η/2η -1∑_I g^I[∫_0^x dx' x'^2(f^I(x')/x')^3/2  +x∫_x^x_mdx' x'(f^I(x')/x')^3/2].


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

    α^I f^I (x)=   -λ^IRx/4/3g^2+
       3η/2η -1∑_Ig̅^I[∫_0^x dx' x'^2 (f̅^I(x')/x')^3/2+ x∫_x^x_m dx' x'(f̅^I(x')/x')^3/2].


Taking first derivatives of Eqs. (<ref>) and (<ref>), we have

    α̅^I df̅^I(x)/dx = -λ̅^IR/4/3g^2 +3η/(2η-1)[∑_Ig^I∫_x^x_maxdx' x'(f^I(x')/x')^3/2],

and


    α^I df^I(x)/dx = -λ^IR/4/3g^2 +3η/(2η-1) [∑_Ig̅^I∫_x^x_mdx' x' (f̅^I(x')/x')^3/2].


Similarly, second derivatives of Eqs. (<ref>) and (<ref>) yield

    α̅^ Id^2f̅^ I(x)/dx^2 = - 3η/(2η -1)1/√(x)∑_I g^If^I(x)^3/2,

and

    α ^I d^2f^I(x)/dx^2 = - 3η/(2η -1)1/√(x)∑_Ig̅^If̅^I(x)^3/2.


Equations (<ref>) and (<ref>) 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=f̅ to reduce the TF differential equations into one incredibly simple form:

    α ^I d^2f^I(x)/dx^2 = - 3η/(2η -1)1/√(x)∑_I g^I f^I(x)^3/2.

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.




§ 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. (<ref>) and apply Eqs. (<ref>), (<ref>) and (<ref>). This gives the kinetic energy as

    T=   ∑_I12/5π( 3πη/2)^1/34/3g^2 ·4/3α_s/aα_I g^I∫_0^x_Idx(f^I(x))^5/2/√(x)
       
    + ∑_I12/5π( 3πη/2)^1/34/3g^2 ·4/3α_s/aα̅_I g̅^I∫_0^x_Idx(f̅^I(x))^5/2/√(x).

The potential energy

    U =    - 9·4/3g^2/( 2η -1 )∑_I,Jg̅^Ig^J∫^r_I∫^r_J d^3r   d^3r'  n^I(r) n^J(r')/|r⃗-r⃗ '|,

becomes

    U=
       -9·4/3g^2/( 2η -1 )η^2/R∑_I,J  g̅^Ig^J[   ∫_0^x_Idx (f̅^I(x))^3/2/√(x)∫_0^xdx' √(x')(f^J(x'))^3/2
        +  ∫_0^x_Idx (f̅^I(x))^3/2√(x)∫_x^x_Jdx'  (f^J(x'))^3/2/√(x')].


The volume energy (E_v) term gives the inward pressure from the vacuum. We assume that<cit.> 

    E_v= 4/3π R^3 x_max^3  B,

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

    E= T+ V+ E_v + η· m_q +η·m̅_q,

where m_q and m̅_q are the mass of the quark and anti-quark respectively.

In Section <ref> we will fit model parameters for a given set of mesons, including m_q and 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.




§ 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.





 §.§ 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 QQ̅, QQ̅QQ̅, QQ̅QQ̅QQ̅ 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<cit.>. Eq. (<ref>) becomes


    d^2f(x)/dx^2 = - 3η/(2η -1)· g_0·1/√(x)f(x)^3/2,

where

    g_0 × N^I=η.

We choose the normalization equation to be

    ∫_0^x_maxdx√(x)f(x)^3/2 =N^I/3η.

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

    (xdf/dx-f )|_x_max=-η/2η -1.


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

    T= 2·12/5π( 3πη/2)^1/34/3g^2 ·4/3α_s/a g_0 ∫_0^x_maxdx(f(x))^5/2/√(x).

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

    T= 24/5π( 3πη/2)^1/34/3g^2 ·4/3α_s/a[ -5/21df(x)/dx|_x_max + 4/7√(x_max)(f(x_max))^5/2 g_0   ].

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. (<ref>) can be written as

    U=-9·4/3g^2/( 2η -1 )η^2/R   (g^I)^2 [   ∫_0^x_maxdx (f(x))^3/2/√(x)∫_0^xdx' √(x')(f(x'))^3/2
        +  ∫_0^x_maxdx (f(x))^3/2√(x)∫_x^x_maxdx'  (f(x'))^3/2/√(x')].

This can be further simplified to
 
    U= 4/π( 3πη/2)^1/34/3g^2 ·4/3α_s/a×[ 4/7df(x)/dx|_x_max - 4/7√(x_max)(f(x_max))^5/2 g_0   ],

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.







 §.§ Z-meson family: Case 2

The constituents of Z-mesons are charm (c) (or bottom(b)), anti-charm (c̅) (or anti-bottom (b̅)), light (u or d) and anti-light quarks. The particles in this family can be represented as QQ̅qq̅, QQ̅qq̅QQ̅qq̅, QQ̅qq̅QQ̅qq̅QQ̅qq̅, 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=f̅ in the TF equations and obtain Eq. (<ref>) 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

    g_1 N_1 + g_2 N_2 =η.

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.,

    f^1(x) = k f^2(x)    for   0≤ x ≤ x_2,  
        f^2(x) = 0           for   x_2 ≤ x ≤ x_1.

 




Eq. (<ref>) can now be written as two set of equations:

    α ^1 d^2f^1(x)/dx^2 = - 3η/(2η -1)1/√(x)( g_1(f^1(x))^3/2  + g_2(f^2(x))^3/2),


    α ^2 d^2f^2(x)/dx^2 =- 3η/(2η -1)1/√(x)( g_1(f^1(x))^3/2  +g_2(f^2(x))^3/2).

To make Eqs. (<ref>) and (<ref>) consistent, we need

    k = α ^2/α ^1,

which is just the inverse ratio of the given masses from (<ref>). The similar step in the case of baryons gives a much more complicated consistency condition<cit.>. The normalization conditions are

    ∫_0^x_2x^1/2(f^2(x))^3/2 dx = N_2/3η,

and

    ∫_0^x_1x^1/2(f^1(x))^3/2 dx = N_1/3η.

In region 0<x<x_2,

    d^2f^1(x)/dx^2 = Q_1 (f^1(x))^3/2/√(x),

where

    Q_1 = - 3η/( 2η -1)α_1( g_1 + g_2/k^ 3/2).

In region x_2<x<x_1,

    d^2f^1(x)/dx^2 = Q_2 (f^1(x))^3/2/√(x),

where

    Q_2 = - 3η/( 2η -1)α_1g_1 .

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

    .(   xdf^2(x)/dx - f^2(x)   )|_x_2 = - N_2 Q_1 k^3/2/3η,


    .(   xdf^1(x)/dx - f^1(x)   )|_x_1 = - η/( 2η -1 ) α_1.



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

    T= 2·12/5π( 3πη/2)^1/34/3g^2 ·4/3α_s/a   [ g_1 α_1 ∫_0^x_1(f^1(x))^5/2/√(x)dx+
        g_2 α_2 ∫_0^x_2(f^2(x))^5/2/√(x)dx].

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:

    T= 2·12/5π   ( 3πη/2)^1/34/3g^2 ·4/3α_s/a[ -5/21α_1 df^1(x)/dx|_x_1
       + 4/7√(x_1)(f^1(x_1))^5/2 g_1 α_1  +  4/7√(x_2)(f^2(x_2))^5/2 g_2 α_2  ].

The expression for potential energy in the same case becomes

    U=-9·4/3g^2/( 2η -1 )η^2/R   ∑_I,Jg^Ig^J[   ∫_0^x_Idx (f^I(x))^3/2/√(x)∫_0^xdx' √(x')(f^J(x'))^3/2
        +  ∫_0^x_Idx (f^I(x))^3/2√(x)∫_x^x_Jdx'  (f^J(x'))^3/2/√(x')].

which can be placed into the form

    U=- 9·4/3g^2/( 2η -1 )η^2/R[   g_1^2 K_1 + g_2^2K_2+ g_1g_2K_12+ g_1g_2K_21] ,

where

    K_1 ≡∫_0^x_1dx (f^1(x))^3/2/√(x)∫_0^xdx' √(x')(f^1(x'))^3/2 +  
    ∫_0^x_1dx (f^1(x))^3/2√(x)∫_x^x_1dx'  (f^1(x'))^3/2/√(x'),


    K_2 ≡∫_0^x_2dx (f^2(x))^3/2/√(x)∫_0^xdx' √(x')(f^2(x'))^3/2 + 
    ∫_0^x_2dx (f^2(x))^3/2√(x)∫_x^x_2dx'  (f^2(x'))^3/2/√(x') ,


    K_12≡∫_0^x_1dx (f^1(x))^3/2/√(x)∫_0^xdx' √(x')(f^2(x'))^3/2  +  
    ∫_0^x_1dx (f^1(x))^3/2√(x)∫_x^x_2dx'  (f^2(x'))^3/2/√(x') ,

and

    K_21≡∫_0^x_2dx (f^2(x))^3/2/√(x)∫_0^xdx' √(x')(f^1(x'))^3/2  +  
    ∫_0^x_2dx (f^2(x))^3/2√(x)∫_x^x_1dx'  (f^1(x'))^3/2/√(x').

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

    K_1   = ∫_0^x_2 dx ( f^1(x) )^5/2/√(x)(-1/Q_1) 
       
        + ∫_0^x_2 dx √(x)( f^1 (x) )^3/2( 1/Q_1. df^1(x)/dx|_x_2 + . 1/Q_2df^1(x)/dx|_x_1 -. 1/Q_2df^1(x)/dx|_x_2) 
       
        +∫_x_2^x_1 dx ( f^1(x))^3/2/√(x)(   N_2 k^3/2/η( 1- Q_1/Q_2)  )
       
            +∫_x_2^x_1 dx ( f^1(x))^5/2/√(x)(  -1/Q_2)
       
                +∫_x_2^x_1 dx ( f^1(x))^3/2√(x)(  .  1/Q_2df^1(x)/dx|_x_1),


    K_2     = ∫_0^x_2 dx ( f^2(x) )^5/2/√(x)(  - 1/Q_1 k^1/2) 
       
         +∫_0^x_2 dx ( f^2(x))^3/2√(x)(  . 1/Q_1 k^1/2df^2(x)/dx|_x_2),

and
  
    K_12    =  ∫_0^x_2 dx ( f^2(x) )^5/2/√(x)(- k/Q_1) 
       
         + ∫_0^x_2 dx ( f^2(x))^5/2√(x)( .k/Q_1df^2(x)/dx|_x_2)
       
         +∫_x_2^x_1( f^1(x) )^3/2/√(x)( N_2/η).

Similarly, the potential energy now simplifies to
 
    U   =-9·4/3g^2/( 2η -1 )η^2/R×
        [ ∫_0^x_2dx(f^2(x))^5/2/√(x)[-k^5/2/Q_1g_1^2 - 1/Q_1 k^1/2g_2^2  - k/Q_1 2 g_1 g_2] . 
       
          +     ∫_0^x_2dx (f^2(x))^3/2√(x)  ×
       [ g_1^2(. k^5/2/Q_1df^2(x)/dx|_x_2 +.k^3/2/Q_2df^1(x)/dx|_x_1 - .k^5/2/Q_2df^2(x)/dx|_x_2) . 
       .+ .g_2^2 1/Q_1 k^1/2df^2(x)/dx|_x_2    +  .2g_1 g_2 k/Q_1df^2(x)/dx|_x_2]  
       
          +∫_x_2^x_1 dx (f^1(x))^3/2/√(x)[ N_2k^3/2/η( 1- Q_1/Q_2) g_1^2 + N_2/η2g_1g_2] 
       .  +∫_x_2^x_1 dx (f^1(x))^3/2√(x)[.1/Q_2df^1(x)/dx|_x_1g_1^2] + .
       .  ∫_x_2^x_1(f^1(x))^5/2/√(x)[-1/Q_2g_1^2] ],

which on further reduction yields
  
    U   = -9·4/3g^2/( 2η -1 )η^2/R×
       [.N_2/ηdf^2(x)/dx|_x_2×[ 5/7k^5/2g_1^2/Q_1 + 5/7g_2^2/Q_1 k^1/2 + 5/7k/Q_1 2g_1g_2  +  k^5/2g_1^2/Q_1 -  k^5/2g_1^2/Q_2..
       .  +g_2^2/Q_1 k^1/2 + k/Q_1 2g_1g_2  -k^5/2/Q_2( 1-Q_1/Q_2) g_1^2 - k/Q_2 2g_1g_2  - 5/7Q_2 Q_1 k^5/2g_1^2/Q_2]  
       
         + df^1(x)/dx×[ N_2/ηk^3/2g_1^2/Q_2+ N_2 k^3/2/η Q_2(1- Q_1/Q_2) g_1^2 +N_2/η Q_2 2g_1g_2. 
       .+ g_1^2/Q_2^2(   -3η/(2η-1)α_1  - N_2  Q_1 k^3/2/η)  -g_1^2/Q_2^2(5/7-3η/(2η-1)α_1)] 
        + 4/7√(x_2)(f^2(x_2))^5/2[    -k^5/2g_1^2/Q_1  -    g_2^2/Q_1 k^1/2 -  k/Q_1 2g_1g_2  +  k^5/2g_1^2/Q_2] 
       . +4/7√(x_1)(f^1(x_1))^5/2[- g_1^2/Q_2]].

Here we observe that the coefficient of N_2/η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

    U = -4/π( 3πη/2)^1/34/3g^2 ·4/3α_s/a×
       [ .-4/7α_1 df^1(x)/dx|_x_1 + 4/7√(x_1)(f^1(x_1))^5/2 g_1 α_1  +4/7√(x_2)(f^2(x_2))^5/2 g_2 α_2  ].


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.







 §.§ 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 Qq̅, Qq̅Qq̅, Qq̅Qq̅Qq̅, 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 f̅, respectively. We assume there is a universal TF function 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

    d^2f(x)/dx^2 = Q_0 (f(x))^3/2/√(x),

where

    Q_0 = - 3η g_0/( 2η -1)· k_0,

and

    k_0=(g_0 α/g̅_0 α̅)^2/5,

is the consistency condition. α is given by Eq. (<ref>), and in this case α̅=1.

In region x_2<x<x_1 with two TF equations, one function is zero and the other becomes

    d^2f̅(x)/dx^2 =  0, 
       f̅(x)=c· x +d ,

where

    c=k_0f'(x_2),

and

    d =k_0(f(x_2)-f'(x_2)· x_2 ) .

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

    T=    12/5π( 3πη/2)^1/34/3g^2 ·4/3α_s/a×
       [ . -10/21αdf(x)/dx|_x_2 + 8/7√(x_2)(f(x_2))^5/2 g_0α   +  g̅_0∫_x_2^x_1(f̅(x))^5/2/√(x)dx],


    U=    - 4/π( 3πη/2)^1/34/3g^2 ·4/3α_s/a×
       [.-4/7αdf(x)/dx|_x_2g̅_0/g_0 + 4/7√(x_2)(f(x_2))^5/2g̅_0α +η g̅_0/2η -1∫_x_2^x_1(f̅(x))^3/2/√(x)dx].


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 η=1 in Case 3 are equal to that of Case 1. Also, the non-degenerate Case 2 with η=2 is the same as degenerate Case 1 with η=2. We confirmed this both in our analytical and numerical results.

 

§ METHOD AND REMARKS



The phenomenological parameters we need for our model are the strong coupling constant α_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<cit.>. In addition, we did a fitting using mesonic states that involves only the charm quarks<cit.>. 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 removethis 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 

    1/4(η_c(1S)) + 3/4(J/Ψ(1S)) = 3069  MeV


    1/4(η_b(1S)) + 3/4(Y(1S)) = 9445  MeV

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<cit.>. The ground state of this set is the χ_c1(3872), with a mass of 3871.7 MeV, discovered by the Belle Collaboration in 2003<cit.>. However, this state appears to be a isospin singlet, although the charge states might not have been seen yet<cit.>. 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^0D̅^*+ 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^*0D̅^*+ 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 χ_c1(3872), Z_c(3900) and X(4020) states has a spin interaction Hamiltonian given by

    H_spin= κ_1 (s_1(s_1+1)-3/2)    + κ_2 (s_2(s_2+1)-3/2) + κ_3 (J(J+1)
       
    -s_1(s_1+1) -s_2(s_2+1)),

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 κ_2<κ_1, κ_3 based on the quark masses. By neglecting κ_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 χ_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<cit.>.

The Case 2 calculations are now simply,

    1/4× Z_c(3900) +3/4× X(4020) = 3990  MeV,

and
 
    1/4× Z_b(10610) +3/4× Z_b(10650) = 10641  MeV,

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

    D( spin zero average)=1/2(D^+) + 1/2(D^0) = 1866.5  MeV,

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

    1/4D( spin zero average) + 3/4D^*( spin one average) = 1973  MeV.

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 √(χ^2)= 86.1 MeV spread over 6 masses. We solved the differential equations using an iterative implementation of NDSolve in Mathematica. We obtained α_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<cit.>. Our baryon paper found α_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 <ref>. Table <ref> 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.


  


§ 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. <ref> and Fig. <ref>, 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. <ref> and Fig. <ref>. 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. <ref> and  Fig. <ref> 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. <ref> and Fig. <ref> 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<cit.> in D and B mesons is a capable tool in particle physics and has been implemented in various studies<cit.>. It arises when the heavy quark appears to the light quark degrees of freedom as a static color source<cit.>. 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. <ref> and <ref>, 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. <ref> and <ref> 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. <ref> and Fig. <ref> 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. (<ref>) we can see that the radius is proportional to the product of η^2/3 and the dimensionless radius x. As the quark content increases, the dimensionless radius, x, becomes smaller, whereas η^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. <ref> and Fig. <ref>. 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. <ref> and Fig. <ref> 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. <ref> and Fig. <ref> 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η), 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. <ref> and <ref> 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. <ref> and <ref>. 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.  has two-point function results for quark flavor content udc̅c̅ with spin 1 and uuc̅c̅ content for spin 0 for dynamical pion masses ranging from 257-688 MeV. Their continuum extrapolations show binding of -23.3± 11.4 MeV for the spin 1 sector and antibinding (resonance energy) of 25.9± 10.9 MeV above threshold in the spin 0 sector. On the other hand, the udb̅b̅ with spin 1 has a reported binding of -143.3± 33.9 MeV and uub̅b̅ has a small spin 0 binding of -5.5± 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<cit.>. The best we can do at this stage is to point out that our Case 3 results from Fig. <ref> 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 ∼ 2σ (double charm case) or ∼ 1σ (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. <cit.> gives

    E_m=8π/39×4/3g^2/(2η -1)∑_ I< Jγ_ Iγ̅_ J∑_i=1^N_ I∑_j=1^N_ J   (S_z^ I)_i(S̅_z^ J)_j  ∫ d^3r (   n^ I_i(r)n̅^ J_j(r)/N_ IN̅_ J) .

The sums are over generalized flavor indices I, J (which include spin) as well as particle indices i,j. The (S_z^ I)_i, (S̅_z^ J)_j are the spin components and the γ_ I,γ̅_ J factors in <ref> above are given by

    γ_ I=g/2m_ Ic, γ̅_ J=g/2m̅_ Jc,

where m_ I,m̅_ J is the quark or antiquark mass and g is the color gyromagnetic factor. Note that the energy shift in <ref> is for each TF configuration; as in <cit.> 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. <cit.>, with newly fit meson parameters will allow us to evaluate systematic uncertainties within the TF quark model approach in future spectroscopic work.





§ 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. (<ref>). 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.









































