text:
# QCD string in light-light and heavy-light mesons

## I Introduction

The description of the mass spectrum of hadrons is one of the fundamental problems of strong interactions. It has been attacked in a sequence of approaches motivated by QCD, but still attracts considerable attention. One of the most intriguing phenomena — namely, the formation of an extended object, the QCD string, between the colour constituents inside hadrons, — plays a crucial role in understanding their properties. In the present paper this role is exemplified by spectra of the mass of light-light and heavy-light mesons. In the former case we study the role played by the QCD string in formation of the straight-line Regge trajectories and discuss the form of the interquark interaction inside light hadrons. For heavy-light mesons we find the masses of several low-lying states in the $`D`$, $`D_s`$, $`B`$, and $`B_s`$ meson spectra including orbitally and radially excited ones.

We calculate and discuss the spin-spin and spin-orbit splittings and compare them to experimental and recent lattice data. Special attention is payed to the role of the proper string dynamics in establishing the correct slope of the Regge trajectories for both light-light and heavy-light states, as opposed to those following from relativistic equations with local potentials.

We remind the reader then that an extra piece of the effective interquark potential, the string correction, which is entirely due to the string-type interaction in QCD , gives a negative contribution to the masses of orbitally excited states. The latter observation allows us to resolve the “mystery” of an extremely narrow $`D(2637)`$ state (and a similar one in the $`B`$-mesonic spectrum)  recently claimed by the DELPHI Collaboration . We present a reasonable fit for the several lowest states in $`D`$\- and $`B`$-mesonic spectra using the standard values for the string tension, the strong coupling constant, and the current quark masses. We also find the correspondence between our model and Heavy Quark Effective Theory, extracting the constants used in the latter approach in the expansion of a heavy-light meson mass in the inverse powers of the heavy quark mass. We find analytical formulae for these constants and compare their numerical estimates with the predictions of other models.

The two main approaches used in the numerical calculations are the quasiclassical method of solving the eigenenergies problem and the variational one based on the einbein field formalism. The accuracy of both methods is tested using exactly solvable equations and found to be about 7% at worst even for the lowest states. Possible improvements of the method are outlined and discussed.

The paper is organized as follows. In Section II we give a brief insight into various aspects of the einbein field formalism. In Section III the exact spectra of relativistic equations are compared to the results of approximate calculations using the quasiclassical and variational einbein field methods, as well as the combined one. In Section IV we discuss the problem of the Regge trajectory slopes as they appear from the relativistic equations with local potentials and from the string-like picture of confinement. Derivation of the Hamiltonian for the spinless quark-antiquark system as well as of the spin-dependent corrections to it is the subject of Section V. The spectra of light-light and heavy-light mesonic states are calculated and discussed in Sections VI and VII, respectively. Section VIII contains our conclusions and outlook.

## II Einbein field formalism

In this section we give a short introduction to the method of einbein fields and its possible applications to relativistic systems. The interested reader can find a more detailed information in  and references therein.

### A Reparametrization invariance and constrained systems

Historically the einbein field formalism was introduced in  to treat the kinematics of the relativistic spinless particles. Later it was generalized to the case of spinning particles  and strings . So the action of a free relativistic particle can be rewritten as***In the path integral formalism this transformation is based on the following relation
$$D\mu (\tau )\mathrm{exp}\left(𝑑\tau \left(\frac{a\mu }{2}+\frac{b}{2\mu }\right)\right)\mathrm{exp}\left(𝑑\tau \sqrt{ab}\right).$$

$$S=_{\tau _i}^{\tau _f}L(\tau ),L=m\sqrt{\dot{x}^2}\frac{m^2}{2\mu }\frac{\mu \dot{x}^2}{2},$$
(1)
where the dot denotes derivative with respect to the proper time $`\tau `$, $`\mu `$ being the einbein fieldUsually $`e=\frac{1}{\mu }`$ is referred to as the einbein .. The original form of the action can be easily restored after solving the Euler–Lagrange equation of motion for the einbein $`\mu `$ which amounts to taking the extremum in the latter. Note that the invariance of the initial action with respect to the change of the proper time,

$$\tau f(\tau ),\frac{df}{d\tau }>0,f(\tau _i)=\tau _i,f(\tau _f)=\tau _f,$$
(2)
is preserved if an appropriate rescaling is prescribed to $`\mu `$:

$$\mu \mu /\dot{f}.$$
(3)
The latter invariance means that one deals with a constrained system. For the free particle the only constraint defines the mass shell,

$$p^2m^2=0,$$
(4)
or, in presence of the einbein field $`\mu `$,

$$\pi =0,H=\frac{p^2m^2}{2\mu },$$
(5)
with $`\pi `$ being the momentum canonically conjugated to $`\mu `$ and $`H`$ being the Hamiltonian function of the system (in case (4) it identically vanishes). The requirement that the constraint $`\pi =0`$ be preserved in time returns one to the mass-shell condition (4):

$$0=\dot{\pi }=\{\pi H\}=\frac{H}{\mu }=\frac{p^2m^2}{2\mu ^2}p^2m^2.$$
(6)
To make things simpler, one can fix the gauge-like freedom (2) identifying the proper time $`\tau `$ with one of the physical coordinates of the particle. The most popular choices are

* the laboratory gauge ($`\tau =x_0`$);
* the proper time gauge ($`\tau =(nx)`$, $`n_\mu =\frac{P_\mu }{\sqrt{P^2}}`$ with $`P_\mu `$ being the total momentum of the system) ;
* the light-cone gauge ($`\tau =\frac{1}{2}(x_0+x_3)=x_+`$),

which lead to quantization of the system on different hypersurfaces.

With the laboratory gauge fixed the Lagrangian function (1) becomes

$$L=\frac{m^2}{2\mu }\frac{\mu }{2}+\frac{\mu \dot{\stackrel{}{x}}^2}{2},$$
(7)
so that the corresponding Hamiltonian function reads

$$H=\frac{\stackrel{}{p}^2+m^2}{2\mu }+\frac{\mu }{2},$$
(8)
and after taking the extremum in $`\mu `$ one ends with the standard relativistic expression

$$H=\sqrt{\stackrel{}{p}^2+m^2}.$$
(9)
### B Einbeins as variational parameters

In the simple example considered above neither the Lagrange nor the Hamilton functions of the system contained $`\dot{\mu }`$, which allowed one to get rid of $`\mu `$ at any stage by taking the extremum in the latter. It is not so for more complicated systems when a change of variables is to be performed which touches upon the einbeins. The velocity corresponding to the original degrees of freedom of the system may mix in a very tangled way with those for einbeins, so that it is not a simple task anymore to follow the lines $`\stackrel{´}{a}`$ la Dirac  to resolve the set of constraints and to get rid of nonphysical degrees of freedom. See, $`e.g.`$,  for several examples when such a resolution can be done explicitly.

Luckily another approach to einbeins is known . They can be treated as variational parameters. Thus one replaces the dynamical function of time $`\mu (\tau )`$ by the parameter $`\mu _0`$ independent of $`\tau `$. The eigenstate problem is solved then, keeping $`\mu _0`$ constant, so that one has the spectrum $`M_{\{n\}}(\mu _0)`$, where $`\{n\}`$ denotes the full set of quantum numbers. Then one is to minimize each eigenenergy independently with respect to $`\mu _0`$Note that solutions for $`\mu _0`$ of both signs appear, but only one of them ($`\mu _0>0`$) is finally left. Neglecting the negative solution is the general lack of the einbein field approach and this leads to the fact that quark Zitterbewegung is not taken into account (see also the discussion in Subsection III D).:

$$\frac{M_{\{n\}}(\mu _0)}{\mu _0}|_{\mu _0=\mu _0^{}}=0,M_{\{n\}}=M_{\{n\}}(\mu _0^{}).$$
(10)
Such an approach has a number of advantages. First, it allows one to avoid the tedious algebra of commuting constraints with one another following the standard Dirac technique . Second, it allows one a very simple and physically transparent interpretation of einbeins. Indeed, in formulae (1) and (7) the einbein $`\mu `$ can be treated as an effective mass of the particle; the dynamics of the system remains essentially relativistic, though being nonrelativistic in form. If $`m`$ is the current quark mass, then $`\mu `$ can be viewed as its constituent mass celebrated in hadronic phenomenology. What is more, the current mass can be even put to zero, whereas the Lagrangian approach remains valid in the presence of the einbeins and the standard Hamiltonian technique can be developed then. The latter observation is intensively used in analytic QCD calculations for glue describing gluonic degrees of freedom in glueballs and hybrids .

An obvious disadvantage of the variational approach to the einbein fields is some loss of accuracy. As a variational method it provides only an approximate solution, giving no hint as to how to estimate the ultimate accuracy of the results. Thus in the next section we test this method, comparing its predictions with exact solutions of some relativistic equations. We consider the accuracy, found to be about 7% at worst, quite reasonable, which justifies our consequent attack on the light-light and heavy-light mesons spectra using this formalism.

## III Testing the method

### A Quasiclassics for the spinless Salpeter equation

We start from the Salpeter equation for the quark-antiquark system with equal masses and restrict ourselves to the zero-angular-momentum case for simplicity:

$$\left(2\sqrt{p_r^2+m^2}+\sigma r\right)\psi _n=M_n^{(ll)}\psi _n,$$
(11)
where the subscript $`(ll)`$ stands for the light-light system.

The quasiclassical quantization condition looks like

$$_0^{r_+}p_r(r)𝑑r=\pi \left(n+\frac{3}{4}\right),n=0,1,2,\mathrm{},r_+=\frac{M_n^{(ll)}2m}{\sigma },$$
(12)
where the integral on the l.h.s. can be worked out analytically, yielding

$$M_n^{(ll)}\sqrt{\left(M_n^{(ll)}\right)^24m^2}4m^2\mathrm{ln}\frac{\sqrt{\left(M_n^{(ll)}\right)^24m^2}+M_n^{(ll)}}{2m}=4\sigma \pi \left(n+\frac{3}{4}\right),$$
(13)
or approximately $`(m\sqrt{\sigma })`$ one has

$$\left(M^{(ll)}\right)^2=4\pi \sigma \left(n+\frac{3}{4}\right)+2m^2\mathrm{ln}\frac{\pi \sigma (n+3/4)}{m^2}+\mathrm{}.$$
(14)
Solution (14) becomes exact in the limit $`m=0`$, whereas for a nonzero mass the leading correction to the linear regime $`(M^{(ll)})^2n`$ behaves like

$$\frac{\mathrm{\Delta }M_n^2}{M_n^2}=O\left(\frac{m^2}{M_n^2}\mathrm{ln}\frac{M_n}{m}\right)\underset{n1}{}\frac{\mathrm{ln}n}{n}.$$
(15)
For a heavy-light system one has the Salpeter equation

$$\left(\sqrt{p_r^2+m^2}+\sigma r\right)\psi _n=M_n^{(hl)}\psi _n,$$
(16)
where $`M_n^{(hl)}`$ denotes the excess over the heavy particle mass $`M`$. Similarly to (13) one finds then

$$M_n^{(hl)}\sqrt{\left(M_n^{(hl)}\right)^2m^2}m^2\mathrm{ln}\frac{\sqrt{\left(M_n^{(hl)}\right)^2m^2}+M_n^{(hl)}}{m}=2\sigma \pi \left(n+\frac{3}{4}\right),$$
(17)
and formula (15) holds true in this case as well.

Comparing the results of the WKB method with the exact solutions of equation (11) (rows $`M_n(\mathrm{WKB})`$ and $`M_n(\mathrm{exact})`$ in Table II), one can see that the error does not exceed 3-4% even for the ground state. See also  where the WKB method is tested for light-light mesons.

### B Quasiclassics for the one-particle Dirac equation

As a next example we discuss the one-particle Dirac equation with linearly rising confining potential :

$$(\stackrel{}{\alpha }\stackrel{}{p}+\beta (m+U)+V)\psi _n=\epsilon _n\psi _n.$$
(18)
The WKB method applied to this equation gives

$$_r_{}^{r_+}\left(p+\frac{\kappa w}{pr}\right)𝑑r=\pi \left(n+\frac{1}{2}\right),n=0,1,2,\mathrm{},$$
(19)
where

$$p=\sqrt{(\epsilon V)^2\frac{\kappa ^2}{r^2}(m+U)^2},$$
(20)
$$w=\frac{1}{2r}\frac{1}{2}\frac{U^{}V^{}}{m+U+\epsilon V},$$

$$|\kappa |=j+\frac{1}{2}.$$

For the most interesting case of purely scalar confinement ($`V=0`$, $`U=\sigma r`$) an approximate quasiclassical solution was found in  ($`m=0`$):

$$\epsilon _n^2=2\sigma \left(2n+j+\frac{3}{2}+\frac{sgn\kappa }{2}+\frac{\kappa \sigma }{\pi \epsilon _n^2}\left(0.38+ln\frac{\epsilon _n^2}{\sigma |\kappa |}\right)+O\left(\left(\frac{\kappa \sigma }{\epsilon _n^2}\right)^2\right)\right).$$
(21)
A detailed comparison of the results of the WKB method and those following from the recursive formula (21) with exact numerical solutions to equation (18) is given in . Here we only note that the coincidence of the three numbers is impressive as even for the lowest states the discrepancy does not exceed 1%.

### C Quasiclassical variational einbein field (combined) method for the spinless Salpeter equation

Finally we combine the two methods discussed above and apply the WKB approximation to the Hamiltonian of a relativistic system with einbeins introduced as variational parameters. Then the resulting quasiclassical spectrum is minimized with respect to the einbeins. Thus we have a powerful method of solving the eigenvalues problem for various relativistic systems which we call “combined.” Let us test the accuracy of this method first.

We start from the Salpeter equation (11) for the light-light system and introduce the parameter $`\mu _0`$ as described in Section II:

$$H_1=2\sqrt{p_r^2+m^2}+\sigma rH_2=\frac{p_r^2+m^2}{\mu _0}+\mu _0+\sigma r.$$
(22)
In what follows we consider the massless case substituting $`m=0`$ into (22).

We give the analytic formulae for the spectrum of the Salpeter equation (11) obtained using the quasiclassical approximation for the Hamiltonian $`H_1`$ (following from equation (13) for $`m=0`$), the exact solution for the Hamiltonian $`H_2`$ minimized with respect to the einbein field and the result of the combined method when the Bohr-Sommerfeld quantization condition is applied to the Hamiltonian $`H_2`$ and the ultimate spectrum is also minimized with respect to $`\mu _0`$.

$$M_n^2(\mathrm{WKB})=4\pi \sigma \left(n+\frac{3}{4}\right),$$
(23)
$$M_n^2(\mathrm{einbein})=16\sigma \left(\frac{\zeta _{n+1}}{3}\right)^{3/2},$$
(24)
$$M_n^2(\mathrm{combined})=\frac{8\pi }{\sqrt{3}}\sigma \left(n+\frac{3}{4}\right),$$
(25)
where $`\zeta _{n+1}`$ is the $`(n+1)`$th zero of the Airy function $`Ai(z)`$ and counting of zeros starts from unity. The extremal values of the einbein field in the latter two cases read

$$\mu _0^{}(\mathrm{einbein})=\sqrt{\sigma }\left(\frac{\zeta _{n+1}}{3}\right)^{3/4},\mu _0^{}(\mathrm{combined})=\sqrt{\frac{\sigma (n+3/4)}{2\sqrt{3}}},$$
(26)
$`i.e.`$, the effective quark mass is $`\mu _0^{}\sqrt{\sigma }`$ and it appears entirely due to the interquark interaction.

In Table II we compare the results of the above three approximate methods of solving the eigenvalues problem for equation (11) with the exact solution. In the last row we give the accuracy of the combined method vs the exact solution. Two conclusions can be deduced from Table II. The first one is that the accuracy of all approximate methods is high enough, including the combined method, which is of most interest for us in view of its consequent applications to the QCD string with quarks at the ends. The other conclusion is that the variational einbein field method gives a systematic overestimation for the excited states which is of order 5-7%.

### D Discussion

Here we would like to make a couple of concluding comments concerning the numerical methods tested in this section, their accuracy and possible ways of their improvement. As stated above the combined quasiclassical variational method is of most interest for us, so we shall concentrate basically on it. The following two remarks are in order here.

From Table II one can see that the relative error is practically constant, tending to the value of 7% for large $`n`$. The reason for such a behaviour will become clear if one compares formulae (23) and (25). Both relations reproduce the same dependence on the radial quantum number $`n`$, whereas the difference comes from different slopes, $`4\pi \sigma `$ in (23) $`vs`$ $`\frac{8\pi }{\sqrt{3}}\sigma `$ in (25). Then for highly excited states the error is practically independent of $`n`$ and can be estimated as

$$\delta =\frac{M_n(\mathrm{combined})M_n(\mathrm{exact})}{M_n(\mathrm{combined})}\frac{M_n(\mathrm{combined})M_n(\mathrm{WKB})}{M_n(\mathrm{combined})}=1\sqrt{\frac{\sqrt{3}}{2}}0.07;$$
(27)
$`i.e.`$, the ultimate accuracy of the quasiclassical variational einbein field method (combined method) appears to be about 7%. Introducing, say, a correcting factor in (25) one could overcome the systematic overestimation and reproduce the spectrum with a better accuracy. We shall return to this observation later on when discussing the spectrum of the heavy-light mesons.

Another source of error in the einbein field approach is neglecting the quark Zitterbewegung (see the footnote on page II B). As stated above we neglect the negative sign solution for the einbein field $`\mu _0`$ expecting its small influence on the spectrum. Let us give some reasoning to justify this action.

It was demonstrated numerically in  that the contribution of the quark Z-graphs into $`M^2`$ is nearly constant for large excitation numbers and is of order 10%, so that the corresponding shift of $`M`$ behaves like

$$\mathrm{\Delta }M\frac{\mathrm{\Delta }M^2}{M}\underset{n1}{}\frac{1}{\sqrt{n}},$$
(28)
so it is somewhat suppressed.

Besides, the good agreement of our numerical results with those provided by the lattice data and taken from the Particle Data Group can also serve as an a posteriori justification of such a neglect. Still some improvements for the einbein field approach are needed to take this effect into account.

## IV Regge trajectories for relativistic equations with local potentials

It was observed long ago that the mesonic Regge trajectories are almost linear if the total momentum or the radial quantum number is plotted vs the mesonic mass squared :

$$M^2(n,J)=c_nn+c_JJ+\mathrm{\Delta }M^2,$$
(29)
where $`c_n`$ and $`c_J`$ are the (inverse) slopes while $`\mathrm{\Delta }M^2`$ denotes corrections to the leading linear regime which come from the self-energy, spin splittings, $`etc`$.

Relations like (29) naturally appear in most of models for confinement, though the (inverse) slopes $`c_n`$ and $`c_J`$ are different for different models.

For example the Salpeter equation for the heavy-light system,

$$\left(\sqrt{p_r^2+p_{\theta ,\phi }^2+m^2}+\sigma r\right)\psi _{nl}=M_{nl}\psi _{nl},$$
(30)
gives

$$c_J^{(hl)}(Salpeter)=4\sigma ,c_n^{(hl)}(Salpeter)=2\sigma ,$$
(31)
where the total momentum $`J`$ coincides with the orbital one $`l`$.

For the light-light system one easily finds, from (31) by a trivial parameters rescaling,

$$c_J^{(ll)}(Salpeter)=8\sigma ,c_n^{(ll)}(Salpeter)=4\sigma .$$
(32)
The one-particle Dirac equation (18) yields different slopes for different natures of the confining force. Thus for the purely vector confinement (potential added to the energy term§§§We disregard the problem of the Klein paradox here.) one finds

$$c_{J\mathrm{vec}}^{(hl)}(Dirac)=4\sigma ,$$
(33)
whereas, for purely scalar confinement (potential added to the mass term),

$$c_{J\mathrm{scal}}^{(hl)}(Dirac)=2\sigma .$$
(34)
In the meantime the spectrum (29) is expected to follow from a string-like picture of confinement which predicts the (inverse) Regge slopes to be

$$c_J^{(ll)}(string)=2\pi \sigma ,c_J^{(hl)}(string)=\pi \sigma .$$
(35)
One can easily see that none of the relativistic equations considered before gives the correct result (35); moreover the discrepancies are rather large (of order 25%). See also  for a discussion of various models of confinement and the corresponding Regge trajectory slopes.

The reason why relativistic Salpeter and Dirac equations fail to reproduce the correct string slope of the Regge trajectories is obvious and quite physically transparent. Indeed, all relativistic equations with local potentials have only a trivial dependence of the interquark interaction on the angular momentum which comes entirely from the quark kinetic energy. Meanwhile, QCD is believed to lead to a string-type interaction between the colour constituents inside hadrons, whereas the QCD string developed between quarks possesses its own inertia and thus it should also contribute to the $`J`$-dependent part of the interaction. It is this extra purely string-type piece of the interquark interaction to give an extra contribution to the Regge trajectory slope and to bring it into the correct form of (35). This statement is proved explicitly in the next section, whereas the string dynamics footprint in the heavy-light mesons spectrum is discussed in detail in Section VII.

## V Hamiltonian of the $`q\overline{q}`$ meson

### A Quark-antiquark Green’s function

We start from the Euclidean Green’s function of the $`q\overline{q}`$ pair in the confining vacuum

$$G_{q\overline{q}}=\mathrm{\Psi }_{q\overline{q}}^{(f)}(\overline{x},\overline{y}|A)^+\mathrm{\Psi }_{q\overline{q}}^{(i)}(x,y|A)_{q\overline{q}A},$$
(36)
where the initial and the final mesonic wave functions

$$\mathrm{\Psi }_{q\overline{q}}^{(i,f)}(x,y|A)=\overline{\mathrm{\Psi }}_{\overline{q}}(x)\mathrm{\Phi }(x,y)\mathrm{\Gamma }^{(i,f)}\mathrm{\Psi }_q(y)$$
(37)
are gauge invariant due to the standard path-ordered parallel transporter

$$\mathrm{\Phi }(x,y)=P\mathrm{exp}\left(ig_y^x𝑑z_\mu A_\mu \right),$$
(38)
$`\mathrm{\Gamma }^{(i,f)}`$ denote the matrices which might be inserted into the initial and final meson-quark-antiquark vertices.

Integrating out the quark fields in (36), one finds, for the mesonic Green’s function,

$$G_{q\overline{q}}=Tr\mathrm{\Gamma }^{(f)}S_q(\overline{x},x|A)\mathrm{\Phi }(x,y)\mathrm{\Gamma }^{(i)}S_{\overline{q}}(y,\overline{y}|A)\mathrm{\Phi }(\overline{y},\overline{x})_A,$$
(39)
where the trace stands for both colour and spinor indices. We have neglected here the $`1/N_C`$-suppressed quark determinant, describing sea quark pairs.

To proceed further we employ the Feynman-Schwinger representation for the one-fermion propagators in the external field, fix the laboratory gauge for both particles,

$$x_{10}=t_1,x_{20}=t_2,$$
(40)
and introduce the einbein fields $`\mu _1`$ and $`\mu _2`$ by means of the following change of variables (see  for details):

$$\mu _i(t_i)=\frac{T}{2s_i}\dot{x}_{i0}(t_i),ds_iDx_{i0}D\mu _i(t_i),i=1,2,$$
(41)
where $`s_{1,2}`$ are the Schwinger times, $`T=\frac{1}{2}(x_0+\overline{x}_0y_0\overline{y}_0)`$.

Then the resulting expression for the mesonic Green’s function reads

$$G_{q\overline{q}}=D\mu _1(t_1)D\mu _2(t_2)D\stackrel{}{x}_1D\stackrel{}{x}_2e^{K_1K_2}Tr[\mathrm{\Gamma }^{(f)}(m_1\widehat{D})\mathrm{\Gamma }^{(i)}(m_2\widehat{D})\times $$
(42)
$$P_\sigma \mathrm{exp}\left(_0^T\frac{dt_1}{2\mu _1(t_1)}\sigma _{\mu \nu }^{(1)}\frac{\delta }{i\delta s_{\mu \nu }(x_1(t_1))}\right)\mathrm{exp}(_0^T\frac{dt_2}{2\mu _2(t_2)}\sigma _{\mu \nu }^{(2)}\frac{\delta }{i\delta s_{\mu \nu }(x_2(t_2))})\mathrm{exp}(\sigma S_{\mathrm{min}})],$$

with $`K_i`$ being the kinetic energies of the quarks

$$K_i=_0^T𝑑t_i\left(\frac{m_i^2}{2\mu _i}+\frac{\mu _i}{2}+\frac{\mu _i\dot{\stackrel{}{x}}_i^2}{2}\right),i=1,2,$$
(43)
$`\sigma _{\mu \nu }=\frac{1}{4i}(\gamma _\mu \gamma _\nu \gamma _\nu \gamma _\mu )`$, and $`\delta /\delta s_{\mu \nu }`$ denotes the derivative with respect to the element of the area $`S`$. We have also used the minimal area law asymptotic for an isolated Wilson loop,

$$TrP\mathrm{exp}\left(ig_C𝑑z_\mu A_\mu \right)_A\mathrm{exp}(\sigma S_{\mathrm{min}}),$$
(44)
which is usually assumed for the stochastic QCD vacuum (see, $`e.g.`$, ) and found on the lattice. Here $`S_{\mathrm{min}}`$ is the area of the minimal surface swept by the quark and antiquark trajectories.

Looking at (42) one can easily recognize the following three main ingredients: the contribution of the quark, the one of the antiquark, and finally the confining interaction given by the string with tension $`\sigma `$. One can write, for the latter,

$$S_{\mathrm{min}}=_0^T𝑑t_0^1𝑑\beta \sqrt{(\dot{w}w^{})^2\dot{w}^2w^2},$$
(45)
with $`w_\mu (t,\beta )`$ being the string profile function chosen in linear form,

$$w_\mu (t,\beta )=\beta x_{1\mu }(t)+(1\beta )x_{2\mu },$$
(46)
thus describing the straight-line string which is a reasonable approximation for the minimal surface .

Finally, synchronizing the quark and the antiquark times $`(t_1=t_2=t)`$ one finds from (42) that in the spinless approximation the quark-antiquark meson can be described by the Lagrangian

$$L(t)=\frac{m_1^2}{2\mu _1}\frac{m_2^2}{2\mu _2}\frac{\mu _1+\mu _2}{2}+\frac{\mu _1\dot{\stackrel{}{x}}_1^2}{2}+\frac{\mu _2\dot{\stackrel{}{x}}_2^2}{2}\sigma r_0^1𝑑\beta \sqrt{1[\stackrel{}{n}\times (\beta \dot{\stackrel{}{x}}_1+(1\beta )\dot{\stackrel{}{x}}_2)]^2},$$
(47)
where $`\stackrel{}{r}=\stackrel{}{x}_1\stackrel{}{x}_2`$ and $`\stackrel{}{n}=\stackrel{}{r}/r`$. Expansion of the surface-ordered exponents in (42) gives a set of spin-dependent corrections to the leading regime (47).

### B Hamiltonian for spinless quarks

Starting from the Lagrangian (47) and introducing an extra einbein field $`\nu (t,\beta )`$ continuously depending on the internal string coordinate $`\beta `$ one can get rid of the square root in (47) arriving at the Hamiltonian of the $`q\overline{q}`$ system in the centre-of-mass frame in the form

$$H=\underset{i=1}{\overset{2}{}}\left(\frac{p_r^2+m_i^2}{2\mu _i}+\frac{\mu _i}{2}\right)+_0^1𝑑\beta \left(\frac{\sigma ^2r^2}{2\nu }+\frac{\nu }{2}\right)+\frac{\stackrel{}{L}^2}{2r^2[\mu _1(1\zeta )^2+\mu _2\zeta ^2+_0^1𝑑\beta \nu (\beta \zeta )^2]},$$

$$\zeta =\frac{\mu _1+_0^1𝑑\beta \nu \beta }{\mu _1+\mu _2+_0^1𝑑\beta \nu }.$$
(48)
Similarly to $`\mu `$’s which have the meaning of the constituent quark masses, the einbein $`\nu `$ can be viewed as the density of the string energy. In the simplest case of $`l=0`$ one easily finds, for the extremal value of $`\nu `$,

$$\nu _0=\sigma r;$$
(49)
$`i.e.`$, the energy distribution is uniform and the resulting interquark interaction is just the linearly rising potential $`\sigma r`$. In the meantime, if $`l0`$, then the two contributions can be identified in the last $`l`$-dependent term in (48). Roughly speaking the first two $`\mu `$-dependent terms in the denominator come from the quark kinetic energy. The last term containing the integral over $`\beta `$ is nothing but the extra inertia of the string discussed before. Rotating string also contributes to the interquark interaction making it essentially nonlocal, so that the very notion of the interquark potential is not applicable to the system anymore.

Note that the Hamiltonian (48) has the form of sum of the “kinetic” and the “potential” parts, but this is somewhat misleading, as extrema in all three einbeins are understood, so that the ultimate form of the Hamiltonian would be extremely complicated and hardly available for further analytical studies.

Expression (48) can be simplified if one expands the Hamiltonian in powers of $`\sqrt{\sigma }/\mu `$. One finds then

$$H=H_0+V_{string},$$
(50)
$$H_0=\underset{i=1}{\overset{2}{}}\left(\frac{\stackrel{}{p}^2+m_i^2}{2\mu _i}+\frac{\mu _i}{2}\right)+\sigma r,$$
(51)
$$V_{string}\frac{\sigma (\mu _1^2+\mu _2^2\mu _1\mu _2)}{6\mu _1^2\mu _2^2}\frac{\stackrel{}{L}^2}{r},$$
(52)
where $`V_{string}`$ is known as the string correction  and this is the term totally missing in the relativistic equations with local potentials. Indeed, the Salpeter equation with the linearly rising potential is readily reproduced from (51) if extrema in $`\mu _{1,2}`$ are taken explicitly, whereas the string correction is lost. Meanwhile, its sign is negative so that the contribution of the string lowers the energy of the system, thus giving a negative contribution to the masses of orbitally excited states, leaving those with $`l=0`$ intact. In Section VI we shall demonstrate how a proper account of the string dynamics in the full Hamiltonian (48) brings the Regge trajectory slope to the correct value (35), whereas in Section VII the string correction (52) will be demonstrated to solve the problem of the identification of the resonance $`D(2637)`$ recently claimed by the DELPHI Collaboration .

### C Spin-dependent corrections

Let us return to the quark-antiquark Green’s function (42) and extract the nonperturbative spin-orbit interaction. Following  one finds

$$V_{so}^{\mathrm{np}}=\frac{\sigma }{2r}\left(\frac{\stackrel{}{S}_1\stackrel{}{L}}{\mu _1^2}+\frac{\stackrel{}{S}_2\stackrel{}{L}}{\mu _2^2}\right).$$
(53)
It follows from  that all potentials $`V_i(r)`$ (in the notation of ) contain both perturbative and nonperturbative pieces given there in explicit form. One can argue that at large distances only the piece (53) is left whereas for light quarks all nonperturbative ones may be important (see ).

Now, to have a full picture of the interquark interaction one is to supply the purely nonperturbative string-type interaction described by the Hamiltonian (48) by the perturbative gluon exchange adding the colour Coulomb potential to the Hamiltonian $`H_0`$ from (51) and calculating the corresponding spin-dependent perturbative terms in addition to the potential (53). The result reads

$$H_0=\underset{i=1}{\overset{2}{}}\left(\frac{\stackrel{}{p}^2+m_i^2}{2\mu _i}+\frac{\mu _i}{2}\right)+\sigma r\frac{4}{3}\frac{\alpha _s}{r}C_0,$$
(54)
where we have also added the overall constant shift $`C_0`$ and

$$V_{sd}=\frac{8\pi \kappa }{3\mu _1\mu _2}(\stackrel{}{S}_1\stackrel{}{S}_2)\left|\psi (0)\right|^2\frac{\sigma }{2r}\left(\frac{\stackrel{}{S}_1\stackrel{}{L}}{\mu _1^2}+\frac{\stackrel{}{S}_2\stackrel{}{L}}{\mu _2^2}\right)+\frac{\kappa }{r^3}\left(\frac{1}{2\mu _1}+\frac{1}{\mu _2}\right)\frac{\stackrel{}{S}_1\stackrel{}{L}}{\mu _1}+\frac{\kappa }{r^3}\left(\frac{1}{2\mu _2}+\frac{1}{\mu _1}\right)\frac{\stackrel{}{S}_2\stackrel{}{L}}{\mu _2}$$

$$+\frac{\kappa }{\mu _1\mu _2r^3}\left(3(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})(\stackrel{}{S}_1\stackrel{}{S}_2)\right)+\frac{\kappa ^2}{2\pi \mu ^2r^3}\left(\stackrel{}{S}\stackrel{}{L}\right)(2\mathrm{ln}(\mu r)\gamma _E),\gamma _E=0.57,$$
(55)
with $`\kappa =\frac{4}{3}\alpha _s`$, $`\mu =\frac{\mu _1\mu _2}{\mu _1+\mu _2}`$. We have also added the term of order $`\alpha _s^2`$ which comes from one-loop calculations and is intensively discussed in the literature . It is important to stress that $`C_0`$ is due to the nonperturbative self-energy of light quarks, which explains the later numerical inputs.

An important comment concerning the expansion (55) is in order. Up to the last term the expression (55) coincides in form with the Eichten–Feiberg–Gromes result , but we have effective quark masses $`\mu _i`$ in the denominators instead of the current ones $`m_i`$. Once $`\mu _i\sqrt{\stackrel{}{p}^2+m_i^2}>m_i`$ or even $`\mu _im_i`$, then the result (55) is applicable to the case of light quark flavours, when the expansion of the interaction in the inverse powers of the quark mass $`m_i`$ obviously fails. The values of $`\mu `$’s are defined dynamically and differ from state to state (see  for details).

The Hamiltonian (54) with spin-dependent terms (55) will be used for explicit calculations for heavy-light mesons. In the case of light-light states one should include additional nonperturbative spin-dependent terms (see  and references therein). The masses of light-light mesons listed in Table II have been calculated from the Regge trajectories which do not take into account spin-dependent terms and we give them for the sake of comparison. A more detailed calculation for the light mesons taking these effects into account can be found in .

## VI Spectrum of light-light mesons

### A Angular-momentum-dependent potential and Regge trajectories

Starting from the Hamiltonian (48) we stick with the case of equal quark masses $`m`$

$$H=\frac{p_r^2+m^2}{\mu }+\mu +\frac{\stackrel{}{L}^2/r^2}{\mu +2_0^1(\beta \frac{1}{2})^2\nu 𝑑\beta }+\frac{\sigma ^2r^2}{2}_0^1\frac{d\beta }{\nu }+_0^1\frac{\nu }{2}𝑑\beta .$$
(56)
The extremal value of the einbein field $`\nu `$ can be found explicitly and reads

$$\nu _0(\beta )=\frac{\sigma r}{\sqrt{14y^2(\beta \frac{1}{2})^2}},$$
(57)
where $`y`$ is the solution of the transcendental equation

$$\frac{L}{\sigma r^2}=\frac{1}{4y^2}(\mathrm{arcsin}yy\sqrt{1y^2})+\frac{\mu y}{\sigma r},$$
(58)
and $`\stackrel{}{L}^2=l(l+1)`$.

For large angular momenta the contribution of the quarks (the last term on the r.h.s. of (58)) is negligible so that the maximal possible value $`y=1`$ is reached, thus yielding the solution for the free open string  (see also the second entry in ).

With the extremal value $`\nu _0`$ from (57) inserted, the Hamiltonian (56) takes the form

$$H=\frac{p_r^2+m^2}{\mu (\tau )}+\mu (\tau )+\frac{\sigma r}{y}\mathrm{arcsin}y+\mu (\tau )y^2,$$
(59)
with $`y`$ defined by equation (58). The last two terms on the r.h.s. of equation (59) can be considered as an effective “potential”

$$U(\mu ,r)=\frac{\sigma r}{y}\mathrm{arcsin}y+\mu (\tau )y^2,$$
(60)
which is nontrivially $`l`$ dependent. In Fig.1 we give the form of the effective potential (60) for a couple of low-lying states (solid line). It has the same asymptotic as the naive sum of the linearly rising potential and the centrifugal barrier coming from the kinetic energy of the quarks (dotted line). In the meantime it differs from the latter at finite distances. The only exception is the case of zero angular momentum which should be treated separately and leads to the linearly rising potential for any interquark separation.

### B Numerical results

Following the variational einbein field method described and tested above, we start from the Hamiltonian (59) and change the einbein field $`\mu `$ for the variational parameter $`\mu _0`$ , so that one has

$$H=\frac{p_r^2+m^2}{\mu _0}+\mu _0+U(\mu _0,r),$$
(61)
$$U(\mu _0,r)=\frac{\sigma r}{y}\mathrm{arcsin}y+\mu _0y^2.$$
(62)
Then the quasiclassical method applied to the Hamiltonian (61) gives

$$_r_{}^{r_+}p_r(r)𝑑r=\pi \left(n+\frac{1}{2}\right),$$
(63)
with

$$p_r(r)=\sqrt{\mu _0(M\mu _0U(\mu _0,r))m^2}.$$
(64)
The eigenvalues $`M_{nl}(\mu _0)`$ for $`m=0`$ were found numerically from (63), (64) and the minimization procedure with respect to $`\mu _0`$ was used then. Results for $`M_{nl}`$ are given in Table II and depicted in Fig.2 demonstrating very nearly straight lines with approximately string slope $`(2\pi \sigma )^1`$ in $`l`$ and twice as small slope in $`n`$. Note that it is the region of intermediate values of $`r`$ that plays the crucial role in the Bohr-Sommerfeld integral (63), $`i.e.`$, the region where the nontrivial dependence of the effective potential $`U(\mu _0,r)`$ on the angular momentum is most important (see Fig.1).

In Table III we give a comparison of the masses of several light-light mesonic states extracted by means of the numerical results from Table II with the experimental data and theoretical predictions taken from . We have fitted our results to the experimental spectrum using the negative constant $`\mathrm{\Delta }M^2`$ (see equation (29)).

### C Discussion

Let us recall the results obtained for the light-light mesons and discuss problems connected to the given approach. The net result of the current section is the $`l`$-dependent effective interquark potential which gives the naive linearly rising interaction only for $`l=0`$. It was observed long ago  that for large angular momenta the quark dynamics is negligible and the slope (35) naturally appears from the picture of open rotating string. In the present paper we find that for massless quarks even the low-lying mesonic states demonstrate nearly straight-line Regge trajectories with string slope (35).

One problem clearly seen from Figs.2,3 is the leading trajectory intercept $`l_0l(M^2=0)`$. To reproduce the experimental intercept around -0.5 (see Fig.3) starting from the theoretical one +0.5 (see left plot in Fig.2) one needs a large negative constant added either to the Hamiltonian (48) (see, $`e.g.`$, $`C_0`$ in equation (54)) or in the form of $`\mathrm{\Delta }M^2`$ directly in (29) (see also Table III). Once the first way might violate the linearity of the Regge trajectories, then one should expect QCD to prefer the second one, though the first way remains more attractive from the practical point of view and will be used in calculations of the heavy-light mesons spectrum in the next section.

Another problem is that one of the most intriguing questions of mesonic spectroscopy, the $`\pi \rho `$ splitting (and a similar problem in the strange sector) cannot be addressed in our model. Taking the exact solution of the spinless Salpeter equation (11) with $`n=l=0`$ (see Table II with an appropriate rescaling from $`\sigma =0.2GeV^2`$ to $`\sigma =0.17GeV^2`$) one finds for the $`\rho `$ mass squared the value of order $`1.7GeV^2`$ which does not violate the linearity of the trajectory (see the circled dot in Fig.2). If the overall negative shift with $`\sqrt{|\mathrm{\Delta }M^2|}=1126MeV^2`$ (see Table III) is applied to this state, then one arrives at a $`\rho `$-meson mass about $`775MeV`$, $`i.e.`$, a value very close to the experimental one. Note that we have practically coinciding constants for the $`\rho `$\- and $`a`$-meson trajectories (see the caption for Table III), which supports the idea that $`\mathrm{\Delta }M^2`$ can be associated with quark self-energies.

Meanwhile, one cannot pretend to describe pions (kaons) in the same framework as their Goldstone nature is not implemented in the current model. In realistic quantum-field-theory-based models each mesonic state possesses two wave functions which describe the motion forward and backward in time of the $`q\overline{q}`$ pair inside the meson . The backward motion is suppressed if at least one of the quarks is heavy, for highly excited states and in the infinite-momentum frame. For the chiral pion, which is expected to be strictly massless in the chiral limit, the two wave functions are of the same order of magnitude (see, $`e.g.`$,  for an explicit pionic solution in QCD2), so that none of them can be neglected. This explains why the naive estimate for the pion mass lies much higher than the experimentally observed value of $`140MeV`$. For the first excited state, $`\rho `$ meson, this effect is already suppressed, though one still has to be careful to neglect the backward motion of the quarks. The progress in this direction was achieved in recent papers by one of the authors (Yu.S.) , where a Dirac-type equation was derived for the heavy-light system and the properties of its solutions were investigated. This new formalism is expected to allow consideration of pionic Regge trajectories as it has the chiral symmetry breaking built in.

## VII Spectrum of heavy-light mesons

All results obtained for the light-light mesons in Section VI can be reproduced for the heavy-light states, so that in the one-body limit the Regge trajectories with the correct string (inverse) slope $`\pi \sigma `$ are readily reproduced. Meanwhile, the aim of this study is to take into account corrections to the leading regime which come from the spin-dependent terms in the interquark interaction as well as those due to the finitness of the heavy quark mass. Corrections of both types are important for establishing the correct spectra of $`D`$ and $`B`$ mesons which are the main target of the present investigation.

### A Spectrum of the spinless heavy-light system

In this subsection we study the spectrum of the heavy-light mesons, disregarding the quark spins. This amounts to solving the Schrödinger-like equation for the Hamiltonian $`H_0`$ from (54). Note that to this end one needs to know the nonrelativistic spectrum in the potential which is the sum of the linearly rising and Coulomb parts :

$$\left(\frac{d^2}{d\stackrel{}{x}^2}+|\stackrel{}{x}|\frac{\lambda }{|\stackrel{}{x}|}\right)\chi _\lambda =a(\lambda )\chi _\lambda ,$$
(65)
where

$$\lambda =\kappa \left(\frac{2\mu }{\sqrt{\sigma }}\right)^{2/3},\kappa =\frac{4}{3}\alpha _s,\mu =\frac{\mu _1\mu _2}{\mu _1+\mu _2}.$$

If solutions of (65) for $`\chi _\lambda `$ and $`a(\lambda )`$ are known as functions of the reduced Coulomb potential strength $`\lambda `$, then one can find the following expressions for the extremal values of the einbeins (constituent quark masses):

$$\mu _1(\lambda )=\sqrt{m_1^2+\mathrm{\Delta }^2(\lambda )},\mu _2(\lambda )=\sqrt{m_2^2+\mathrm{\Delta }^2(\lambda )},\mu (\lambda )=\frac{1}{2}\sqrt{\sigma }\left(\frac{\lambda }{\kappa }\right)^{3/2},$$
(66)
with $`\mathrm{\Delta }(\lambda )`$ given by

$$\mathrm{\Delta }^2(\lambda )=\frac{\sigma \lambda }{3\kappa }\left(a+2\lambda \left|\frac{a}{\lambda }\right|\right).$$

The definition of the reduced einbein field $`\mu `$ via $`\mu _1`$ and $`\mu _2`$ leads to the equation defining $`\lambda `$

$$\mu (\lambda )=\frac{\mu _1(\lambda )\mu _2(\lambda )}{\mu _1(\lambda )+\mu _2(\lambda )}.$$
(67)
Technically this means that one should generate self-consistent solutions to equations (65) and (67) which are subject to numerical calculations . In Table IV we give such solutions for several radial and orbital excitations in $`D`$, $`D_s`$, $`B`$, and $`B_s`$mesonic spectra. We use the standard values for the string tension, the strong coupling constant, and the current quarks masses. Note that $`\alpha _s`$ is chosen close to its frozen value  and it does not change a lot between $`D`$ and $`B`$ mesons. The reason is that in both cases one has a light quark moving in the field of a very heavy one, so that the one-gluon exchange depends on the size of the system, rather than on its total mass. Once the difference in size between $`D`$ and $`B`$ mesons is not that large, the difference between the two values of the strong coupling constant is also small (see Table IV).

The $`\psi `$ function at the origin given in the last column of Table IV and which will be used later on for spin-spin splittings is defined for radially excited states as

$$\left|\psi (0)\right|^2=\frac{2\mu \sigma }{4\pi }\left(1+\lambda x^2\right),$$
(68)
where

$$r^N=(2\mu \sigma )^{N/3}x^N=(2\mu \sigma )^{N/3}_0^{\mathrm{}}x^{N+2}\left|\chi _\lambda (x)\right|^2𝑑x,N>32l,$$
(69)
which immediately follows from the properties of equation (65) and the corresponding redefinition of variables.

### B Spin-spin and spin-orbit splittings. The string correction

In this subsection we calculate the spin-dependent corrections to the results given in Table IV as well as those due to the proper string dynamics and which were intensively discussed before.

The eigenstates of the Hamiltonian $`H_0`$ from (54), which we consider to be the zeroth approximation, can be specified in the form of terms $`n^{2S+1}L_J`$ ($`n`$ being the radial quantum number) as the angular momentum $`\stackrel{}{L}`$, the total spin $`\stackrel{}{S}`$, and the total momentum $`\stackrel{}{J}=\stackrel{}{L}+\stackrel{}{S}`$ are separately conserved by $`H_0`$. The corresponding matrix elements for various operators present in (55) read as follows:

$$\begin{array}{cccc}{}_{}{}^{2S+1}P_{J}^{}\hfill & & & \\ & ^1P_1|\stackrel{}{S}_1\stackrel{}{L}|^1P_1=0,\hfill & ^1P_1|\stackrel{}{S}_2\stackrel{}{L}|^1P_1=0,\hfill & ^1P_1|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^1P_1=\frac{1}{4},\hfill \\ & & & \\ & ^3P_0|\stackrel{}{S}_1\stackrel{}{L}|^3P_0=1,\hfill & ^3P_0|\stackrel{}{S}_2\stackrel{}{L}|^3P_0=1,\hfill & ^3P_0|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3P_0=\frac{1}{4},\hfill \\ & & & \\ & ^3P_1|\stackrel{}{S}_1\stackrel{}{L}|^3P_1=\frac{1}{2},\hfill & ^3P_1|\stackrel{}{S}_2\stackrel{}{L}|^3P_1=\frac{1}{2},\hfill & ^3P_1|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3P_1=\frac{1}{4},\hfill \\ & & & \\ & ^3P_2|\stackrel{}{S}_1\stackrel{}{L}|^3P_2=\frac{1}{2},\hfill & ^3P_2|\stackrel{}{S}_2\stackrel{}{L}|^3P_2=\frac{1}{2},\hfill & ^3P_2|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3P_2=\frac{1}{20};\hfill \end{array}$$
(70)
$$\begin{array}{cccc}{}_{}{}^{2S+1}D_{J}^{}\hfill & & & \\ & ^1D_2|\stackrel{}{S}_1\stackrel{}{L}|^1D_2=0,\hfill & ^1D_2|\stackrel{}{S}_2\stackrel{}{L}|^1D_2=0,\hfill & ^1D_2|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^1D_2=\frac{1}{4},\hfill \\ & & & \\ & ^3D_1|\stackrel{}{S}_1\stackrel{}{L}|^3D_1=\frac{3}{2},\hfill & ^3D_1|\stackrel{}{S}_2\stackrel{}{L}|^3D_1=\frac{3}{2},\hfill & ^3D_1|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3D_1=\frac{1}{12},\hfill \\ & & & \\ & ^3D_2|\stackrel{}{S}_1\stackrel{}{L}|^3D_2=\frac{1}{2},\hfill & ^3D_2|\stackrel{}{S}_2\stackrel{}{L}|^3D_2=\frac{1}{2},\hfill & ^3D_2|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3D_2=\frac{1}{4},\hfill \\ & & & \\ & ^3D_3|\stackrel{}{S}_1\stackrel{}{L}|^3D_3=1,\hfill & ^3D_3|\stackrel{}{S}_2\stackrel{}{L}|^3D_3=1,\hfill & ^3D_3|(\stackrel{}{S}_1\stackrel{}{n})(\stackrel{}{S}_2\stackrel{}{n})|^3D_3=\frac{1}{28}.\hfill \end{array}$$
(71)
The interaction $`V_{sd}`$ given by (55) mixes orbitally excited states with different spins, so that the transition matrix elements are given by

$$\begin{array}{cc}^1P_1|\stackrel{}{S}_1\stackrel{}{L}|^3P_1=\frac{1}{\sqrt{2}},\hfill & ^1P_1|\stackrel{}{S}_2\stackrel{}{L}|^3P_1=\frac{1}{\sqrt{2}},\hfill \\ & \\ ^1D_2|\stackrel{}{S}_1\stackrel{}{L}|^3D_2=\sqrt{\frac{3}{2}},\hfill & ^1D_2|\stackrel{}{S}_2\stackrel{}{L}|^3D_2=\sqrt{\frac{3}{2}},\hfill \end{array}$$
(72)
which lead to mixing within $`|^1P_1`$, $`|^3P_1`$ and $`|^1D_2`$, $`|^3D_2`$ pairs so that the physical states are subject to matrix equations of the following type:

$$\left|\begin{array}{cc}E_1E& V_{12}\\ V_{12}^{}& E_2E\end{array}\right|=0.$$
(73)
Another important ingredient is the string correction given by (52) which leads to an extra negative shift for orbitally excited states

$$\delta M_l\frac{\sigma (\mu _1^2+\mu _2^2\mu _1\mu _2)}{6\mu _1^2\mu _2^2}l(l+1)r^1.$$
(74)
Thus the model is totally fixed and the only remaining fitting parameter is the overall spectrum shift $`C_0`$ which finally takes the following values:

$$C_0(D)=212MeV,C_0(D_s)=124MeV,C_0(B)=203MeV,C_0(B_s)=124MeV.$$
(75)
Note that $`C_0`$ does not depend on the heavy quark ($`C_0(D)C_0(B)`$, $`C_0(D_s)C_0(B_s)`$) and is completely defined by the properties of the light one. For states with two light quarks one would have an overall negative shift $`2C_0`$, which gives the contribution $`4C_0(M_nC_0)`$ to $`M_n^2`$. In case of the $`\rho `$ meson this provides a negative constant of order $`1GeV`$, $`i.e.`$, the right value needed to bring the theoretical intercept $`l_0`$ into the correct experimental one (see $`|\mathrm{\Delta }M^2|`$ in Table III).

### C Comparison with the experimental and lattice data. “Mystery” of the $`D(2637)`$ state

In Tables VI and VI we compare the results of our numerical calculations for the spectrum and splittings with experimental and recent lattice data as well as with the theoretical predictions from  and . The underlined figures in Table VI are considered as the most probable candidates for the experimentally observed values. One can see good agreement between our theoretical predictions and the experimental values, as well as with the lattice calculations . To demonstrate the relevance of the corrections due to the heavy mass we consider a simplified system containing one infinitely heavy particle and the light one having its real mass. The best fits for the experimental spectra with the results for such simplified systems are also given in Table VI in the column entitled $`M_{hl}^{(0)}`$ for comparison. One can easily see that corrections in the inverse powers of the heavy mass are strongly needed to reproduce the experimental spectrum with a reasonable accuracy and to remove degeneracy of $`S`$ states.

Now we are in the position to resolve the “mystery” of the $`D(2637)`$ state (and a similar one in the $`B`$-meson spectrum). This state was claimed recently by the DELPHI Collaboration , but once its quantum numbers were not defined, then there was a problem of the identification of this state. In most quark models (see also Table VI) the first radial excitation $`J^P=0^{}`$ lies approximately in the desired region of mass, but estimates of the width of such a state lead to a confusion, as all such estimates give values much larger than the width of about $`15MeV`$ reported by DELPHI. The only would-be way out of the problem is to identify this narrow state with orbital excitations with $`J^P`$ being $`2^{}`$ or $`3^{}`$. In spite of the fact that orbitally excited states are really narrower than the radially excited ones and can have a width compatible with the experimental value, the following two objections can be made : i) quark models predict orbitally excited mesons to be at least $`50MeV`$ heavier than needed, and ii) a neighboring slightly more massive state should be observed as well.

It follows from Table VI that we can remove both objections mentioned above (see also ). Indeed, one can easily see that orbitally excited states $`2^{}`$ and $`3^{}`$ lie even somewhat lower than the radial excitation $`0^{}`$. The reason for that is the negative string correction (74) for the orbitally excited states which comes from the proper dynamics of the string. Besides that, the single $`D`$-wave $`3^{}`$ state is an even more probable candidate for the role of the observed $`D(2637)`$ resonance than the lightest one from the pair of $`2^{}`$ states, so that the problem of the “missing state” is also avoided.

Note that our predictions $`D(2654)`$, $`D(2663)`$, and $`D(2664)`$ give larger masses compared to the experimental one. This must be a reflection of the general lack of the variational einbein field method ($`\mu `$ technique) discussed before, which gives slightly overestimated values for the spectrum of excited states (see Table II and the discussion in subsection III D).

### D A bridge to the Heavy Quark Effective Theory

In this subsection we discuss the correspondence between our model and the Heavy Quark Effective Theory (HQET) approach widely discussed in the literature (see  and references therein). We use the standard parametrization for the heavy-light meson mass,

$$M_{hl}=m_Q+\overline{\mathrm{\Lambda }}\frac{1}{2m_Q}(\lambda _1+d_H\lambda _2)+O\left(\frac{1}{m_Q^2}\right),$$
(76)
where $`m_Q`$ is the mass of the heavy constituent, and the coefficient $`d_H`$ describes the hyperfine splitting,

$$d_H=\{\begin{array}{cc}+3,\hfill & \mathrm{for}\mathrm{\hspace{0.33em}0}^{}\mathrm{states},\hfill \\ 1,\hfill & \mathrm{for}\mathrm{\hspace{0.33em}1}^{}\mathrm{states},\hfill \end{array}$$
(77)
whereas $`\overline{\mathrm{\Lambda }}`$, $`\lambda _1`$, and $`\lambda _2`$ are free parameters which are subject to theoretical investigation. The parameter $`\lambda _2`$ is directly connected to the splitting between $`1^3S_1`$ and $`1^1S_0`$ states and can be estimated from the experimental $`B`$-mesonic spectrum to be

$$\lambda _2=\frac{1}{4}(M_B^{}M_B)0.12GeV^2.$$
(78)
From Table VI one can easily find our prediction for $`\lambda _2`$

$$\lambda _20.16GeV^2,$$
(79)
which being slightly overestimated is still in reasonable agreement with the experimental value (78).

In the meantime our model allows direct calculation of the parameters $`\overline{\mathrm{\Lambda }}`$, $`\lambda _1`$, and $`\lambda _2`$ based on the Hamiltonian (54). We apply the variational procedure described above to the idealized system with $`m_1m_Q\mathrm{}`$ and $`m_20`$. This yields

$$\overline{\mathrm{\Lambda }}=\sqrt{\sigma }\left[\sqrt{\frac{\kappa }{\lambda _0}}a(\lambda _0)+\sqrt{\frac{\lambda _0}{12\kappa }\left(a(\lambda _0)+2\lambda _0\left|\frac{a}{\lambda _0}\right|\right)}\right],$$
(80)
where $`\kappa =\frac{4}{3}\alpha _s`$, and $`a(\lambda )`$ is the dimensionless eigenvalue introduced in (65) (see subsection VII A); $`\lambda _0`$ is a solution to the equation

$$\lambda ^2=\frac{4}{3}\kappa ^2\left(a+2\lambda \left|\frac{a}{\lambda }\right|\right).$$
(81)
Then one can extract the coefficient $`\lambda _2`$ from the first term in equation (55):

$$\lambda _2=\frac{4\pi \kappa }{3\mu _2}|\psi (0)|^2=\frac{2}{3}\sigma \kappa \left(1+\lambda _0x^2\right).$$
(82)
An analytical formula for $`\lambda _1`$ is also available, but it is rather bulky and we do not give it here. For $`\alpha _s=0.39`$ one can find the numerical solution to equation (81) to be $`\lambda _0=1.175`$. The corresponding values for $`\overline{\mathrm{\Lambda }}`$, $`\lambda _1`$, and $`\lambda _2`$ are given in the first column of Table VII where they are compared with the results of other approaches.

Another way to estimate the discussed constants is to find the best fit of the form

$$M_{fit}=m_Q+\overline{\mathrm{\Lambda }}+C_0\frac{\lambda _1}{2m_Q},$$
(83)
with $`\overline{\mathrm{\Lambda }}`$ and $`\lambda _1`$ being the fitting parameters and $`C_0=203MeV`$ taken from (75), for eigenvalues of the Hamiltonian (54) with $`m_Q`$ varied around the bottom quark mass $`m_b=4.8GeV`$ (we use the region $`4GeV<m_Q<6GeV`$). The coefficient $`\lambda _2`$ can be found using formula (82) with $`\lambda _0`$ changed for the exact solution for $`\lambda `$ taken from Table IV. Results are listed in the second column of Table VII.

One can see our figures to be in general agreement with those found in other approaches among which we mention the QCD sum rules method , the inclusive semileptonic $`B`$-meson decays , and the Dyson-Schwinger equation for the system of a light quark and a static antiquark . We find the parameter $`\lambda _1`$ to be rather sensitive to the strong coupling constant $`\alpha _s`$. For example, for $`\alpha _s=0.3`$ one has $`\lambda _1=0.38GeV^2`$ which should be confronted with the value $`\lambda _1=0.506GeV^2`$ from the first column of Table VII found for $`\alpha _s=0.39`$.

All our predictions for $`\lambda _2`$ exceed the value given by equation (78). The reason is the slightly overestimated value of $`\psi (0)`$ given by the variational einbein field approach.

One should appreciate the advantage of the einbein field method which allows one to obtain relatively simple analytical formulae for various parameters and to investigate their dependence on the strong coupling constant $`\alpha _s`$ (the dependence on the only dimensional parameter $`\sigma `$ is uniquely restored, giving $`\overline{\mathrm{\Lambda }}\sqrt{\sigma }`$ and $`\lambda _1\lambda _2\sigma `$).

## VIII Conclusions

In conclusion let us briefly recall the main results obtained in the present paper.

We use the model for the QCD string with quarks at the ends to calculate the spectra of light-light and heavy-light mesons. There are two main points in which we differ from other approaches to the same problem based on various relativistic Hamiltonians and equations with local potentials. The first point is that we do not introduce the constituent mass by hand. On the contrary, starting from the current mass we naturally arrive at the effective quark masses which appear due to the interaction. Moreover, the resulting effective mass is large enough even for the lightest quarks and lowest states in the spectrum, so that the spin-dependent terms in the interquark interaction can be treated as perturbations in most cases (except for pions and kaons) and thus accounted for in this way.

The second advantage of the method is that the dynamics of the QCD string naturally enters the game and it can be studied systematically. A proper account of this dynamics allows one to resolve several problems of the mesonic spectroscopy; namely, one can find that the rotating string lowers the masses of orbitally excited states, bringing the (inverse) slope of Regge trajectories to their correct values ($`2\pi \sigma `$ for light-light and $`\pi \sigma `$ for heavy-light states). In addition, this allows one to resolve the problem of the identification of the $`D(2637)`$ state recently claimed by the DELPHI Collaboration and which is known to lead to a contradiction between its small width, incompatible with the decay modes for the radial excitation, and its mass lying considerably lower than the values predicted by the quark models for orbitally excited states. In the meantime, taking into account the negative string correction contributing into the masses of orbitally excited mesons readily resolves this “mystery” for the $`D(2637)`$ state as well as for its counterpart in the spectrum of $`B`$ mesons.

For the heavy-light system we extract the constants $`\overline{\mathrm{\Lambda }}`$, $`\lambda _1`$, and $`\lambda _2`$ used in the framework of the Heavy Quark Effective Theory, for which we derive analytical formulae. We find our numerical results to be in agreement with those obtained from the experimental data and calculations in other approaches like QCD sum rules, inclusive semileptonic $`B`$-meson decays, and the relativistic Dyson-Schwinger equation for the $`q\overline{Q}`$ system.

We also conclude that the string-like interaction favoured by QCD invalidates the very notion of any local interquark potential, still leaving room to the einbein field method for the Hamiltonian approach to the bound states of quarks and gluons in QCD. Being rather accurate, this method still requires improvements to increase the accuracy and to have the full control over it.

The authors are grateful to A.M.Badalian, P.Bicudo, R.N.Faustov, A.B.Kaidalov, V.S.Popov and E.Ribeiro for useful discussions and valuable comments and to V.L.Morgunov and B.L.G.Bakker for providing numerical solutions to some equations. One of the authors (A.N.) would like to thank the staff of the Centro de Física das Interacções Fundamentais (CFIF-IST) for cordial hospitality during his stay in Lisbon.

Financial support of RFFI grants 00-02-17836 and 00-15-96786, INTAS-RFFI grant IR-97-232 and INTAS CALL 2000, project # 110 is gratefully acknowledged. One of the authors (A.N.) is also supported via RFFI grant 01-02-06273.


format: text

text:
# Hadronization properties of b quarks compared to light quarks in 𝒆⁺𝒆⁻→𝒒𝒒̄ from 𝟏𝟖𝟑 to 𝟐𝟎𝟎 GeV

## 1 Introduction

The study of the properties of the fragmentation of heavy quarks compared to light quarks offers new insights in perturbative QCD. Particularly important is the difference in charged particle multiplicity between light quark and heavy quark initiated events in $`e^+e^{}`$ annihilations.

In a first approximation one could expect that the multiplicity of hadrons produced in addition to the possible decay products of the primary quark-antiquark is a universal function of the available invariant mass; this would give a difference in charged particle multiplicity between light quark and heavy quark initiated events decreasing with the centre-of-mass energy $`E_{cm}`$ . QCD predicts, somehow counter-intuitively, that this difference is energy independent; this is motivated by mass effects on the gluon radiation (see  and  for a recent review).

The existing experimental tests were not conclusive (see  and references therein, ). At LEP 2 energies, however, the difference between the QCD prediction and the model ignoring mass effects is large, and the experimental measurement can firmly distinguish between the two hypotheses.

## 2 Analysis and Results

A description of the DELPHI detector can be found in ; its performance is discussed in .

Data corresponding to a luminosity of 54 pb-1 collected by DELPHI at centre-of-mass (c.m.) energies around 183 GeV during 1997, to 158 pb-1 collected around 189 GeV during 1998, and to 187 pb-1 collected between 192 and 200 GeV during 1999, were analysed.

The 1999 data were taken at different energies: 25.8 pb-1 at 192 GeV, 77.4 pb-1 at 196 GeV and 83.8 pb-1 at 200 GeV. Each energy was analyzed separately and the results were then combined as described later and attributed to a c.m. energy of 200 GeV.

A preselection of hadronic events was made, requiring at least 10 charged particles with momentum $`p`$ above 100 MeV/$`c`$ and less than 1.5 times the beam energy, with an angle $`\theta `$ with respect to the beam direction between 20 and 160, a track length of at least 30 cm, a distance of closest approach to the interaction point less than 4 cm in the plane perpendicular to the beam axis and less than (4/$`\mathrm{sin}\theta `$) cm along the beam axis, a relative error on the momentum measurement $`\mathrm{\Delta }p/p<1`$, and a total transverse energy of the charged particles above 0.2$`E_{cm}`$.

The influence of the detector on the analysis was studied with the full DELPHI simulation program, DELSIM . Events were generated with PYTHIA 5.7 and JETSET 7.4 , with parameters tuned to fit LEP1 data from DELPHI . The Parton Shower (PS) model was used. The particles were followed through the detailed geometry of DELPHI giving simulated digitisations in each subdetector. These data were processed with the same reconstruction and analysis programs as the real data.

The hadronic cross-section for $`e^+e^{}`$ interactions above the Z peak is dominated by radiative q$`\overline{\mathrm{q}}\gamma `$ events; the initial state radiated photons (ISR photons) are generally aligned along the beam direction and not detected. In order to compute the hadronic c.m. energy, the procedure described in  was used. In this procedure particles are clustered into jets and the effective centre-of-mass energy of the hadronic system, $`\sqrt{s^{}}`$, is computed as being the invariant mass of the system recoiling against an ISR photon, possibly unseen.

Events with reconstructed hadronic c.m. energy ($`\sqrt{s^{}}`$) above 0.9$`E_{cm}`$ were used. The selected 1997 (1998, 1999) data sample consisted of 1699 (4583, 4881) hadronic events.

For each year’s data, two samples enriched in (1) $`b`$ events and in (2) $`uds`$ events were selected from the $`b`$ tagging variable $`y`$ defined as in Ref. ; this variable represents essentially the probability that none of the tracks in the event comes from a vertex separated from the primary one. To select the samples of the type (2), it was required in addition that the narrow jet broadening $`B_{min}`$ is smaller than 0.065, to reduce the background due to WW and ZZ events. $`B_{min}`$ is defined as follows. The event is separated into two hemispheres $`H_1`$ and $`H_2`$ with respect to the thrust axis, defined by the thrust unit vector $`\widehat{t}`$. Then, calling $`\stackrel{}{p_k}`$ the momentum vector of the $`k`$-th particle,

$$B_{min}=\underset{i=1,2}{\mathrm{min}}\frac{_{kH_i}\left|\stackrel{}{p_k}\times \widehat{t}\right|}{2_k\left|\stackrel{}{p_k}\right|}.$$

The contamination from non-q$`\overline{\mathrm{q}}`$ events in the samples of type (1) was 7% (8%, 15%), while it was 13% (17%, 20%) in the samples of type (2). After applying the event selection criteria and the cuts to reduce the WW and ZZ background, the purities were approximately 91% (90%, 90%) ($`b`$ events) over the total q$`\overline{\mathrm{q}}`$ in sample (1), and 79% (79%, 79%) ($`uds`$ events) over the total q$`\overline{\mathrm{q}}`$ in sample (2). The fractions of $`q`$-type quarks in the $`\left(i\right)`$-th sample, $`f_q^{\left(i\right)}`$, were determined from the simulation. The sample (1) consisted of 103 (326, 416) events; the sample (2) of 590 (1450, 1652) events.

The average charge multiplicity was measured in the samples (1) and (2), after subtracting the background bin-by-bin by means of the simulation. It should be noted that the average multiplicity for a given flavour $`q`$ in each sample is equal to $`C_q^{\left(i\right)}\times n_{q\overline{q}}`$, with $`C_q^{\left(i\right)}1`$ in general. The factors $`C_q^{\left(i\right)}`$ account for biases introduced by the application of the $`b`$ probability and the jet broadening cuts, as well as for detector effects; these factors were computed by means of the simulation.

A third sample (3) was taken into account by considering the measurement of multiplicity described in . This measurement was performed from a sample of 1297 (3444, 3648) hadronic events, with a contamination of 11% (14%, 18%) after applying all the selection criteria; the remaining background mostly comes from the hadronic decay of W and Z pairs. The values $`n^{\left(3\right)}`$ shown in Table 1 are fully corrected for these backgrounds and for detector effects with their statistical errors; hence the nominal quark flavour ratios appear in the equation (3) below. The systematic errors are reported as the last contribution in Table 3.

The measured mean multiplicities together with the event probability cuts and the factors $`f_q^{\left(i\right)}`$ and $`C_q^{\left(i\right)}`$ are shown in Table 1. For the 1999 data, the values only at $`\sqrt{s}=200`$ GeV are tabulated.

In each of the three samples, the average multiplicity $`n`$ is a linear combination of the unknowns $`n_{b\overline{b}}`$, $`n_{l\overline{l}}`$ and $`n_{c\overline{c}}`$. One can thus formulate a set of three simultaneous equations to compute these unknowns:

$`n^{\left(1\right)}`$ $`=`$ $`f_b^{\left(1\right)}C_b^{\left(1\right)}n_{b\overline{b}}+f_{uds}^{\left(1\right)}C_{uds}^{\left(1\right)}n_{l\overline{l}}+f_c^{\left(1\right)}C_c^{\left(1\right)}n_{c\overline{c}},`$ (1)
$`n^{\left(2\right)}`$ $`=`$ $`f_b^{\left(2\right)}C_b^{\left(2\right)}n_{b\overline{b}}+f_{uds}^{\left(2\right)}C_{uds}^{\left(2\right)}n_{l\overline{l}}+f_c^{\left(2\right)}C_c^{\left(2\right)}n_{c\overline{c}},`$ (2)
$`n^{\left(3\right)}`$ $`=`$ $`f_b^{\left(3\right)}n_{b\overline{b}}+f_{uds}^{\left(3\right)}n_{l\overline{l}}+f_c^{\left(3\right)}n_{c\overline{c}}.`$ (3)
Solving the above equations gave the following mean charge multiplicities at 183 GeV:

$`n_{b\overline{b}}\left(183\mathrm{GeV}\right)`$ $`=29.79\pm 1.11,`$
$`n_{c\overline{c}}\left(183\mathrm{GeV}\right)`$ $`=29.41\pm 4.05,`$
$`n_{l\overline{l}}\left(183\mathrm{GeV}\right)`$ $`=25.25\pm 1.35,`$
$`\delta _{bl}\left(183\mathrm{GeV}\right)`$ $`=4.55\pm 1.31,`$
with correlation coefficient of $`0.45`$ between $`n_{b\overline{b}}`$ and $`n_{l\overline{l}}`$, and at 189 GeV:

$`n_{b\overline{b}}\left(189\mathrm{GeV}\right)`$ $`=30.53\pm 0.70,`$
$`n_{c\overline{c}}\left(189\mathrm{GeV}\right)`$ $`=28.63\pm 2.81,`$
$`n_{l\overline{l}}\left(189\mathrm{GeV}\right)`$ $`=26.10\pm 0.97,`$
$`\delta _{bl}\left(189\mathrm{GeV}\right)`$ $`=4.43\pm 0.85,`$
with correlation coefficient of $`0.52`$ between $`n_{b\overline{b}}`$ and $`n_{l\overline{l}}`$.

From the 1999 data, the results obtained for each energy are tabulated in Table 2. The values were scaled to 200 GeV using JETSET and then a weighted average was calculated using the inverse of the square of the statistical error as weight. One obtains

$`n_{b\overline{b}}\left(200\mathrm{GeV}\right)`$ $`=29.38\pm 0.65,`$
$`n_{c\overline{c}}\left(200\mathrm{GeV}\right)`$ $`=29.89\pm 2.92,`$
$`n_{l\overline{l}}\left(200\mathrm{GeV}\right)`$ $`=25.99\pm 1.03,`$
$`\delta _{bl}\left(200\mathrm{GeV}\right)`$ $`=3.39\pm 0.89,`$
with average correlation coefficient of $`0.52`$ between $`n_{b\overline{b}}`$ and $`n_{l\overline{l}}`$. The difference between the average of the values rescaled to 200 GeV and the average of the values without the scaling was added in quadrature to the final systematic error. This difference is anyway small (0.16 units for $`n_{b\overline{b}}`$ and less than 0.01 units for $`\delta _{bl}`$).

The relatively large uncertainty of the measured mean multiplicities for charm stems from the inability of the $`P_E`$ variable to extract a $`c`$-enriched sample of events.

It should be noted that the transition between particle and detector level measurements in equations (1) and (2) is done via multiplicative factors $`C`$ applied to the mean value of the distributions. The validity of this procedure requires that the simulation used to compute the $`C`$ values reproduces the real data well; the $`chi^2`$/DF at centre of mass energies of 183, 189 and 200 GeV are respectively 0.81, 1.17 and 0.67 for the sample (1) and 0.93, 1.44 and 1.36 for the sample (2).

The analysis was repeated with different cuts applied to the $`b`$-tag probability, $`P_E`$, and the results for the $`\delta _{bl}`$ were found to be quite stable (see Figure 1). A systematic error was evaluated as half of the difference between the greatest and the smallest multiplicity values obtained from varying the cut on $`P_E`$ from $`0.5\times 10^5`$ to $`1.5\times 10^5`$.

The uncertainty due to the event selection in sample (2) was investigated by repeating the analysis after variation of the narrow jet broadening cut, from 0.05 to 0.08. Half of the differences between the greatest and the smallest multiplicities were added in quadrature to the systematic error previously calculated. The propagated systematic error in the total multiplicity in equation (3) from  was also added in quadrature to the systematic error. Uncertainties arising from the modelling of short-lived particles in the simulation were considered. The main physics sources of these uncertainties come from the assumed lifetime of B-hadrons ($`\tau _B=1.564\pm 0.014`$ ps) , and the D$`^+,`$ D0 lifetimes and production rates . Also a variation in the modelling of the $`b`$ fragmentation was investigated, by allowing the average fractional energy of a $`B`$ hadron to vary by 1.5%. The same relative uncertainty was assumed as in .

The effect of a variation of 1% in the fraction $`R_b`$ of $`b\overline{b}`$ events and of 3% in the fraction $`R_c`$ of $`c\overline{c}`$ events was found to be negligible. Since the multiplicity difference, $`\delta _{bl}`$, was found to be independent of energy within errors, the effect of the modelling of the initial state radiation is also expected to be negligible.

The contributions to the systematic error are summarized in Table 3.

The final mean values of the event multiplicity in b events are $`n_{b\overline{b}}\left(183\mathrm{GeV}\right)`$ = $`29.79\pm 1.11\left(stat\right)\pm 0.28\left(syst\right)`$, $`n_{b\overline{b}}\left(189\mathrm{GeV}\right)`$ = $`30.53\pm 0.70\left(stat\right)\pm 0.34\left(syst\right)`$, and $`n_{b\overline{b}}\left(200\mathrm{GeV}\right)`$ = $`29.38\pm 0.65\left(stat\right)\pm 0.50\left(syst\right)`$. The multiplicity difference between $`b\overline{b}`$ and light quark-antiquark events measured at the different energies is:

$`\delta _{bl}\left(183\mathrm{GeV}\right)`$ $`=4.55\pm 1.31\left(stat\right)\pm 0.73\left(syst\right),`$ (4)
$`\delta _{bl}\left(189\mathrm{GeV}\right)`$ $`=4.43\pm 0.85\left(stat\right)\pm 0.61\left(syst\right),`$ (5)
$`\delta _{bl}\left(200\mathrm{GeV}\right)`$ $`=3.39\pm 0.89\left(stat\right)\pm 1.01\left(syst\right).`$ (6)
These values include the products of K$`{}_{S}{}^{}{}_{}{}^{0}`$ and $`\mathrm{\Lambda }`$ decays. The uncertainties on the modelling of the detector largely cancel out in the difference.

Our results on $`\delta _{bl}`$ are plotted in Figure 2 and compared with previous results in the literature.

## 3 Comparison with Models and QCD Predictions

Flavour-Independent Fragmentation — In a model in which the hadronization is independent of the mass of the quarks, one can assume that the non-leading multiplicity in an event, i.e., the light quark multiplicity which accompanies the decay products of the primary hadrons, is governed by the effective energy available to the fragmentation system following the production of the primary hadrons . One can thus write:

$`\delta _{bl}\left(E_{cm}\right)`$ $`=`$ $`2n_B^{\left(decay\right)}+{\displaystyle _0^1}𝑑x_Bf_{E_{cm}}\left(x_B\right){\displaystyle _0^1}𝑑x_{\overline{B}}f_{E_{cm}}\left(x_{\overline{B}}\right)n_{l\overline{l}}\left(\left(1{\displaystyle \frac{x_B+x_{\overline{B}}}{2}}\right)E_{cm}\right)`$ (7)
$``$ $`n_{l\overline{l}}\left(E_{cm}\right),`$
where $`n_B^{\left(decay\right)}`$ is the average number of charged particles coming from the decay of a $`B`$ hadron, $`x_B`$ ($`x_{\overline{B}}`$) is the fraction of the beam energy taken by the $`B`$ ($`\overline{B}`$) hadron, and $`f_{E_{cm}}\left(x_B\right)`$ is the $`b`$ fragmentation function.

We assumed $`2n_B^{\left(decay\right)}=11.0\pm 0.2`$ , consistent with the average $`n_B^{\left(decay\right)}=5.7\pm 0.3`$ measured at LEP . For $`f_{E_{cm}}\left(x_B\right)`$, we assumed a Peterson function with hardness parameter $`ϵ_p=0.0047_{0.0008}^{+0.0010}`$ , evolving with energy as in  to take into account the effects of scaling violations. The value of $`n_{l\overline{l}}\left(E\right)`$ was computed from the fit to a perturbative QCD formula  including the resummation of leading (LLA) and next-to-leading (NLLA) corrections, which reproduces well the measured charged multiplicities , with appropriate corrections to remove the effect of heavy quarks  and leading particles.

The prediction of the model in which the hadronization is independent of the quark mass is plotted in Figure 2. The reason for the drop with collision energy is that the heavy quark system carries away a large fraction of the available energy, approximately (i.e., neglecting scaling violations) linear with $`\sqrt{s}`$, while the multiplicity growth with $`\sqrt{s}`$ is less than linear. There are several variations of this model in the literature, leading to slightly different predictions (see  and references therein). The result from substituting in Eq. (7) $`n_{l\overline{l}}\left(\left(1\frac{x_B+x_{\overline{B}}}{2}\right)E_{cm}\right)`$ with $`n_{l\overline{l}}\left(E_{cm}\sqrt{\left(1x_B\right)\left(1x_{\overline{B}}\right)}\right)`$ as in , or approximating the Peterson fragmentation function with a Dirac delta function at $`x_B`$, are within the errors. Also by using for $`n_{l\overline{l}}`$ the expression in  one stays within the band in Figure 2. The prediction as plotted in Figure 2 agrees with the one calculated in .

QCD Calculation — The large mass of the $`b`$ quark, in comparison to the scale of the strong interaction, $`\mathrm{\Lambda }0.2`$ GeV, results in a natural cut off for the emission of gluon bremsstrahlung. Furthermore, where the c.m. energy greatly exceeds the scale of the $`b`$ quark mass, the inclusive spectrum of heavy quark production is expected to be well described by perturbative QCD in the Modified Leading Logarithmic Approximation (MLLA, ).

The value of $`\delta _{bl}`$ has been calculated in perturbative QCD:

$$\delta _{bl}=2n_B^{\left(decay\right)}n_{l\overline{l}}\left(\sqrt{s}=e^{1/2}m_b\right)+O\left(\alpha _s\left(m_b\right)\right)n_{l\overline{l}}\left(\sqrt{s}=m_b\right).$$
(8)
The reason for the appearance of the $`e^{1/2}`$ factor in the above expression is discussed in detail in . The calculation of the actual value of $`\delta _{bl}`$ in  on the basis of the first two terms in (8) gives a value of $`5.5\pm 0.8`$. A different calculation of $`\delta _{bl}`$ gives 3.68 . These two calculations assume $`m_b=5`$ GeV/$`c^2`$ and $`m_b=4.8`$ GeV/$`c^2`$ respectively, and different parametrizations for the function $`n_{l\overline{l}}\left(\sqrt{s}\right)`$. The dependence of the perturbative part in Eq. (8) on $`m_b`$ is such that moving the $`m_b`$ value from 5 GeV/$`c^2`$ to 4 GeV/$`c^2`$ induces a change of +0.6 units of multiplicity.

The difference of the results in  and in  demonstrates the importance of the contribution proportional to $`\alpha _s\left(m_b\right)`$. A less restrictive condition is the calculation of upper limits: an upper limit $`\delta _{bl}<4.1`$ is given in , based on the maximization of the nonperturbative term; $`\delta _{bl}<4`$ is obtained from phenomenological arguments in Ref. .

Although the presence of the last term in the equation limits the accuracy in the calculation of $`\delta _{bl}`$, QCD tells that $`\delta _{bl}`$ is fairly independent of $`E_{cm}`$. In this article the average of the experimental values of $`\delta _{bl}`$ up to $`m_Z`$ included, $`\delta _{bl}=2.96\pm 0.20`$ (dominated by the LEP 1 data), is taken as the high energy prediction from QCD. The accuracy of the measurement at the Z is thus used to constrain the theoretical prediction.

Our measurement of $`\delta _{bl}`$, as seen in Figure 2, is consistent with the prediction of energy independence based on perturbative QCD, and more than three standard deviations larger than predicted by the naive model presented in the beginning of this section.

## 4 Conclusions

The difference $`\delta _{bl}`$ between the average charged particle multiplicity $`n_{b\overline{b}}`$ in $`e^+e^{}b\overline{b}`$ events and the multiplicity in generic light quark $`l=u,d,s`$ events has been measured at centre-of-mass energies of 183, 189 and 200 GeV:

$`\delta _{bl}\left(183\mathrm{GeV}\right)`$ $`=`$ $`4.55\pm 1.31\left(stat\right)\pm 0.73\left(syst\right)`$
$`\delta _{bl}\left(189\mathrm{GeV}\right)`$ $`=`$ $`4.43\pm 0.85\left(stat\right)\pm 0.61\left(syst\right)`$
$`\delta _{bl}\left(200\mathrm{GeV}\right)`$ $`=`$ $`3.39\pm 0.89\left(stat\right)\pm 1.01\left(syst\right).`$
This difference is in agreement with QCD predictions, while it is inconsistent with calculations assuming that the multiplicity accompanying the decay of a heavy quark is independent of the mass of the quark itself.

### Acknowledgements

We are grateful to Jorge Dias de Deus, Vladimir Petrov, Alexander Kisselev, Valery Khoze and Torbjörn Sjöstrand for useful discussions.
We are greatly indebted to our technical collaborators, to the members of the CERN-SL Division for the excellent performance of the LEP collider, and to the funding agencies for their support in building and operating the DELPHI detector.
We acknowledge in particular the support of
Austrian Federal Ministry of Science and Traffics, GZ 616.364/2-III/2a/98,
FNRS–FWO, Belgium,
FINEP, CNPq, CAPES, FUJB and FAPERJ, Brazil,
Czech Ministry of Industry and Trade, GA CR 202/96/0450 and GA AVCR A1010521,
Danish Natural Research Council,
Commission of the European Communities (DG XII),
Direction des Sciences de la Mati$`\stackrel{`}{\text{e}}`$re, CEA, France,
Bundesministerium f$`\ddot{\text{u}}`$r Bildung, Wissenschaft, Forschung und Technologie, Germany,
General Secretariat for Research and Technology, Greece,
National Science Foundation (NWO) and Foundation for Research on Matter (FOM), The Netherlands,
Norwegian Research Council,
State Committee for Scientific Research, Poland, 2P03B06015, 2P03B1116 and SPUB/P03/178/98,
JNICT–Junta Nacional de Investigação Científica e Tecnol$`\stackrel{´}{\text{o}}`$gica, Portugal,
Vedecka grantova agentura MS SR, Slovakia, Nr. 95/5195/134,
Ministry of Science and Technology of the Republic of Slovenia,
CICYT, Spain, AEN96–1661 and AEN96-1681,
The Swedish Natural Science Research Council,
Particle Physics and Astronomy Research Council, UK,
Department of Energy, USA, DE–FG02–94ER40817.


text:
# Thermostatistics of the G-on Gas

## 1 Derivation of the distribution function of the g-on gas system

### 1.1 Variational method

Let us consider an ideal g-on gas system which consists of N particles. Let the system be in heat and particle baths that evolves towards equilibrium, that is, the total energy E and the total number of g-ons N fluctuate, but, these quantities are conserved on the average. The entropy of the ideal g-on system which is not in equilibrium has been obtained from the statistical weight and is given by ,

$$S_g=\underset{k=1}{\overset{K}{}}g_k[1+n_k(1g)]\mathrm{ln}[1+n_k(1g)(1gn_k)\mathrm{ln}[1gn_k]n_k\mathrm{ln}n_k$$
(1)
where the Boltzmann constant is set to unity, $`g_k`$ is the degeneracy, $`n_k`$ is the occupation number, K is the total number of states; g is the parameter that characterizes the complete or partial action of the exclusion principle i.e. makes interpolation between quantum gases . Following the variational method outlined in the ref., the following equation is obtained

$$[\mathrm{ln}1+n_k(1g)]^{\left(1g\right)}+(1g)+\mathrm{ln}(1gn_k)^g+g\mathrm{ln}n_k1\alpha \beta \epsilon _k=0$$
(2)
where $`\alpha =\frac{\mu }{T}`$

, $`\beta =\frac{1}{T}`$ , $`\epsilon _k`$ is is the energy of the k. state and $`\mu `$ is the chemical potential. After performing some algebra Eq.(2) becomes;

$$\frac{\left[1+n_k\left(1g\right)\right]^{1g}\text{ }\left(1gn_k\right)^g}{n_k}=e^{x_k}$$
(3)
where

$$x_k=\beta \left(\epsilon _k\mu \right).$$
(4)
Eq.(3) is known as Haldane statistics.

Let us write Eq.(3) in a more compact form

$$\frac{1}{y}\left(1y\right)^g=e^x$$
(5)
where y is defined as

$$y=\frac{1}{\frac{1}{n_k}+1g}.$$
(6)
In order to get a real x on the right hand side of Eq.(5), $`y1`$. This inequality imposes a condition

$$gn_k1.$$
(7)
Note that for g=1 inequality(7) gives $`n_k1`$ , for g=0; $`n_{k\text{ }}`$can take any values and for the intermediate values $`n_{k\text{ }}`$is restricted as $`n_k\frac{1}{g}.`$ Therefore, inequality(7) could be interpreted as the generalization of the Pauli Exclusion principle.

Since it is difficult to solve Eq.(5), an approximation is needed for y . Expanding the factor $`(1y)^g`$ in a Taylor series about y =0:

$$\left(1y\right)^g1gy+O\left(y^2\right).$$
(8)
Thus, Eq.(5) becomes

$$n_k=\frac{1}{e^{\beta \left(\epsilon _k\mu \right)}+2g1}.$$
(9)
The last approximation is a linear approximation to Eq.(3) and therefore does not satisfy inequality (7), if g is not close to 0 or 1. Interesting enough, Eq.(9) turns out to be the Maxwell distribution when $`g=\frac{1}{2}`$ .

### 1.2  Steepest Descent Method

In order to get another solution of distribution function $`n_k`$ we are willing to calculate it using the steepest descent optimization method . Two of the authors of the present study FB and DD recently introduced a partition function $`𝒵_g`$ for a grandcanonical ensemble of a g-on gas . Let us write it for the extensive systems as

$$𝒵_g=\underset{k=1}{\overset{K}{}}\left\{\underset{k=0}{\overset{1}{}}e^{x_k}+\left(1g\right)\underset{k=2}{\overset{\mathrm{}}{}}e^{x_k}\right\}$$
(10)
where $`x_k`$ is given by Eq.(4). The given partition function leads to the following partition function for the canonical ensemble

$$Z_g=\underset{k=1}{\overset{K}{}}\left(1+y_k\right)+\left(1g\right)\left(y_k^2+y_k^3+\mathrm{}\right)$$
(11)
where $`y_k=e^{\beta \epsilon _k}`$ and k denoting the quantum state $`\epsilon _k`$. Here g is a parameter that unifies the partition functions of boson and fermion systems under g-on gas system. For g=0 the partition function in Eq.(11) becomes the partition function of a boson system, for g=1 it becomes that of a fermion system. Let us find the distribution function of a g-on system using this partition function and the steepest descent method.

In order to use the steepest descent method one has to express the partition function as a complex integral around a closed contour. For this purpose let us define the function

$$f_g\left(z\right)=\underset{k=1}{\overset{K}{}}\left\{\left(1+zy_k\right)+\left(1g\right)\left(z^2y_k^2+z^3y_k^3+\mathrm{}\right)\right\}.$$
(12)
Using this function, the partition function given in Eq.(11) can be written as

$$Z_g=\frac{1}{2\pi i}\underset{C}{}\frac{f_g\left(z\right)}{z^{N+1}}𝑑z$$
(13)
or

$$Z_g=\frac{1}{2\pi i}\underset{C}{}\underset{k=1}{\overset{K}{}}\left[\left(\frac{1+zy_k}{z^{N+1}}+\frac{\left(1g\right)z^2y_k^2\left(1+zy_k+z^2y_k^2+\mathrm{}\right)}{z^{N+1}}\right)\right]dz.$$
(14)
This integral (or partition function) can be approximated by using the method of steepest descent. In order to see it let us write

$$e^{h\left(z\right)}=\frac{f_g\left(z\right)}{z^{N+1}}$$
(15)
where $`f_g\left(z\right)`$ is given by Eq.(12). Using this definition one can rewrite the $`h(z)`$ as

$$h\left(z\right)=\mathrm{ln}\underset{k=1}{\overset{K}{}}\left[\left(\frac{1+zx_k}{z^{N+1}}+\frac{\left(1g\right)z^2y_k^2\left(1+zy_k+z^2y_k^2+\mathrm{}\right)}{z^{N+1}}\right)\right]$$
(16)
or as

$$h(z)=\underset{k=1}{\overset{\mathrm{}}{}}\left\{\mathrm{ln}\left[\left(1+zy_k\right)+\left(1g\right)z^2y_k^2\left(1+zy_k+z^2y_k^2+\mathrm{}\right)\right]\right\}\mathrm{ln}z^{N+1}.$$
(17)
The function $`h(z)`$ has a minimum at $`z_o`$ and this minimum is very steep because the function $`f_g\left(z\right)`$ consists of infinite number of multiplication and $`N+1`$ is very large. Thus by using the approximation formula of the steepest descent method, one gets an expression for the partition function of a g-on gas as

$$Z_g\left[z_o^{\mathrm{\_}N+1}f_g\left(z_0\right)2\pi f^{\prime \prime }\left(z_0\right)\right]^{1/2}.$$
(18)
In order to find the distribution function one applies the steepest descent condition to h(z), that is

$$h^{}(z_0)=\frac{N+1}{z_0}+\underset{k=0}{\overset{\mathrm{}}{}}\frac{y_k+\left(1g\right)\left(2z_0y_k^2+3z_0^2y_k^3+\mathrm{}\right)}{1+z_0y_k+\left(1g\right)z_0^2y_k^2\left(1+z_0y_k+z_0^2y_k^2+\mathrm{}\right)}=0.$$
(19)
Since $`N1,`$ one can write $`N+1N.`$ Using this approximation and arranging Eq.(18) one obtains

$$N\underset{k=0}{\overset{\mathrm{}}{}}\frac{z_0y_k+\left(1g\right)\left(2z_0^2y_k^2+3z_0^3y_k^3+\mathrm{}\right)}{1+z_0y_k+\left(1g\right)z_0^2y_k^2\left(1+z_0y_k+z_0^2y_k^2+\mathrm{}\right)}.$$
(20)
The factor $`\frac{1}{z_0}`$ can be interpreted as fugacity $`e^{\beta \mu }`$ . For the sake of simplicity let us define a new variable as

$$a_k=e^{\beta \left(\epsilon _k\mu \right)}=e^{x_k}=z_0y_k.$$
(21)
Using the series

$$1+a_k+a_k^2+a_k^3+\mathrm{}=(1a_k)^1$$
(22)
$$1+2a_k+3a_k^2+\mathrm{}=(1a_k)^2$$
(23)
where $`a_k<1`$. By means of these series and the definition (21) one can write Eq.(20) as

$$N=\underset{k=1}{\overset{\mathrm{}}{}}\frac{a_k}{a_k+\frac{ga_k^22a_k+1}{12ga_k+ga_k^2}}=\underset{k=1}{\overset{\mathrm{}}{}}\frac{a_k\left(12ga_k+ga_k^2\right)}{\left(1a_k\right)\left(1ga_k^2\right)}.$$
(24)
or

$$N=\underset{k=1}{\overset{\mathrm{}}{}}\frac{e^{x_k}\left(12ge^{x_k}+ge^{2x_k}\right)}{\left(1e^{x_k}\right)\left(1ge^{2x_k}\right)}.$$
(25)
Since the total number of particles of a system is given as

$$N=\underset{k=1}{\overset{\mathrm{}}{}}\overline{n_k},$$
(26)
the terms in Eq.(25) gives the number of particles which occupy the kth state. Hence for the distribution of g-on particles it is possible to write

$$\overline{n_k}=\frac{e^{x_k}\left(12ge^{x_k}+ge^{2x_k}\right)}{\left(1e^{x_k}\right)\left(1ge^{2x_k}\right)}.$$
(27)
Eq.(27) which unifies the distribution of the quantum gases recovers the well known Bose and Fermi distributions for g=0 and g=1 respectively.

In order to find a more elegant result which could be used moreeasily in the calculations the term in Eq.(24) can be written as

$$\frac{ga_k^22a_k+1}{12ga_k+ga_k^2}1+(2g2)a_{k.}$$
(28)
By taking into account Eq.(28), Eq.(27) reads,

$$\overline{n_k}=\frac{a_k}{1+a_k\left(2g1\right)}$$
(29)
which is valid about g=1 and g=0. Let us write Eq.(29) explicitly;

$$\overline{n_k}=\frac{1}{e^{\beta \left(\epsilon _k\mu \right)}+2g1}.$$
(30)
It is observed that Eq.(30) is same as Eq.(9) which is obtained through an approximation to Haldane formula.

### 1.3  Statistical Method

The partition function of g-ons which is given by Eq.(10) can be used to find the probability of the kth state being occupied by $`n_k`$ g-ons. If the state k with an energy $`\epsilon _k`$ is occupied by $`n_k`$ g-ons then the other states are occupied by N-$`n_k`$ g-ons, where N is the total particle number. Therefore the ratio of the partition function for the states except kth state times $`e^{n_kx_k}`$ to the partition function gives us the probability of finding $`n_{k\text{ }}`$particles in the k. state. In the thermodynamic limit i.e. as $`N\mathrm{}`$ this probability can be written as

$$P_g^k\left(n_k\right)=\frac{e^{n_kx_k}{\displaystyle \underset{\begin{array}{c}i=1\\ ik\end{array}}{}}\left\{\underset{n_i=0}{\overset{1}{}}e^{n_ix_i}+\left(1g\right)\underset{n_i=2}{\overset{\mathrm{}}{}}e^{n_ix_i}\right\}}{{\displaystyle \underset{i=1}{}}\left\{\underset{n_i=0}{\overset{1}{}}e^{n_ix_i}+\left(1g\right)\underset{n_k=2}{\overset{\mathrm{}}{}}e^{n_ix_i}\right\}}.$$
(31)
or writing in a simplified form:

$$P_g^k\left(n_k\right)=\frac{e^{n_kx_k}}{\underset{n_k=0}{\overset{1}{}}e^{n_kx_i}+\left(1g\right)\underset{n_k=2}{\overset{\mathrm{}}{}}e^{n_kx_i}}.$$
(32)
After performing the summation Eq.(32) turns out to be

$$P_g^k\left(n_k\right)=\frac{e^{n_kx_k}\left(1ge^{x_k}\right)}{1ge^{2x_k}}.$$
(33)
It could be mentioned that Eq.(33) is a unified form of the probabilities of occupation number of the quantum gases simply bosons and fermions. One may write the distribution functions by means of Eq.(33) for g-ons as

$$\overline{n_k}=\underset{n_k=0}{\overset{1}{}}n_kP_g^k\left(n_k\right)+\left(1g\right)\underset{n_k=2}{\overset{\mathrm{}}{}}n_kP_g^k\left(n_k\right).$$
(34)
After performing the summations it reads

$$\overline{n_k}=\frac{e^{x_k}\left(1e^{x_k}\right)}{1ge^{2x_k}}\left[g+\frac{1g}{\left(1e^{x_k}\right)^2}\right]=\frac{e^{x_k}(ge^{2x_k}2ge^{x_k}+1)}{\left(1ge^{2x_k}\right)\left(1e^{x_k}\right)}$$
(35)
or more explicitly;

$$\overline{n_k}=\frac{e^{\beta \left(\epsilon _k\mu \right)}\left(12ge^{\beta \left(\epsilon _k\mu \right)}+ge^{2\beta \left(\epsilon _k\mu \right)}\right)}{\left(1e^{\beta \left(\epsilon _k\mu \right)}\right)\left(1ge^{2\beta \left(\epsilon _k\mu \right)}\right)}$$
(36)
which is the same as Eq.(27) as expected.

## 2 Fluctuations in the occupation number

In the preceding section we have calculated the probability $`P_g^k\left(n_k\right)`$ of occupying a g-on of the kth state in the Fock space and the distribution function $`\overline{n_k}`$ .

The deviations of mean values in the distribution function are given by the standard deviations of distributions,

$$\sigma ^2=\overline{n_k^2}\overline{n_k}^2.$$
(37)
Therefore at first $`\overline{n_k^2}`$ has to be calculated. For this purpose let us write the following equation

$$\overline{n_k^2}=\underset{n_k=0}{\overset{1}{}}n_k^2P_g^k\left(n_k\right)+\left(1g\right)\underset{n_k=2}{\overset{\mathrm{}}{}}n_k^2P_g^k\left(n_k\right).$$
(38)
Substituting the expression for $`P_g^k\left(n_k\right)`$ into Eq.(38) and performing the summations one finds

$$\overline{n_k^2}=\frac{e^{x_k}\left(1e^{x_k}\right)}{1ge^{2x_k}}\left[g+\left(1g\right)\frac{1+e^{x_k}}{\left(1e^{x_k}\right)^3}\right].$$
(39)
By means of Eqs.(35),(37) and (38) and one can calculates the variance as

$$\sigma ^2=\frac{e^{x_k}(g^2e^{4x_k}4ge^{3x_k}+6ge^{2x_k}4ge^{x_k}+1)}{\left(1e^{x_k}\right)^2\left(1ge^{2x_k}\right)^2}.$$
(40)
and the ratio $`\frac{\sigma ^2}{\overline{n_k}^2}`$ as

$$\frac{\sigma ^2}{\overline{n_k}^2}=\frac{e^{x_k}(g^2e^{4x_k}4ge^{3x_k}+6ge^{2x_k}4ge^{x_k}+1)}{\left(ge^{2x_k}2ge^{x_k}+1\right)^2}$$
(41)
which reduces to the standard results.

## 3 Thermodynamics of the g-on gas

In this section, as a concrete application of g-on statistics we would like to investigate the thermodynamics of an indistinguishable g-on gas. The ideal g-on gas is a model system, where the influence of quantum statistical effects can be very well studied. Our aim is to calculate the thermodynamical potential of the g-on gas $`\mathrm{\Phi }_g(T,V,\mu )`$.

$$\mathrm{\Phi }_g(T,V,\mu )=k_BT\mathrm{ln}𝒵_g(T,V,z)$$
(42)
where summation approximately is over all of the quantum states. By taking into account equation (10) one derives approximately

$$\mathrm{\Phi }_g(T,V,\mu )=k_BT\underset{k=1}{\overset{K}{}}\mathrm{ln}\left[1+\frac{1}{\omega _k}\right]$$
(43)
where

$$\omega _k\frac{1}{e^{x_k}+g1}.$$
(44)
Since it is difficult to handle with exact formula for $`\overline{n_k}`$ its approximation Eq.(10) or Eq.(30) is appropriate

$$\overline{n_k}=\frac{1}{\omega _k+g}.$$
(45)
Given the generalized thermodynamic potential

$$\mathrm{\Phi }_g(T,V,\mu )=p_gV$$
(46)
the pressure of the g-on gas is obtained as

$$p_gV=k_BT\underset{k=1}{\overset{K}{}}\mathrm{ln}\left[1+\frac{1}{\omega _k}\right]$$
(47)
while the total number of particles is found as

$$N_g=\underset{k}{}\frac{1}{\omega _k+g}.$$
(48)
The state of the gas is determined by Eqs. (47) and (48). But one can also calculate the energy of the gas which is given by

$$E_g=\underset{k}{}\overline{n_k}\epsilon _k=\underset{k}{}\frac{\epsilon _k}{\omega _k+g}.$$
(49)
The logarithm of the grandpartition function $`Q_g`$($`T`$,$`V`$,$`z`$) is given as

$$Q_g(T,V,z)=\mathrm{ln}𝒵_g(T,V,z)=\underset{k=1}{\overset{K}{}}\mathrm{ln}\left[1+\frac{1}{\omega _k}\right].$$
(50)
Let one of the g-on energies $`\epsilon _{k\text{ }}`$are those of the energies of a free quantum mechanical particle in a box of volume V. The chemical potential or fugacity is not fixed, but the total number of g-ons are conserved on the average and z has to be determined from Eq.(48). For a large volume the sum over all single g-on states can be given as the single g-on density states, where spin is not taken into account. Hence,

$$D\left(\epsilon \right)=A\epsilon ^{1/2}$$
(51)
can be written where

$$A=\frac{2\pi V}{h^3}\left(2mh\right)^{3/2}.$$
(52)
Thus the summation in Eq.(50) becomes

$$Q_g(T,V,z)=\underset{0}{\overset{\mathrm{}}{}}\mathrm{ln}\left[1+\frac{1}{\omega _k}\right]\epsilon ^{1/2}𝑑\epsilon $$
(53)
where Eq.(51) is taken into account. When the integration is performed by parts, Q(T,V,z) is reduced to:

$$Q_g(T,V,z)=\frac{2}{3}A\beta \underset{0}{\overset{\mathrm{}}{}}\epsilon ^{3/2}n𝑑\epsilon .$$
(54)
After substituting n from Eq.(45) and then rearranging

$$Q_g(T,V,z)=\frac{2}{3}A\beta z\underset{0}{\overset{\mathrm{}}{}}\frac{\epsilon ^{3/2}}{e^{\beta \epsilon }\left(12g\right)z}𝑑\epsilon $$
(55)
is obtained, where z is the usual fugacity i.e. $`z=e^{\beta \mu }.`$ The integral is calculated using the related formula pp. 326 of reference. Thus one concludes that for a g-on system

$$Q_g(T,V,z)=\frac{2}{3}A\beta zF_{5/2}(g,z)$$
(56)
where the definition

$$F_n(g,z)=\frac{\mathrm{\Gamma }\left(n\right)}{z\left(12g\right)\beta ^n}\underset{k=1}{\overset{K}{}}\frac{\left[z\left(12g\right)\right]^k}{k^n}$$
(57)
is adopted in which $`\mathrm{\Gamma }\left(n\right)`$ is the usual gamma function. Then the state of the g-on system is determined by the following equation:

$$p_gV=k_BTQ_g(T,V,z).$$
(58)
Thus, the equation of state of fermions and bosons is unified

In thermodynamics limit, that is in the limit of infinite volume with particle density held fixed,particles in the ground state, having no contribution of the kinetic energy to the pressure, then, in terms of the thermal wavelength, the pressure of the Bose system becomes, as it is expected,

$$p_0=\frac{kT}{\lambda ^3}\underset{k=1}{\overset{\mathrm{}}{}}\frac{z^k}{k^{5/2}}$$
(59)
where $`\mathrm{\Gamma }\left(\frac{5}{2}\right)=3\frac{\sqrt{\pi }}{4}`$ is taken and $`\lambda `$is the thermal wavelength .

In a similar manner,the thermodynamical properties of a Fermi system follows immediately from the logarithm of the grandpartition function of a g-on system, which is given by Eq. (43). In terms of the thermal wavelength, the pressure of the Fermi system is

$$p_1=g_s\frac{k_BT}{\lambda ^3}\underset{k=1}{\overset{\mathrm{}}{}}\frac{\left(1\right)^{k1}z^k}{k^{5/2}}$$
(60)
where $`g_s`$ is a factor due to spin .

For a large volume, in accordance with equation (49), the total number of g-ons could be written as

$$N_g(T,V,z)=Az\underset{0}{\overset{\mathrm{}}{}}\frac{\epsilon ^{1/2}}{e^{\beta \epsilon }\left(12g\right)z}𝑑\epsilon $$
(61)
assuming that single number of g-on states can be calculated in terms of integrals. Calculation of the integral in Eq.(61) using the related integrals in reference pp326 leads to

$$N_g(T,V,z)=AzF_{3/2}(g,z)$$
(62)
where $`F_{3/2}(g,z)`$ is defined by Eq.(58) . If g=0 and g=1 are substituted in Eq.(62) the formula for the total number of Bose system $`N_0(T,V,z)`$ and Fermi system $`N_1(T,V,z)`$ are respectively recovered .

For a Bose system the total number of particles $`N_0(T,V,z)`$ including $`\epsilon =0`$ energy states, is found to be:

$$N_0(T,V,z)=AzF_{3/2}(z)+\frac{z}{1z}$$
(63)
where the last term represents the contribution of the energy level $`\epsilon =0`$ to the mean particle number.

In a similar manner, for the total number of fermions one gets

$$N_1(T,V,z)=g_s\frac{V}{\lambda ^3}\underset{k=1}{\overset{\mathrm{}}{}}\frac{\left(1\right)^{k1}z^k}{k^{1/2}}$$
(64)
where $`\mathrm{\Gamma }\left(3/2\right)=\sqrt{\pi /2}`$ is taken. This is an expected standard result

Conclusions

This study has been initiated from the point of view bosons and fermions, which obey quantum statistics, could be regarded as g-ons which obey fractional statistics. In this paper, Haldane statistics is rederived by variational method where entropy of the g-on gas is involved. Due to the importance of the distribution function, by means of steepest descent and statistical methods, mean occupation number is derived where partition function of g-on gas is used. It is concluded that mean occupation number which is found by these methods lead to same expressions. Within this context the probabilities of occupation number of g-on gas is also calculated which unifies the probabilities of quantum gases.

The approximate solutions of the distribution function of a g-on gas system, which is found from the variational, steepest descent and statistical methods are developed to an elegant expression that recovers FD distribution for g=1, and BE distribution for g=0. It should be remarked that the partition function which is used in this study is appropriate for the g-on gas since it leads to the same approximate expression for the mean occupation number with Haldane statistics.

As an application of the mean occupation number, the thermodynamical formulation of the g-on gas is derived which gives the opportunity to unify the thermostatistical formulation of quantum gases, simply fermions and bosons.

Authors would like to thank Ege University Research Fund for their partial support under the Project Number 2000 Fen 19.

References

1. Büyükkılıç F., Demirhan D., Eur.Phys.Jour.B. 14 (2000) 704.
2. G.Su, M.Suzuki, Eur Phys.Jour.B. 5 (1998) 577.
3. Y.Su Wu, Phys.Rev.Lett.73 (1994) 922
4. F.D.M. Haldane, Phys.Rev.Lett. 66 (1991) 1529.
5. C.B.Arfken, H.J.Weber, Mathematical Methods for Physicists (Academic Press) 1966
6. C.Kittel, Elementary Statistical Physics (John Wiley, New York) 1958.
7. I.S.Grandshteyn, I.M.Ryzhik, Table of Integral Series and Products (Academic Press, New York) 1965.
8. W.Greiner, L. Neise, H Stoecker, Thermodynamics and Statistical Mechanics (Springer Verlag) 1995.