text:
# Self-Organized Criticality and Universality in a Nonconservative Earthquake Model

## I Introduction

Earthquakes may be the most dramatic example of self-organized criticality (SOC)  that can be observed by humans on earth. Most of the time the crust of the earth is at rest, or quiescent. These periods of stasis are punctuated by sudden, thus far unpredictable bursts, or earthquakes. According to the empirical Gutenberg-Richter (GR) law , the distribution of earthquake events is scale-free over many orders of magnitude in energy. The GR scaling extends from the smallest measurable earthquakes, which are equivalent to a truck passing by, to the most disastrous that have been recorded where hundreds of thousands of people have perished.

The relevance of SOC to earthquakes was first pointed out by Bak and Tang , Sornette and Sornette , as well as Ito and Matsusaki . According to this theory, plate tectonics provides energy input at a slow time scale into a spatially extended, dissipative system that can exhibit breakdown events via a chain reaction process of propagating instabilities in space and time. The GR law arises from the system of driven plates building up to a critical state with avalanches of all sizes. These above-mentioned authors used a spatially extended (but conservative) cellular automata model as a prototype resembling earthquake dynamics which gave a power law distribution of avalanches, or earthquakes. This was followed by a study using block-spring models  by Carlson and Langer , who found characteristic earthquake sizes, rather than asymptotic critical behavior. Later studies of a continuous “train” block-spring model by de Sousa Vieira  recovered criticality. The train model describes a driven elastic chain sliding over a surface with friction. It has been conjectured to be in the same universality class as interface depinning and a model of avalanches in granular piles, which agrees with numerical simulation results .

Several groups made lattice representations of the block-spring model. These models were nonconservative  and were driven uniformly, but did not display SOC. Then Olami, Feder, and Christensen (OFC) introduced a nonconservative model on a lattice that displayed SOC . In this simplified earthquake model, sites on a lattice are continuously loaded with a force. After a threshold force is reached, the sites transfer part of their force to their local neighborhood when discharging. Each discharge event is accompanied by a local loss in accumulated force from the system, when the force on each element is reset to zero. A uniform driving force is slowly applied to all the elements and the model is completely deterministic. This conceptually simple and seemingly numerically tractable model reproduces some of the qualitative phenomenology of the statistics of earthquake events such as power law behavior over a range of sizes, intermittency or clustering of large events , and lack of predictability .

Although SOC and this type of modelling approach has been more or less accepted as a reasonable description of the phenomena of earthquakes (see for example Ref.  and references therein), the OFC model itself has had a controversial existence, both on the theoretical  and numerical side . The initial numerical studies found that the distribution of earthquake sizes obeyed finite size scaling (FSS) over the range of system sizes that could be studied at the time . This placed the nonconservative model into the framework of standard critical behavior. However, these simulations also indicated that there was no universality. In particular the exponents characterizing the power law distributions appeared to vary with both the dissipation parameter, $`\alpha `$, and the form of boundary condition. If this were the case then the OFC model would be very different from familiar critical systems where most microscopic details are irrelevant and have no effect on critical coefficients. In fact an argument was made  that one should not expect universal behavior in far from equilibrium critical phenomena. If this were correct, it would drastically limit the application of any known theoretical tools to these problems.

Another strange aspect was that the dimension $`D`$ characterizing the scaling of the cutoff in the earthquake size distribution was found numerically to be larger than two. This is inconsistent with the fact that each site can only discharge a finite number of times in an earthquake event, requiring $`D2`$ for the two dimensional model . This last result together with the strange lack of universality cast some doubt on whether the OFC model was actually critical or just close to being critical, with some large as yet undetermined length scale beyond which the earthquake distribution would always be cut off. Hwa and Kardar as well as Grinstein and collaborators postulated that conservation of the quantity being transported was required for criticality , but the theoretical arguments made do not take into account SOC phenomena such as avalanches and long-term memory associated with the self-organization process (for more details see ). The fact that the random neighbor version of the nonconservative OFC model is never critical but has an essential singularity as the conservative limit is approached  has added to the mystery surrounding the behavior of the lattice model.

In a previous large scale numerical simulation study of the model discussed here, Grassberger  also claimed that the model was critical but found that some of the conclusions of OFC “have to be modified considerably.” He argued that the probability distribution of earthquake sizes does not show ordinary FSS over the larger range of systems he was able to study, and “conjecture(d) that the cutoff of $`P(s)`$ becomes a step function for $`L\mathrm{}`$,” although he did not present direct numerical evidence of this.

There are a number of important, unresolved questions about the behavior of the model, which have enormous implication for any type of eventual theoretical understanding. Is the nonconservative model on a lattice (or for fixed connectivity matrix) critical? If so, is the critical behavior of the model universal over a range of values of $`\alpha `$, or for different boundary conditions? Is it described by a power-law distribution at all? Are there any other universal quantities? What type of data analysis technique besides FSS would be useful to extract the large scale behavior of the nonconservative model? Our numerical study and analysis will address those issues and answer those questions.

### A Summary

In the first section we review the definition of the OFC model and present some numerical data for the distribution of earthquake sizes using standard FSS. For the range of lattice sizes we have simulated (up to linear size $`L=512`$), our data confirm the lack of apparent FSS in the model, particularly in the cutoff region. In Section III, we present an extensive set of results using a multiscaling method. We analyze the rescaled probability distribution, $`\frac{\mathrm{log}P(s)}{\mathrm{log}L}`$ , in terms of the quantity $`D_{av}\mathrm{log}s/\mathrm{log}L`$, with $`s`$ being the size of an earthquake. We observe that there are no avalanches with $`D_{av}`$ larger than two, consistent with the bound imposed on the cutoff $`s_{co}`$ (see previous discussion). By analyzing how this distribution behaves for different values of the nonconservation parameter, $`\alpha `$, and system size, $`L`$, we show how the multiscaled probability distribution tends to converge to a universal curve as $`L`$ increases. The direction of convergence on increasing $`L`$ changes as $`\alpha `$ varies enabling us to put fairly firm limits on the asymptotic curve. The model appears not be to described at all by FSS. However, for $`s<s_{co}`$ the distribution converges toward a power law with a universal exponent $`\tau 1.8`$ over a range of $`\alpha `$ values. Moreover the cutoff in this distribution becomes very sharp as $`L`$ increases and its behavior indicates that $`s_{co}\mathrm{const}(\alpha )L^2`$ as $`L\mathrm{}`$. In Section IV we summarize our main conclusions.

## II Definition of the Model

We consider a two-dimensional square lattice of $`L\times L`$ sites. At each site $`i`$, a force $`F_i`$ is assigned to be a real variable. Initially, the force at each site is chosen randomly from the uniform distribution between 0 and 1. The dynamics proceeds by two steps in the limit of infinite time scale separation between the slow drive, representing motion of the tectonic plates, and the earthquake process .

1. Increase the force at all sites: Find the largest force $`F_{max}`$ in the system and increase the force at all sites by the same amount $`1F_{max}`$.
2. Relax all unstable sites, i.e. sites with $`F_i1`$: The force of an unstable site is reset to zero, $`F_i0`$, and a fraction of it, $`\alpha F_i`$, is distributed to each of its four nearest neighbours, $`F_{nn}F_{nn}+\alpha F_i`$. This step is repeated in a parallel update until there are no unstable sites left.

This two step rule is iterated indefinitely. The sequence of toppling events (step 2) between application of the uniform drive (step 1) defines an avalanche or earthquake. Since only a fraction, $`4\alpha `$, of the force is redistributed in each toppling, the model is nonconservative for $`\alpha <1/4`$.

To completely define the model we need to specify the boundary conditions, by defining the dissipation parameter, $`\alpha `$, at the sites on the boundaries ($`\alpha _{bc}`$). As in OFC, we consider both “free” and “open” boundary conditions. The sites on the boundaries of the system can be considered to be bounded by fictitious sites with $`F_i=\mathrm{}`$, which can never discharge and only absorb force from the boundary sites. In the case of open boundary conditions, the sites at the boundary have the same $`\alpha `$ parameter as all other sites in the bulk ($`\alpha _{bc}=\alpha `$). In the case of free boundary conditions, the sites at the boundary have their $`\alpha `$ parameters modified in order to have the same level of dissipation per reset event as for sites in the bulk. This latter condition implies $`\alpha _{bc}=\alpha /(1\alpha )`$, except at corner sites where $`\alpha _{bc,c}=\alpha /(12\alpha )`$. The model with periodic boundary conditions is not critical , and will not be discussed. It is probably worthwhile to underline at this point that the model is completely deterministic, the only possible source of randomness coming from the initial conditions.

After a transient period of many earthquakes, the model settles into a statistically stationary state. One way to characterize this state is to measure statistical properties of the earthquakes. The size of an earthquake, $`s`$, is defined as the number of resets of the local force $`F_i0`$ in the system in between applications of the uniform force. One can also measure the temporal duration, $`t`$, in terms of the parallel update, or the radius of gyration $`r`$ of the sites which participated in the earthquake event.

### A Finite Size Scaling

We focus on the probability distribution of earthquake sizes, $`s`$, in a system of size $`L`$, $`P_L(s)`$. If the model is critical then this distribution will have no scale other than the physical extent $`L`$ and the lattice constant, which is set to unity. One ansatz that can describe critical behavior is the FSS ansatz that was previously used by OFC:

$$P_L(s)L^\beta G(\frac{s}{L^D})$$
(1)
where $`G`$ is a suitable scaling function, and $`\beta `$ and $`D`$ are critical exponents describing the scaling of the distribution function. As shown in Fig. 1, the model does not exhibit FSS. In this figure we chose $`D=2`$ as the largest possible allowed value. Still the cutoff in the “collapsed” probability distribution moves to the right as $`L`$ increases. Nevertheless, for earthquake sizes smaller than the cutoff, this figure shows what appears to be a convergence to a well-defined power law, $`P_L(s)s^\tau `$, as $`L`$ increases with the power law exponent $`\tau =\beta /D1.8`$ for both $`\alpha =0.18`$ and $`\alpha =0.21`$, and possibly also for $`\alpha =0.15`$.

## III Multiscaling Analysis

The FSS ansatz is only one possible description of critical behavior. As pointed out some time ago by Kadanoff and co-workers, some SOC phenomena are better described by a multifractal ansatz, rather than FSS . This form has recently been used to clarify the behavior of the Bak-Tang-Weisenfeld  model by Stella and co-workers, who have measured different moments associated with the distribution . For the OFC model, it appears to us that a clearer picture can be obtained by simply examining the probability distribution directly.

The multiscaling ansatz postulates for the probability distribution function $`P_L(s)`$ a form

$$\frac{\mathrm{log}P_L(s/s_o)}{\mathrm{log}(L/l_o)}=F\left(\frac{\mathrm{log}s/s_o}{\mathrm{log}L/l_o}\right),$$
(2)
where $`s_o`$ and $`l_o`$ are parameters which typically (but not always) reflect phenomena at small scales associated with the lattice . Usually, a multiscaling analysis consists of choosing these two parameters in order to get the best collapse for different system sizes. This is quite different than FSS where the critical exponents themselves, reflecting behavior at large scales, must be chosen in an ad hoc manner in order to obtain the “best” collapse. We do not attempt to collapse the data using the multiscaling form of Eq. 2. Instead, we define $`l_o=s_o=1`$ to represent the smallest earthquake which occurs at only one site and involves only one discharge. Moreover, we define the dimension of an earthquake of size $`s`$ in a systems of size $`L`$ as

$$D_{av}=\mathrm{log}s/\mathrm{log}L,$$

and we denote with $`D_{max}`$ the largest value of $`D_{av}`$.

In Fig. 2 we show the multiscaled probability distribution according to Eq. 2, for different system sizes for $`\alpha =0.21`$. One observes immediately that all avalanches have dimension $`D_{av}<2`$, as required, with the largest dimension $`D_{max}`$ approaching two as the system size increases. Also it is clear that the shape of the cutoff function is sharpening as the system size increases. Since the largest dimension $`D_{max}`$ cannot be larger than two, we can infer from this that the cutoff region narrows to the region $`D_{av}2`$ and becomes increasingly sharp.

Note that the increase of $`D_{max}`$ as $`L`$ increases is totally inconsistent with the notion that the OFC model is slightly off criticality, because in that case one would expect that the relative size of the largest dissipating events with respect to the maximum total force allowed in the system would decrease in larger systems. In fact what happens is exactly the opposite. In larger systems a larger fraction of the total energy can be dissipated in the largest event that occurs, and $`D_{max}`$ increases with $`L`$. This result is completely consistent with the nonconservative model being critical, rather than slightly off criticality.

In order to get more explicit visual information on the probability distribution we try to subtract out the leading asymptotic term, which we propose is a power law as suggested by Fig. 2. This is

$$\frac{\mathrm{log}P_L(s)}{\mathrm{log}(L)}=F(D_{av})=(\tau D_{av})+F_{cutoff}(D_{av}),$$
(3)
where $`F_{cutoff}`$ in the limit $`L\mathrm{}`$ should be constant up to a cutoff near $`D_{max}`$.

We get a consistent picture for a range of $`\alpha `$ values by choosing $`\tau =1.8`$, as shown in Fig. 3. Although it appears for small system sizes that smaller values of $`\alpha `$ give a steeper distribution with a larger value of $`\tau `$ (so the lines tend to decrease from left to right rather than remaining horizontal), it is clear from this figure that as the system size increases all the different values of $`\alpha `$ appear to reach the same value of $`\tau =1.8`$, corresponding to a completely horizontal line in this figure. The deviation for small systems is more pronounced for $`\alpha =0.15`$, and less again for $`\alpha =0.18`$, becoming the smallest for $`\alpha =0.21`$. For $`\alpha `$ closer to the conservative limit, as shown for $`\alpha =0.23`$, the approach to the asymptotic horizontal behavior changes direction. Namely smaller systems appear to have a smaller exponent $`\tau `$ than larger systems, at least for small avalanches. Thus rather than having the slope increase as $`L`$ increases, for $`\alpha >.21`$, the apparent slope decreases as $`L`$ increases, as clearly demonstrated in this figure.

We ascribe the change of direction to a crossover effect of the conservative fixed point. For $`\alpha `$ close to 1/4, smaller avalanches behave as avalanches would in the conservative system. It is only the larger avalanches that are affected by nonconservative dissipation. This is associated with the fact that as $`\alpha `$ approaches 1/4, each site can topple more and more times in a single earthquake. For any finite value of dissipation ($`14\alpha `$), the maximum number of times that a site can reset in an earthquake is determined by this dissipation and is finite in the limit of large $`L`$. However for small systems, the largest avalanches are not large enough to be effected by dissipation, and the cutoff in the number of times that a site can topple is not determined by ($`14\alpha `$) but by $`L`$. Then the cutoff in the number of topplings at a given site is the same as in a conservative system of the same size.

The results we have described thus far are for the model with open boundary conditions. Fig. 4 shows that the same behavior occurs for the OFC model with free boundary conditions, with the same value $`\tau =1.8`$. In this case, the asymptotic behavior as $`L`$ increases is approached for decreasing apparent slope for both $`\alpha =0.17`$ and $`\alpha =0.20`$, as is plainly evident in the figure. Again the cutoff appears to sharpen as $`L`$ increases and approach $`D_{max}=2`$.

Our analysis indicates that the OFC model exhibits SOC with a universal power law distribution of events $`P(s)s^{1.8}`$ for $`\mathrm{log}s<(D_{max}\mathrm{log}L)`$. As shown in Fig. 3, the cutoff gets sharper and sharper as $`L`$ increases with $`D_{max}`$ approaching 2 from below as $`L`$ increases for all values of $`\alpha `$ we have studied. As indicated above, the incorrect results obtained by OFC were due to the fact that there is a strong system size dependence that varies with $`\alpha `$ and is not described by FSS. In addition to a systematic change in apparent $`\tau `$ as $`L`$ increases, the largest dimension of earthquakes, $`D_{max}`$, is increasing only slowly towards two, and probability distribution of avalanche dimensions is getting sharper at the cutoff. This means that large avalanches are suppressed in small systems relative to the total amount of force in the system as compared to a larger system, and that a FSS “fit” for any $`L`$ range will always give an apparent $`D>2`$, which is not allowed. In this sense the model appears to violate $`FSS`$ for all values of $`L`$.

## IV Conclusions

The main results of this paper are as follows: The nonconservative model on a two dimensional lattice self-organizes into a critical state. The critical state is robust and universal over a range of values of the dissipation parameter, $`\alpha `$, and for different boundary conditions. The model does not exhibit finite-size scaling. The cutoff becomes sharper as $`L`$ increases with largest earthquakes of dimension $`D_{max}`$ approaching two from below. Nevertheless, the probability distribution of earthquake sizes is a power law with a universal exponent $`\tau 1.8`$. This value can be identified with an exponent $`B=\tau 1`$ for the distribution of energy dissipated in earthquakes. According to the Gutenberg-Richter law this is a power law with $`B=0.8\mathrm{to}1`$, completely consistent with our result.

We thank K. Christensen and H. J. Jensen for helpful conversations, and P. Bak for comments on the manuscript. S.L. was supported by the EPSRC through a post-doctoral research fellowship.


format: text

text:
# References

Today it seems evident that quark potential model gives a rather good description of spin-averaged mass spectrum of hadrons, considered as a system of quarks (see e.g.  and literature cited there). Usually non-relativistic Cornell as well as oscillator potential with Coulumb-like one-gluon exchange or other power-law confinement terms is used. In this work we try to extend this approach to incorporate the second order spin-terms in two-quark Fermi-Breit equation for evaluating the spin-orbit splitting. Instead of calculating as usual spin-terms in first-order perturbation approximation the expansion of the total wave function into a basic set of non-perturbative solutions up to the fifth order of configurationally interacting states is carried out.

The influence of a specific form of potential on a spin-effects (like Cornel, oscillator, screened etc.) which generally give the same results for spin-averaged data, the way of accounting of both long-range non-perturbative and short-range (perturbative) parts in these calculations, the Lorentz structure of different parts of potentials all these questions remain to be unsolved yet.

Our main goal is to clarify some aspects of these questions in framework of configuration interaction (CI) approach . This method does not need the assumption that the coupling constant is to be small, the assumption which is required by perturbation method. Since in quark potential case it is not so the use of perturbation method looks very dubious.

Let us suggest that the static quark potential has vector and scalar property of Lorentz transform:

$$V_{NR}\left(r\right)=V_V\left(r\right)+V_S\left(r\right).$$
(1)
Index ”v” means that potential is a 4-th component of operator $`\widehat{p}_\mu `$, index ”s” means that the potential is scalar.

Following many authors we assume the admixture of vector-scalar screened potential

$$V_V\left(r\right)=V_{OGE}+\epsilon V_{Conf},$$
(2)
$$V_S\left(r\right)=(1\epsilon )V_{Conf},$$
(3)
where

$`V_{OGE}`$ is the one gluon exchange potential,

$`V_{Conf}`$ is the confinement part of the potential.

where $`\epsilon `$ is mixing constant. Here the Lorentz nature of the one-gluon and the confining potential is different, the one-gluon potential being totally vector while the confining potential is a vector-scalar mixture. This choice seems to be reasonable since non-perturbative vertex corrections are important at small $`q^2`$ i.e. large distances and naturally are coupled only with the long-range term in potential. Very interesting rewiew was done in the work of Nora Brambilla conserning the choice of the interaction potential in .

In our case the choice of potential itself is dictated by the consideration of most accurate description of averaged mass spectrum, and here Chikovani-Jenkovszky-Paccanoni (CJP) type potential .

Our approach is based on the model of a non-perturbative gluon propagator, which was recently proposed by Chikovani, Jenkovszky and Paccanoni (CJP) . In this model the propagator has the form

$$D\left(q^2\right)=\frac{c}{\left(q^2\mu ^2\right)}\frac{1}{q^2M^2}$$
(4)
and the corresponding potential assumes the form

$$V\left(r\right)=\frac{g^2}{6\pi \mu }\left(1\mathrm{exp}\left[\mu r\right]\right)\frac{16\pi }{25}\frac{\mathrm{exp}\left[Mr\right]}{r\mathrm{ln}\left(b+1/\left(\mathrm{\Lambda }r\right)^2\right)}.$$
(5)
In fact for accelerating numerical calculations we use simple form $`\frac{\alpha _s}{r}`$ for one-gluon-exchange type term with QCD-like of asymptotic freedom

$$\alpha _s\left(r\right)=\frac{12\pi }{332N}\frac{1}{\mathrm{ln}\left(1/\left(\stackrel{~}{\mathrm{\Lambda }}r\right)^2\right)},$$
(6)
where $`\stackrel{~}{\mathrm{\Lambda }}`$ is taken to be equal to $`\stackrel{~}{\mathrm{\Lambda }}=0.14`$ $`GeV`$. The initial value of $`\alpha _s`$ was defined via calculating $`q\overline{q}`$ masses.

At the finally, we present scalar and vector parts of the potential (2,3) in form

$$V_{OGE}\left(r\right)=\frac{\alpha _S}{r},$$
(7)
$$V_{Conf}\left(r\right)=\frac{g^2}{6\pi }\frac{\left(1e^{\mu r}\right)}{\mu }$$
(8)
Recently screened potential of type (7,8) was successfully applied to a description of spin-averaged meson and baryon mass-spectra . Recently, the study of the quark confining potential from QCD has been well developed. For example, the consideration of the gluon condense may give the confining potential a theoretical description. By considering the loop-diagrams of fermion in the lattice gauge calculation, E. Laermann et al  showed that the quark confinement potential goes lower than the linear confinement potential, when r, the distance between quarks, becomes larger. It reflects the screening between valence quarks and quark sea enhances, consequently, it suppresses the strength of color confinement between valance quarks. According to quantum mechanics, this feature may reduce the energies of highly excited baryons and mesons and help us to understand the experimental data.

Timo Thonhauser  approached the problem of the screened potential from opposite site, namely of Few-Body approach to baryon structure. Few-Body potential for baryons usually has confinement  part plus screened part. Very detailed description of nucleon-nucleon interaction is obtained within this approach. Thonhauser applied this potential to meson structure and obtained quite reasonable values of averaged meson masses. Unfortunately he did not report the results concerning fine splitting. The question of comparison of two approaches needs further study.

In addition, (see e. g. Gerasimov ) the spin-orbit term has to be of short range, as is indicated by quantum chromodinamics (QCD), the condition which is evidently satisfied by CJP potential.

In what follows we shall use the screened potential , which proved to be very good in describing the spin-averaged mass-spectrum of both bosons and baryons as quark systems ,- and which secures the necessary fall-of the spin-dependent forces.

Let us start with two-body Fermi-Breit equation. We shall use nuclear system of units $`\mathrm{}=c=1`$,$`1GeV=\frac{5.068}{1Fm}.`$ The Hamiltonian of the system has the form:

$$\widehat{H}=\widehat{H}_0+\widehat{W},$$
(9)
where

$$\widehat{H}_0=\frac{1}{2m}\mathrm{}+\left(\frac{\alpha _S}{r}+\frac{g^2}{6\pi }\frac{\left(1e^{\mu r}\right)}{\mu }\right),$$
(10)
$`m`$ is the reduced mass,

$$\widehat{W}=\widehat{H}_{LS}+\widehat{H}_{ST},$$
(11)
where spin-dependent potentials are given by :

$``$ Spin-orbit interaction

$$\begin{array}{c}\widehat{H}_{LS}=\frac{1}{4m_1^2m_2^2}\frac{1}{r}\{[((m_1+m_2)^2+2m_1m_2)\stackrel{}{L}\stackrel{}{S_+}+(m_2^2m_1^2)\stackrel{}{L}\stackrel{}{S_{}}]\frac{dV_V}{dr}\\ [(m_1^2+m_2^2)\stackrel{}{L}\stackrel{}{S_+}+(m_2^2m_1^2)\stackrel{}{L}\stackrel{}{S_{}}]\frac{dV_S}{dr}\},\\ \stackrel{}{L}=\stackrel{}{r}\times \stackrel{}{p,}\stackrel{}{S_+}\stackrel{}{S_1}+\stackrel{}{S_{2,}}\stackrel{}{S_{}}\stackrel{}{S_1}\stackrel{}{S_2},\\ \stackrel{}{L}\stackrel{}{S}=\frac{1}{2}\left[j\left(j+1\right)l\left(l+1\right)S\left(S+1\right)\right].\end{array}$$
(12)
Tensor terms

$$\begin{array}{c}\widehat{H}_T=\frac{1}{12m_1m_2}\left[\frac{1}{r}\frac{dV_V}{dr}\frac{d^2V_V}{dr^2}\right]S_{12},\\ S_{12}=\frac{4}{\left(2l+3\right)\left(2l1\right)}\left[\stackrel{}{S}^2\stackrel{}{L}^2\frac{3}{2}\stackrel{}{L}\stackrel{}{S}3\left(\stackrel{}{L}\stackrel{}{S}\right)^2\right].\end{array}$$
(13)
All notations are considered to be familiar and coincide with those which were used by Lucha and Schoberl in their well-known monograph .

Lucha and Shoberl  has indicated that the term $`\frac{1}{r}\frac{dV_V}{dr}`$ causes serious dimensional trouble because of particle falling on the center. But in a more consistent approach based on using Dirac equation we would have instead of previous term the following one $`\frac{1}{EV+mc^2}\frac{1}{r}\frac{dV_V}{dr}`$, which for $`r0`$ will behave like $`\frac{1}{r^2}`$ and the problem is removed. In matrix elements, we have calculated we obtained that, there are no essential differences between these two results. Previously we have studied these questions in .

The ether possibility to solve this problem is the regularisation of the potential . The same idea is based on the using of $`V_{reg}\frac{1}{r+r_0}`$ in the meaning of the regularized potential; where $`r_0`$ is free parameter, which has been obtain from the comparison with the experimental data. So in the case of $`r=0`$ the problem of particle fooling is gone.

Based on this reason we believe that our numerical results are not up set by this problem.

However there are more ”exact” QCD - motivated potentials, which are based on two loop back diagrams . But by more detailed consideration it appears, that these accounts were carried out only for one-gluon exchange part. As races for systems considered by us, the very essential role is played by confinement. Besides use of this potential results in increase of number of parameters, with 1 in our case up to 3. Moreover analysis the obtained results with two loopback potential show that of essential improvement of concurrence theoretical and experimental data is not observed.

In addition to above, indicated terms presented in Fermi-Breit Hamiltonian are the $`\stackrel{}{S_1}\stackrel{}{S_2}`$ (spin-spin) and relativistic correction terms (of $`p^4`$-order). Some authors (like ) indicated that these terms are important for calculating mass spectra, other authors (like ) argue that these results are purely constant. In our case we believe that for LS-mass differences they will not play an essential role.

Now, we consider Fermi-Breit equation

$$\left(\frac{1}{2m}\mathrm{}+\left(\frac{\alpha _S}{r}+\frac{g^2}{6\pi }\frac{\left(1e^{\mu r}\right)}{\mu }\right)+\widehat{W}\right)\mathrm{\Psi }\left(\stackrel{}{r}\right)=\widehat{E}\mathrm{\Psi }\left(\stackrel{}{r}\right).$$
(14)
There are certain difficulties in solving equation (14) in the case of spin-orbit-coupling. Using the potential (11) lead us to the appearance of term of order $`1/r^3`$. Usually authors use one of two possibilities:

a) use the perturbation method;

b) use the numerical computation.

Both of these possibilities are inferior: first one is unacceptable because in some cases the ”fine” splitting turns out to be not so ”fine” at all, being up to 50% contribution to the final mass, the second one needs introducing the cut-off parameter which is highly undesirable.

Here we suggest to use CI approach which was previously very successfully applied in atomic physics . The essence of this approximation is that the total wave function $`\mathrm{\Psi }\left(\stackrel{}{r}\right)`$ is expanded in set of eigen functions $`\phi _n`$of the unperturbated Hamiltonian $`\widehat{H}_0`$, that is

$$\mathrm{\Psi }\left(\stackrel{}{r}\right)=a_n\phi _n\left(\stackrel{}{r}\right).$$
(15)
After substituting (15) into (14) and using eigen value $`E_n^0`$ we obtain the system of linear equations for $`a_n`$ which have to be truncated for reasonably large $`n`$.

$$\begin{array}{c}a_1\left(EE_1^0W_{11}\right)a_2W_{12}a_3W_{13}\mathrm{}a_nW_{1n}=0\\ a_1W_{21}+a_1\left(EE_2^0W_{22}\right)a_3W_{23}\mathrm{}a_nW_{2n}=0\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\\ a_1W_{n1}a_2W_{n2}a_3W_{n3}\mathrm{}+a_n\left(EE_n^0W_{nn}\right)=0\end{array},$$
(16)
where

$$W_{ij}=\phi _i\left|\widehat{W}\right|\phi _j.$$
(17)
Both the basic functions $`\phi _i`$ and matrix elements $`W_{ij}`$ are calculated numerically. Nontrivial solution there will be only if the determinant of this system vanishes can be solved by diagonalizing of the matrix for $`E`$. The system (16) is called CI. This procedure goes far outside of the perturbation method.

The CI approximation turned out to be extremely successful in atomic physics. In calculating atomic structure it allowed to increase the precision of calculating energy levels by one order. In the scattering processes it allowed to reveal fine resonance structure in scattering cross-sections due to formation of auto-ionizing states. So we expect that its applications will be even more important in strong interaction, where the perturbation method is evidently not correct. The technique of application of CIA is quite complicated, since it needs to handle the matrices of large dimensions. In current work we used the code elaborated by O. Zatsarinny .

In this work we have applied the above described method for calculating $`P`$ wave ”fine”-splitting of $`b\overline{b}`$, $`c\overline{c}`$ and $`u\overline{u}`$ systems to $`{}_{}{}^{3}P_{0}^{}`$, $`{}_{}{}^{3}P_{1}^{}`$ and $`{}_{}{}^{3}P_{2}^{}`$ levels. As well as for quarkonium constituted out of different masses like $`s\overline{u}`$, etc. In our case the corresponding operators $`\widehat{H}_{LS}`$, $`\widehat{H}_T`$ will be has the form

$$\widehat{W}_{LS}=\frac{1}{2m^2}\frac{1}{r}\left[3\frac{\alpha _S}{r^2}+\left(4\epsilon 1\right)\frac{g^2}{6\pi }e^{\mu r}\right]\stackrel{}{L}\stackrel{}{S},$$
(18)
$$\widehat{W}_T=\frac{1}{12m^2}\left[3\frac{\alpha _S}{r^3}+\left(\frac{1}{r}\mu \right)\epsilon \frac{g^2}{6\pi }e^{\mu r}\right]S_{12}.$$
(19)
It is important that all parameters except $`\epsilon `$ are taken from ,, where excellent description of bottomonium and charmonium spectra were obtained. Moreover as it was shown in  the same parameters gives good $`\rho `$ meson trajectories. Actually the values $`\frac{g^2}{6\pi }=0.3`$ $`GeV^2`$, $`\mu =0.054`$ $`GeV`$ were taken and $`\alpha _S`$ was taken in accordance to QCD. The only adjustable parameter was $`\epsilon `$. The experimental values are taken from . As mentioned above all calculations were carried out numerically. Special code was constructed for this purpose. The calculations were extended to fifth order in (15) ( see tables, rows 1-8), i.e. until the differences between the results did not go below several $`MeV`$ level.

Table 1. $`b\overline{b}`$-system, $`\alpha _s=0.3`$, $`\epsilon =0.5`$, $`m_b=5.05`$ $`GeV.`$

| State | $`𝚫𝐌_{TH}`$ | order 1 | order 2 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b1}\left(0^+\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`21.32`$ | $`21.23`$ | $`21.3\pm 1.3`$ |
| $`\chi _{b1}\left(0^+\left(1^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`24.25`$ | $`25.51`$ | $`32.1\pm 2`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`45.57`$ | $`46.74`$ | $`53.4\pm 1.9`$ |
| $`R`$ |  | $`0.88`$ | $`0.83`$ | $`0.66`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b1}\left(0^+\left(1^{++}\right)\right)`$ | $`2^3P_22^3P_1`$ | $``$ | $`17.72`$ | $`13.3\pm 0.9`$ |
| $`\chi _{b1}\left(0^+\left(1^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`2^3P_12^3P_0`$ | $``$ | $`19.09`$ | $`23.1\pm 1.1`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`2^3P_22^3P_0`$ | $``$ | $`36.81`$ | $`36.4\pm 1.0`$ |
| $`R`$ |  | $``$ | $`0.93`$ | $`0.57`$ |

| State | $`𝚫𝐌_{TH}`$ | order 3 | order 4 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b1}\left(0^+\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`21.20`$ | $`21.19`$ | $`21.3\pm 1.3`$ |
| $`\chi _{b1}\left(0^+\left(1^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`26.08`$ | $`26.44`$ | $`32.1\pm 2`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`47.29`$ | $`47.63`$ | $`53.4\pm 1.9`$ |
| $`R`$ |  | $`0.81`$ | $`0.8`$ | $`0.66`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b1}\left(0^+\left(1^{++}\right)\right)`$ | $`2^3P_22^3P_1`$ | $`17.61`$ | $`17.58`$ | $`13.3\pm 0.9`$ |
| $`\chi _{b1}\left(0^+\left(1^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`2^3P_12^3P_0`$ | $`20.37`$ | $`21`$ | $`23.1\pm 1.1`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`2^3P_22^3P_0`$ | $`37.98`$ | $`38.58`$ | $`36.4\pm 1.0`$ |
| $`R`$ |  | $`0.86`$ | $`0.84`$ | $`0.57`$ |

Table 2. $`b\overline{b}`$-system, $`\alpha _s=0.3`$, $`\epsilon =0.45`$, $`m_b=5.05`$ $`GeV.`$

| State | $`𝚫𝐌_{TH}`$ | order 1 | order 2 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b1}\left(0^+\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`20.00`$ | $`19.92`$ | $`21.3\pm 1.3`$ |
| $`\chi _{b1}\left(0^+\left(1^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`23.37`$ | $`24.57`$ | $`32.1\pm 2`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`43.37`$ | $`44.49`$ | $`53.4\pm 1.9`$ |
| $`R`$ |  | $`0.85`$ | $`0.81`$ | $`0.66`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b1}\left(0^+\left(1^{++}\right)\right)`$ | $`2^3P_22^3P_1`$ | $``$ | $`16.68`$ | $`13.3\pm 0.9`$ |
| $`\chi _{b1}\left(0^+\left(1^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`2^3P_12^3P_0`$ | $``$ | $`18.49`$ | $`23.1\pm 1.1`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`2^3P_22^3P_0`$ | $``$ | $`35.25`$ | $`36.4\pm 1.0`$ |
| $`R`$ |  | $``$ | $`0.9`$ | $`0.57`$ |

| State | $`𝚫𝐌_{TH}`$ | order 3 | order 4 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b1}\left(0^+\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`19.89`$ | $`19.87`$ | $`21.3\pm 1.3`$ |
| $`\chi _{b1}\left(0^+\left(1^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`25.13`$ | $`25.48`$ | $`32.1\pm 2`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`45.02`$ | $`45.36`$ | $`53.4\pm 1.9`$ |
| $`R`$ |  | $`0.79`$ | $`0.78`$ | $`0.66`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b1}\left(0^+\left(1^{++}\right)\right)`$ | $`2^3P_22^3P_1`$ | $`16.76`$ | $`16.65`$ | $`13.3\pm 0.9`$ |
| $`\chi _{b1}\left(0^+\left(1^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`2^3P_12^3P_0`$ | $`19.73`$ | $`20.34`$ | $`23.1\pm 1.1`$ |
| $`\chi _{b2}\left(0^+\left(2^{++}\right)\right)\chi _{b0}\left(0^+\left(0^{++}\right)\right)`$ | $`2^3P_22^3P_0`$ | $`36.38`$ | $`36.95`$ | $`36.4\pm 1.0`$ |
| $`R`$ |  | $`0.85`$ | $`0.82`$ | $`0.57`$ |

Table 3. $`c\overline{c}`$-system, $`\alpha _s=0.386`$, $`\epsilon =0.3`$, $`m_c=1.675`$ $`GeV.`$

| State | $`𝚫𝐌_{TH}`$ | order 1 | order 2 | order 3 | order 4 |
| --- | --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\chi _{c2}\left(0^+\left(2^{++}\right)\right)\chi _{c1}\left(0^+\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`52.20`$ | $`51.66`$ | $`51.49`$ | $`51.41`$ |
| $`\chi _{c1}\left(0^+\left(1^{++}\right)\right)\chi _{c0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`69.31`$ | $`77.27`$ | $`81.55`$ | $`84.51`$ |
| $`\chi _{c2}\left(0^+\left(2^{++}\right)\right)\chi _{c2}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`121.51`$ | $`128.93`$ | $`133.04`$ | $`135.93`$ |
| $`R`$ |  | $`0.75`$ | $`0.66`$ | $`0.63`$ | $`0.61`$ |

| State | $`𝚫𝐌_{TH}`$ | order 5 | order 6 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\chi _{c2}\left(0^+\left(2^{++}\right)\right)\chi _{c1}\left(0^+\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`51.37`$ | $`51.33`$ | $`45.64\pm 0.25`$ |
| $`\chi _{c1}\left(0^+\left(1^{++}\right)\right)\chi _{c0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`86.72`$ | $`88.48`$ | $`95.43\pm 1.12`$ |
| $`\chi _{c2}\left(0^+\left(2^{++}\right)\right)\chi _{c2}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`138.09`$ | $`139.82`$ | $`141.07\pm 1.13`$ |
| $`R`$ |  | $`0.59`$ | $`0.58`$ | $`0.48`$ |

Table 4. $`c\overline{c}`$-system, $`\alpha _s=0.386`$, $`\epsilon =0.4`$, $`m_c=1.675`$ $`GeV.`$

| State | $`𝚫𝐌_{TH}`$ | order 1 | order 2 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\chi _{c2}\left(0^+\left(2^{++}\right)\right)\chi _{c1}\left(0^+\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`66.81`$ | $`66.17`$ | $`45.64\pm 0.25`$ |
| $`\chi _{c1}\left(0^+\left(1^{++}\right)\right)\chi _{c0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`79.36`$ | $`88.57`$ | $`95.43\pm 1.12`$ |
| $`\chi _{c2}\left(0^+\left(2^{++}\right)\right)\chi _{c2}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`146.17`$ | $`154.74`$ | $`141.07\pm 1.13`$ |
| $`R`$ |  | $`0.84`$ | $`0.75`$ | $`0.48`$ |

| State | $`𝚫𝐌_{TH}`$ | order 3 | order 4 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\chi _{c2}\left(0^+\left(2^{++}\right)\right)\chi _{c1}\left(0^+\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`65.99`$ | $`65.92`$ | $`45.64\pm 0.25`$ |
| $`\chi _{c1}\left(0^+\left(1^{++}\right)\right)\chi _{c0}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`93.45`$ | $`96.75`$ | $`95.43\pm 1.12`$ |
| $`\chi _{c2}\left(0^+\left(2^{++}\right)\right)\chi _{c2}\left(0^+\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`159.45`$ | $`162.67`$ | $`141.07\pm 1.13`$ |
| $`R`$ |  | $`0.70`$ | $`0.68`$ | $`0.48`$ |

Table 5. $`u\overline{u}`$-system, $`\alpha _s=0.52`$, $`\epsilon =0.14`$, $`m_u=0.33`$ $`GeV.`$

| State | $`𝚫𝐌_{TH}`$ | order 1 | order 2 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`a_2\left(1^{}\left(2^{++}\right)\right)a_1\left(1^{}\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`15.34`$ | $`10.26`$ | $`88.2\pm 40.7`$ |
| $`a_1\left(1^{}\left(2^{++}\right)\right)a_0\left(1^{}\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`241.19`$ | $`316.91`$ | $`246.5\pm 40.9`$ |
| $`a_2\left(1^{}\left(2^{++}\right)\right)a_0\left(1^{}\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`256.53`$ | $`327.17`$ | $`334.6\pm 1.6`$ |
| $`R`$ |  | $`0.063`$ | $`0.032`$ | $`0.357`$ |

| State | $`𝚫𝐌_{TH}`$ | order 3 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`a_2\left(1^{}\left(2^{++}\right)\right)a_1\left(1^{}\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`8.57`$ | $`88.2\pm 40.7`$ |
| $`a_1\left(1^{}\left(2^{++}\right)\right)a_0\left(1^{}\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`374.09`$ | $`246.5\pm 40.9`$ |
| $`a_2\left(1^{}\left(2^{++}\right)\right)a_0\left(1^{}\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`382.66`$ | $`334.6\pm 1.6`$ |
| $`R`$ |  | $`0.023`$ | $`0.357`$ |

Table 6. $`u\overline{u}`$-system, $`\alpha _s=0.52`$, $`\epsilon =0.145`$, $`m_u=0.33`$ $`GeV.`$

| State | $`𝚫𝐌_{TH}`$ | order 1 | order 2 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`a_2\left(1^{}\left(2^{++}\right)\right)a_1\left(1^{}\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`24.41`$ | $`19.15`$ | $`88.2\pm 40.7`$ |
| $`a_1\left(1^{}\left(2^{++}\right)\right)a_0\left(1^{}\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`247.52`$ | $`324.96`$ | $`246.5\pm 40.9`$ |
| $`a_2\left(1^{}\left(2^{++}\right)\right)a_0\left(1^{}\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`271.92`$ | $`344.10`$ | $`334.6\pm 1.6`$ |
| $`R`$ |  | $`0.098`$ | $`0.059`$ | $`0.357`$ |

| State | $`𝚫𝐌_{TH}`$ | order 3 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`a_2\left(1^{}\left(2^{++}\right)\right)a_1\left(1^{}\left(1^{++}\right)\right)`$ | $`1^3P_21^3P_1`$ | $`17.45`$ | $`88.2\pm 40.7`$ |
| $`a_1\left(1^{}\left(2^{++}\right)\right)a_0\left(1^{}\left(0^{++}\right)\right)`$ | $`1^3P_11^3P_0`$ | $`383.24`$ | $`246.5\pm 40.9`$ |
| $`a_2\left(1^{}\left(2^{++}\right)\right)a_0\left(1^{}\left(0^{++}\right)\right)`$ | $`1^3P_21^3P_0`$ | $`400.69`$ | $`334.6\pm 1.6`$ |
| $`R`$ |  | $`0.046`$ | $`0.357`$ |

Table 7. $`s\overline{u}`$-system, $`\alpha _s=0.421`$, $`\epsilon =0.1875`$, $`m_s=0.5`$ $`GeV,m_u=0.33`$ $`GeV.`$

| State | $`𝚫𝐌_{TH}`$ | order 1 | order 2 | order 3 | order 4 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\mathrm{?}`$ | $`1^3P_21^3P_1`$ | $`43.16`$ | $`40.23`$ | $`39.29`$ | $`38.87`$ | $``$ |
| $`\mathrm{?}`$ | $`1^3P_11^3P_0`$ | $`187.78`$ | $`230.62`$ | $`259.86`$ | $`282.63`$ | $``$ |
| $`\mathrm{?}`$ | $`1^3P_21^3P_0`$ | $`230.94`$ | $`270.85`$ | $`299.15`$ | $`321.50`$ | $``$ |
| $`R`$ |  | $`0.23`$ | $`0.17`$ | $`0.15`$ | $`0.14`$ | $``$ |

Table 8. $`s\overline{u}`$-system, $`\alpha _s=0.421`$, $`\epsilon =0.2`$, $`m_s=0.5`$ $`GeV,m_u=0.33`$ $`GeV.`$

| State | $`𝚫𝐌_{TH}`$ | order 1 | order 2 | order 3 | $`𝚫𝐌_{EXP}[18]`$ |
| --- | --- | --- | --- | --- | --- |
|  |  | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ | $`\mathrm{𝐌𝐞𝐕}`$ |
| $`\mathrm{?}`$ | $`1^3P_21^3P_1`$ | $`60.13`$ | $`56.95`$ | $`56.00`$ | $``$ |
| $`\mathrm{?}`$ | $`1^3P_11^3P_0`$ | $`199.43`$ | $`244.92`$ | $`331.74`$ | $``$ |
| $`\mathrm{?}`$ | $`1^3P_21^3P_0`$ | $`259.56`$ | $`301.84`$ | $`387.75`$ | $``$ |
| $`R`$ |  | $`0.30`$ | $`0.23`$ | $`0.17`$ | $``$ |

Let us make the following conclusions:

1. The results for heavy quarkonium are quite good for values $`\epsilon =0.30.45`$. Which is coincide with results obtained by Lai-Him Chan  For light quarkonium the results are worse, which means that relativistic effects have to be taken into account more carefully. We believe, that the difference between $`\epsilon `$ in case of $`b\overline{b}`$ and $`u\overline{u}`$ system exactly reflects this fact.

2. The value of indicates that confinement has prevailingly scalar character. This conclusion does not contradict other authors .

3. As it follows from (18) at $`\epsilon =0.25`$ $`\left(4\epsilon 1\right)=0`$ the contribution of confinement vanishes totally. May be exactly this circumstances was the reason that some authors stated the pure one-gluon character of LS-splitting.

4. The first column in tables corresponds to pure perturbation approach. It is clearly seen that this approach gives only rough qualitative estimate, but the results are drastically improving with switch on the CIA expansion. We believe that the using of CIA approach in quark physics will has bright future. Very convenient for fine-splitting characteristic is the coefficient

$$R=\frac{M\left({}_{}{}^{3}P_{2}^{}\right)M\left({}_{}{}^{3}P_{1}^{}\right)}{M\left({}_{}{}^{3}P_{1}^{}\right)M\left({}_{}{}^{3}P_{0}^{}\right)}.$$
(20)
As cited in  the experimental values of this parameter are $`R=0.66`$ for $`b\overline{b}`$, $`R=0.48`$ for $`c\overline{c}`$ and also $`R=0.21`$ for $`u\overline{u}`$. It is interesting to note that firstly this parameter was introduced in atomic physics  where its value was established for atoms with two outer electrons (like Mg) as being $`R=2`$. Since the fine mass splitting in atoms is totally due to one-photon exchange this value gives good indication as to the nature of the Lorentz character of $`q\overline{q}`$ potential. As our calculations show the one-gluon term is totally vector while non-perturbative confinement term is of about 75-80% scalar.

Confinement gives 20-25% plus to the results, obtained with using only the one gluon exchange term, for the description of the $`LS`$, $`ST`$ interaction. Finally we obtain description of fine-effects at the quark systems with the precision 80-90%.

These conclusions confirms with the qualitative estimation of Lucha and Shöberl . This result approximately confirms with the conclusion of Franzini  in the part, that the confinement is pure scalar. We show that the results are very sensitive to the exact value of $`\epsilon `$ and small change of $`\epsilon `$ can destroy the agreement with the experiment what actually occurred in the icase of Franzini.

We believe that exactly this feature will be successfully exploited in near future not only in hadronic structure, but also explaining resonances in $`\pi \pi `$ or $`\pi N`$ resonances.


text:
# Flow equation analysis of the anisotropic Kondo model

## I Introduction

The Kondo model, originally introduced by Zener , is the most basic Hamiltonian describing the physics of dilute magnetic impurities in metals. Despite its apparent simplicity, it has become a paradigm for complex many–body effects, due to its non–trivial strong–coupling behaviour at low temperatures. Unprecedented theoretical activity during the last four decades was triggered by the seminal work of Kondo , who showed that perturbation theory for an antiferromagnetic exchange coupling diverges at low temperatures. Within perturbative scaling  the coupling constant grows continuously and eventually diverges. This indicates that the impurity spin is screened in the ground state. However, the diverging coupling constant also implies the breakdown of the perturbative scaling approach in the strong–coupling phase. Therefore perturbative scaling does not provide a systematic and controlled expansion describing the weak to strong–coupling crossover. This is a typical situation in many *strong–coupling* problems. The development of new, nonperturbative methods like Wilson’s numerical renormalization group (NRG)  and the Bethe ansatz  eventually led to a quantitative understanding of the strong–coupling regime in the Kondo model. Despite these successes, however, a simple *analytical* description of the crossover from a free spin to the Fermi liquid regime in an RG type framework was still missing.

In the present work, we will attempt to fill this gap by applying Wegner’s *flow equation method* , which was recently employed successfully to diagonalize a related strong–coupling problem, the quantum sine-Gordon model . We will mainly focus on the calculation of equilibrium dynamical correlation functions at zero temperature as this provides one of the most interesting new perspectives for this new approach. The calculation of correlation functions throughout the crossover region is notoriously difficult with exact methods building on the integrability of the model. There is considerable interest in theoretical tools that allow their determination in strong–coupling models, which provides the main motivation for our work . Also the integrability of the Kondo model is not used in our method and nonintegrable perturbations can be studied as well.

We consider the *anisotropic* Kondo Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{k}{}}ϵ_k^{}c_{k\alpha }^{}c_{k\alpha }^{}+J_{}c_{0\alpha }^{}\sigma _{\alpha \beta }^zc_{0\beta }^{}S^z`$ (2)
$`+{\displaystyle \frac{J_{}}{2}}(c_{0\alpha }^{}\sigma _{\alpha \beta }^+c_{0\beta }^{}S^{}+h.c.).`$
Here $`\stackrel{}{S}`$ is a spin-1/2 degree of freedom and $`c_{0\alpha }^{()}=_kc_k^{()}/\sqrt{L}`$ are creation and annihilation operators for the localized electron state at the spin site. Notice that the anisotropic Kondo Hamiltonian is equivalent to the dissipative two-state system (spin boson model) . For the noninteracting conduction band we assume a linear dispersion. We will treat the model in its bosonized form : One introduces the bosonic spin density modes $`\sigma (p)=\frac{1}{\sqrt{2|p|}}_q\left(c_{p+q}^{}c_q^{}c_{p+q}^{}c_q^{}\right)`$ with the commutator $`[\sigma (q),\sigma (q^{})]=\delta _{qq^{}}L/2\pi `$ for $`q,q^{}>0`$. Here $`L`$ is the system size and all other nonvanishing commutators can be derived using $`\sigma ^{}(q)=\sigma (q)`$. Notice $`\sigma (q)|\mathrm{\Omega }>=0`$ on the ground state for $`q>0`$. The charge density modes in (2) decouple completely and one only has to look at the spin density part of (2)

$$H=H_0\frac{J_{}}{\sqrt{2}2\pi }_x\mathrm{\Phi }(0)S^z+\frac{J_{}}{4\pi a}(e^{i\sqrt{2}\mathrm{\Phi }(0)}S^{}+h.c.).$$

Here $`H_0=\frac{2\pi }{L}v_F_{q>0}q\sigma (q)\sigma (q)`$. In the following we set $`v_F=1`$. The bosonic field is defined as

$$\mathrm{\Phi }(x)=\frac{2\pi i}{L}\underset{q0}{}\frac{\sqrt{|q|}}{q}e^{iqxa|q|/2}\sigma (q).$$

The parameter $`a>0`$ generates the UV–regularization of our model. We are interested in its universal properties at energies $`|E|a^1`$. The longitudinal spin coupling can be eliminated by a unitary transformation $`U=\mathrm{exp}\left(i\mu S^z\mathrm{\Phi }(0)\right)`$ with $`\mu =J_{}/\sqrt{2}2\pi `$. We arrive at

$$H=H_0+g_0\left(V(\lambda _0,0)S^{}+V(\lambda _0,0)S^+\right)$$
(3)
with the coupling constant $`g_0=J_{}/4\pi a`$ and $`\lambda _0=\sqrt{2}J_{}/\sqrt{2}2\pi `$. Here *vertex operators* have been introduced

$$V(\lambda ,x)\mathrm{exp}\left(i\lambda \mathrm{\Phi }(x)\right),$$

which are complicated many-body objects for a general value of the *scaling dimension* $`\lambda `$ . Notice that our $`V(\lambda ,x)`$ are not normal–ordered. We will also need the Fourier-transformed normalized vertex operators $`C_p^{}=\alpha _p^1L^{1/2}𝑑xe^{ipx}V(\lambda ,x)`$ with $`\alpha _p^2=2\pi a|pa|^{\lambda ^21}/\mathrm{\Gamma }(\lambda ^2)`$. Notice $`C_p^{}|\mathrm{\Omega }=C_p^{}|\mathrm{\Omega }=0`$ for $`p>0`$, and $`\mathrm{\Omega }|C_p^{}C_p^{}^{}|\mathrm{\Omega }=\delta _{pp^{}}\theta (p)`$, $`\mathrm{\Omega }|C_p^{}C_p^{}^{}|\mathrm{\Omega }=\delta _{pp^{}}\theta (p)`$ for $`|ap|1`$.

Eq. (3) with general initial couplings $`\lambda _0`$ and $`g_0`$ is our starting point for the flow equation approach. The universal low–energy properties of (3) are equivalent to those of the spin–density part of (2). As proposed by Wegner , we apply a sequence of infinitesimal unitary transformations to (3) in order to make it successively more diagonal. One can set this up in a differential formulation

$$_BH(B)=[\eta (B),H(B)]$$
(4)
with an anti–Hermitian generator $`\eta (B)`$. $`H(B=0)`$ is our initial Hamiltonian (3) and all $`H(B)`$ as obtained by the solution of (4) are unitarily equivalent to it. In general we will have to truncate our system of equations, so this equivalence will only be approximate. In order to generate a stable sequence of approximations, we choose $`\eta (B)`$ such that first (for small $`B`$) offdiagonal matrix elements are removed which couple states with large energy differences (of order $`B^{1/2}`$), and later more degenerate states. This is reminiscent of the energy scale separation underlying perturbative scaling. Similar ideas have independently been developed by Głazek and Wilson (*similarity renormalization scheme*)  that contain Wegner’s method as a special case.

As the flow parameter $`B`$ grows, the effective spin flip interaction in the interaction part $`H_{\mathrm{int}}(B)`$ will become increasingly nonlocal and the scaling dimension of the vertex operator will begin to flow. We take this into account by writing

$$H(B)=H_0+dxg(B,x)(V(\lambda (B),x)S^{}+\mathrm{h}.\mathrm{c}.)$$
(5)
with $`g(B,x)=L^{1/2}_pg_p(B)e^{ipx}`$. Wegner has suggested the generic choice $`\eta ^{(1)}=[H_0,H_{\mathrm{int}}]`$ for the generator, as this always removes the off–diagonal matrix elements in the above energy–scale separated way . Our generator has the following more general form

$`\eta `$ $`=`$ $`{\displaystyle \underset{p}{}}\eta _p^{(1)}\left(C_p^{}S^{}C_p^{}S^+\right)`$ (7)
$`+{\displaystyle \underset{p,q}{}}\eta _{pq}^{(2)}\left(C_p^{}C_q^{}C_q^{}C_p^{}\right).`$
Here $`\eta _p^{(1)}=p\alpha _pg_p`$ follows from Wegner’s canonical choice. The second generator part has been introduced since it will later allow us to generate transformed Hamiltonians with a particularly simple structure. The precise form of $`\eta _{pq}^{(2)}`$ will be discussed below.

In general, new many-particle interactions are generated in (4). Our truncation criterion is like in  the *operator product expansion* (OPE) for vertex operators. For example in the anticommutator of vertex operators the OPE leads to

$`\{V(\lambda ,x),V(\lambda ,y)\}=\left(1+i\lambda (xy)_x\mathrm{\Phi }(x)+\mathrm{}\right)`$ (9)
$`\times \left({\displaystyle \frac{1}{[1+i(xy)/a]^{\lambda ^2}}}+{\displaystyle \frac{1}{[1i(xy)/a]^{\lambda ^2}}}\right).`$
We only keep the leading nonvanishing term of this expansion in our calculation. This amounts to neglecting terms with larger scaling dimensions (more irrelevant terms) and is therefore reminiscent of an RG approach. Higher order terms can be successively taken into account in a systematic expansion. However, the present approximation will already turn out to lead to very accurate results. Within our approximation, the primary effect of the flow (4) is to generate a new interaction term

$$H_{\mathrm{new}}^{(1)}=S^z𝑑xf(x)_x\mathrm{\Phi }(x)$$
(10)
in (5) with a function $`f(x)`$ that depends on the parameters $`g_p`$. (10) is identical to the initial longitudinal spin coupling. Our strategy is to eliminate $`H_{\mathrm{new}}^{(1)}`$ after it is generated infinitesimally by again performing a unitary transformation $`U=\mathrm{exp}[iϵS^z\mathrm{\Phi }(0)]`$ with a suitable infinitesimal $`ϵ`$. Similar to the initial value of the scaling dimension that is determined by the longitudinal Kondo coupling, we find a further shift in $`\lambda (B)`$ due to the elimination of $`H_{\mathrm{new}}^{(1)}`$: one derives the following flow equation

$$_B\lambda ^2=\frac{8\pi a\lambda ^2(1\lambda ^2)}{\mathrm{\Gamma }(\lambda ^2)}\underset{p}{}g_pg_p|pa|^{\lambda ^21}.$$
(11)
In the sequel we will only look at the strong–coupling phase ($`|J_{}|>J_{}`$ for small couplings $`|J_{}|1`$) since the weak coupling regime is trivial. In the strong–coupling phase the flow is always directed towards the *Toulouse line* $`\lambda =1`$ , which constitutes the *strong–coupling fixed point* of our approach (see fig. 1). For $`\lambda =1`$ the vertex operators in (3) obey fermionic anticommutation relations  and the Hamiltonian is quadratic (equivalent to a noninteracting resonant level model ): The OPE series (9) terminates after the leading term and the flow equations close *exactly* since they can readily diagonalize a quadratic Hamiltonian. We flow to a point where our method becomes exact, thereby avoiding the usual strong–coupling divergence. This is the *key improvement* of our approach as compared to perturbative scaling, which does not make use of the special properties at $`\lambda =1`$.

The second new interaction that is generated is a potential scattering term. Only diagonal terms appear

$$H_{\mathrm{new}}^{(2)}(B)=\underset{p}{}\omega _p(B)\left(C_p^{}C_p^{}C_p^{}C_p^{}\right)$$

with the following choice of the generator $`\eta _{pq}^{(2)}(B)`$

$$\eta _{pq}^{(2)}=\frac{1}{2}\frac{(p+q)\alpha _p\alpha _qg_pg_q}{pq+2(\omega _p\omega _q)}.$$

$`H_{\mathrm{new}}^{(2)}(B=\mathrm{})`$ is a term in order $`1/L`$ since $`\omega _p(B)=O(1/L)`$ and thus represents the only ”impurity effect” in the Hamiltonian after elimination of the Kondo coupling for $`B\mathrm{}`$ in our present approximation

$$H(B=\mathrm{})=H_0+H_{\mathrm{new}}^{(2)}(B=\mathrm{}).$$

Together with (11) we thus obtain a closed set of flow equations

$`_Bg_p`$ $`=`$ $`p^2g_p+2{\displaystyle \underset{q}{}}{\displaystyle \frac{\alpha _q}{\alpha _p}}\eta _{pq}^{(2)}g_q`$ (13)
$`+{\displaystyle \frac{1}{2}}g_p\mathrm{ln}\left({\displaystyle \frac{B}{a^2}}\right)\lambda _B\lambda 2pg_p\omega _p`$
$`_B\omega _p`$ $`=`$ $`p\alpha _p^2g_p^2.`$ (14)
Neglecting the momentum dependence induced by $`\eta ^{(2)}`$ and the potential scattering terms $`\omega _p`$ (which become important only at energies of the order of the Kondo scale), we can parametrize the solution of (13) as $`g_p(B)=\stackrel{~}{g}(B)e^{p^2B}`$ with the *running coupling* $`\stackrel{~}{g}(B)`$. $`\stackrel{~}{g}(B)`$ is in general strongly renormalized. In the low energy limit the scaling dimension approaches $`\lambda =1`$ and we can then deduce the Kondo scale from the effective Toulouse Hamiltonian where $`T_K\stackrel{~}{g}(B=\mathrm{})^2/a`$. For the symmetric case $`J_{}=J_{}=J`$ and a small antiferromagnetic coupling $`0<J1`$ one e.g. shows $`T_K\frac{J^\tau }{a}\mathrm{exp}\left(\frac{2\pi }{J}\right),`$ where $`\tau =(1+\gamma \mathrm{ln}2)/30.295`$. NRG and third order scaling results by Wilson  predict the same leading exponential behaviour but a slightly different prefactor with $`\tau =1/2`$.

Our main focus in this work is the calculation of the dynamical impurity susceptibility $`\chi (t)=i\theta (t)[S^z(t),S^z(0)]`$ at zero temperature. The key ideas for the calculation of correlation functions for impurity systems within the flow equation approach have been developed in : we need to evaluate matrix elements of the observable $`O=S^z`$ with respect to the eigenstates of $`H(B=\mathrm{})`$, *which is simplified by the fact that the final Hamiltonian is diagonal.* In order to make use of this simplification, we need to perform the same series of infinitesimal unitary transformations on $`O`$ as on $`H(B)`$

$$_BO(B)=[\eta (B),O(B)].$$
(15)
Usually, keeping the Hamiltonian flow (4) simple leads to a complicated structure of the transformed observable. In the present case, however, it turns out that the OPE can be used as a unifying truncation criterion for *both* $`H(B)`$ and $`O(B)`$. Using $`S^z=[S^+,S^{}]/2`$ and the ansatz

$$S^+(B)=h(B)S^++S^z\underset{p}{}d_p(B)C_p^{}$$

for the flow of the observable  another closed set of equations is obtained:

$`_Bh`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{p}{}}p\alpha _pd_pg_p`$ (16)
$`_Bd_p`$ $`=`$ $`2pg_p\alpha _ph+2{\displaystyle \underset{q}{}}\eta _{pq}^{(2)}d_q.`$ (17)
We have numerically solved this system of equations in combination with (13). It is known from the theory of dissipative quantum systems  that in the thermodynamic limit the impurity observable $`S^+`$ decays completely, i.e. $`h(B=\mathrm{})=0`$. With $`S^z(B=\mathrm{})`$ completely expressed in terms of band operators, it is then trivial to calculate the dynamical susceptibility: since $`H_{\mathrm{new}}^{(2)}(B=\mathrm{})`$ is of order $`1/L`$, it does not contribute to the dynamics and we need only study the time evolution of $`S^z(B=\mathrm{})`$ under $`H_0`$. The imaginary part of the Fourier-transform of $`\chi (t)`$ follows easily

$`\chi ^{\prime \prime }(\omega )`$ $`=`$ $`{\displaystyle \frac{\pi }{16}}{\displaystyle \underset{p,q>0}{}}d_p^2(B=\mathrm{})d_q^2(B=\mathrm{})`$ (19)
$`\times \left(\delta (\omega ϵ_pϵ_q)\delta (\omega +ϵ_p+ϵ_q)\right).`$
Results are shown in fig. 2: The curves display a broad maximum at an energy of the order of the Kondo scale and a power law decay $`\chi ^{\prime \prime }(\omega )\omega ^{(32\alpha )}`$ at high frequencies, consistent with results obtained for the spin boson model . At low frequencies one finds $`\chi ^{\prime \prime }(\omega )\omega `$, where the corresponding slope scales as $`T_K^2`$ with a prefactor of order one . The normalization condition $`_0^{\mathrm{}}𝑑\omega \chi ^{\prime \prime }(\omega )=\pi /4`$ is fulfilled due to the sum rule $`h^2(B)+_pd_p^2(B)/4=1`$ that holds exactly .

The real part $`\chi ^{}(\omega )`$ can be obtained by a Kramers-Kronig transformation. In particular, the *static* susceptibility in response to a local field is given by $`\chi _0=\chi ^{}(\omega =0)`$ and therefore contains spectral information on all energy scales. In the limit of small spin flip coupling $`J_{}`$ we have numerically verified the power law behaviour $`\chi _0J_{}^{1/(\alpha 1)}`$ known from the ohmic spin boson model with the dissipation parameter $`\alpha =(1J_{}/4\pi )^2`$ .

Summing up, we have applied Wegner’s flow equations to the anisotropic Kondo problem. The Hamiltonian is written in terms of vertex operators with scaling dimensions that flow to the Toulouse point. In contrast to the perturbative scaling approach, our flow equations close exactly at the Toulouse point with finite couplings and a nontrivial strong–coupling behaviour. Similar to the analysis of the sine–Gordon model  no strong–coupling divergence of the coupling constants is encountered. We are thus able to give an analytic description of the crossover from weak to strong coupling in a systematic expansion that can be improved by taking higher orders in the OPE into account. The stability of the strong–coupling fixed point cannot be endangered in this expansion.

In our solution of the flow equations we have focussed on the equilibrium dynamics of the impurity at zero temperature. As an example we have calculated the universal scaling forms of the local dynamical susceptibility $`\chi (\omega )`$ and found good agreement with results known from the spin-boson model. This exemplifies the usefulness of our approach for the calculation of correlation functions in strong–coupling problems, and should be of interest in other models as well. Notice that we have not used the integrable structure of the Kondo model in our method. Future work will concentrate on the effect of finite temperature, and on the calculation of static impurity contributions where the newly generated potential scattering terms give the leading contribution (in contrast to the dynamical susceptibility where these terms play no role) and our present approximation must be improved.

The authors acknowledge many valuable discussions with T.A. Costi, D. Fisher, A. Mielke and D. Vollhardt. This work was supported by the Deutsche Forschungsgemeinschaft (DFG), by the SFB 484 of the DFG and by the National Science Foundation (NSF) under grants DMR 9630064, DMR 9976621 and DMR 9981283.


text:
# NRQCD, Effective Field Theories and Potential Models

## Abstract

We review some recent developments on non-relativistic effective field theories for heavy quark-antiquark systems and how they can bring a well founded connection between QCD and potential models.

The study of the heavy quark–antiquark system is an old topic (see  for earlier references). Here we will concentrate on recent developments based on effective field theories. For large enough masses, these systems can be considered to be non-relativistic (NR) and are then characterized by, at least, three widely separated scales: hard (the mass $`m`$, of the heavy quarks), soft (the relative momentum of the heavy-quark–antiquark $`|𝐩|mv`$, $`v1`$), and ultrasoft (US, the typical kinetic energy $`Emv^2`$ of the heavy quark in the bound state system). In 1986, NRQED , an effective field theory for non-relativistic leptons, was presented, providing the first and decisive link in a chain of developments that is still growing. NRQED is obtained from QED by integrating out the hard scale $`m`$ (see  for some applications). It is characterized by an ultraviolet cut-off much smaller than the mass $`m`$ and much larger than any other scale. NRQCD  was born soon afterwards. NRQCD has proved to be extremely successful in studying $`Q`$-$`\overline{Q}`$ systems near threshold. The Lagrangian of NRQCD can be organized in powers of $`1/m`$, thus making explicit the non-relativistic nature of the physical systems.

In order for effective field theories to be useful, a power counting is needed. The first power counting (organized in powers of $`v`$ and $`\alpha _s`$) used to assess the relative importance of the NRQCD matrix elements followed from arguments valid in the perturbative regime . In particular, it was assumed that ΛQCD<

mv2subscriptΛ𝑄𝐶𝐷<

𝑚superscript𝑣2\Lambda_{QCD}{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 6.0pt\vbox{\hbox{$\sim$}}}}\ }mv^{2}. This implies that dynamics is mainly perturbative, with $`v\alpha _s`$, dictated by the solution of the Schrödinger equation with the Coulomb potential $`C_f\alpha _s/r`$. Therefore, it is debatable whether it could be applied for the charmonium system or for higher states of the bottomonium systems. Moreover, even in the perturbative case, because of the different scales involved in the problem, the matrix elements do not have a unique power counting but they also contribute to subleading orders in the $`v`$ counting, whereas the power counting in Ref.  only provided the leading order.

At this stage, there were two major questions (somewhat related) in the formulation of NRQCD: 1) the first was the non-existence of explicit power counting rules to systematically incorporate subleading effects, even in the perturbative situation; 2) the second was that the regularization procedure was based on cut-off regularization. Attempts to perform the matching between QCD and NRQCD using dimensional regularization had the drawback that the naive incorporation of the kinetic term in the quark propagator jeopardizes the power counting rules.

The solution to the latter problem was first given in Ref. . There, it was argued that the matching between QCD and NRQCD in dimensional regularization should be performed just as in HQET –namely, the kinetic term must be treated as a perturbation and then the matching conditions and their Lagrangian in both theories are the same (in particular the matching computation in the effective theory is zero)– and performed for the terms bilinear in the fermions. For the four-fermion terms the matching along the same lines was worked out in Refs. . The key point in the above derivations was that, in order to carry out the matching, it is not so important to know the power counting of each term in the effective theory as to know that the remaining dynamical scales of the effective theory are much lower than the mass. The power counting tells us the relative importance between different operators, but this does not change the value of the matching coefficients. That is, we only need

$$m|𝐩|,E,\mathrm{\Lambda }_{QCD}.$$

As a matter of fact, this defines NRQCD in dimensional regularization (if it exists) irrespective of the relative size between $`\mathrm{\Lambda }_{QCD}`$ and the soft and ultrasoft scale, in particular, even if $`\mathrm{\Lambda }_{QCD}`$ were the next relevant scale.

The first problem, to obtain a complete power counting, was first studied within a mainly perturbative framework and triggered two lines of research:

A) On the one hand, by trying to classify the different momentum regions existing in a pure perturbative version of NRQCD and/or to reformulate NRQCD in ways where some of these regions were explicitly displayed by introducing new fields in the NRQCD Lagrangian. In particular, in Ref. , within a QED context, the existence of the ultrasoft region was first made explicit, besides the soft/potential one (of which the author only considered the leading order contribution somewhat missing the soft region as defined in ), and rules were given to estimate their size. Subsequent work  also tried to get complete power counting rules for NRQCD, but they also missed the soft region. Within this philosophy, in Ref. , the most complete classification of (perturbative) momentum regions to date was made by a rigorous diagrammatic study called the threshold expansion, including the soft region missed in these works. It was then realized that the hard region corresponded to the matching between NRQCD and QCD as described in Refs. . This allowed the computation of the electromagnetic current matching coefficient at two-loops . Finally, some work was also done on formulating the threshold expansion of NRQCD within an effective Lagrangian formalism .

B) On the other hand, in parallel, a different approach was worked out in detail in  to deal with the power-counting problem. It tried to answer the question: How would we like the effective theory for $`\overline{Q}`$$`Q`$ systems near threshold to be? The first observation was that we did not want to describe all the degrees of freedom included in NRQCD, but rather only those with US energy. Therefore, the unwanted degrees of freedom should be integrated out. Moreover, we wanted to get a closer connection with a Schrödinger-like formulation for these systems (see also ). The idea was to connect NRQCD with potential models also, eventually, in the non-perturbative regime (more on that later). Roughly speaking we wanted our effective theory to be something like:

$`\begin{array}{cc}& \left(i_0{\displaystyle \frac{𝐩^2}{2m}}V_0(r)\right)\mathrm{\Phi }(𝐫)=0\hfill \\ & +\mathrm{corrections}\mathrm{to}\mathrm{the}\mathrm{potential}\hfill \\ & +\mathrm{interaction}\mathrm{with}\mathrm{other}\mathrm{low}\hfill \\ & \mathrm{energy}\mathrm{degrees}\mathrm{of}\mathrm{freedom}\hfill \end{array}\}\mathrm{pNRQCD}`$
where $`V_0(r)=C_f\alpha _s/r`$ in the perturbative case and $`\mathrm{\Phi }(𝐫)`$ is the $`\overline{Q}`$$`Q`$ wave-function.

The resulting effective field theory was called potential NRQCD (pNRQCD). This is obtained after integrating out the soft scale, understood in the following way: pNRQCD is defined by two ultraviolet (UV) cut-offs $`\mathrm{\Lambda }_1`$ and $`\mathrm{\Lambda }_2`$. The former fulfils the relation $`mv^2\mathrm{\Lambda }_1mv`$ and is the cut-off of the energy of the fermions and of the energy and the three-momentum of the gluons, whereas the latter fulfils $`mv\mathrm{\Lambda }_2m`$ and is the cut-off of the relative three-momentum of the heavy fermion–antifermion system, $`𝐩`$. This choice of the cutoffs can be motivated as follows. We are only interested in the degrees of freedom with US energy and, as a general philosophy, we would like to keep the UV cut-offs of the effective theory as low as possible. Therefore, any degree of freedom with larger energy should be integrated out. This fixes $`\mathrm{\Lambda }_1`$ (we also fix the three-momentum cut-off of the gluons to be $`\mathrm{\Lambda }_1`$ since the gluons satisfy a relativistic dispersion relation). The motivation for the choice of $`\mathrm{\Lambda }_2`$ is that fermions with US energy have soft three-momentum owing to their non-relativistic dispersion relation. In short, we have only integrated out the degrees of freedom where a perturbative, order by order in $`\alpha `$, matching calculation can be performed (note that heavy fermions with soft three-momentum and US energy cannot be integrated out perturbatively since they are the responsible for producing the bound state). We believe that pNRQCD, defined in this way, has several advantges:

i) The matching between NRQCD and pNRQCD can be done order by order in $`\alpha _s`$ and in $`1/m`$ (by analytically expanding the NRQCD Feynman diagrams in terms of the remaining dynamical scales of pNRQCD there is also an expansion in terms of the ultrasoft scale). In fact, this basically corresponds to obtaining the potential of the Schrödinger-like equation.

ii) The final effective theory resembles very much a Schrödinger-like equation. We believe that pNRQCD as described above provides a rigorous conection between NRQCD (quantum field theories) and potential models (non-relativistic quantum-mechanics formulation) in the situation where ΛQCD<

mv2subscriptΛ𝑄𝐶𝐷<

𝑚superscript𝑣2\Lambda_{QCD}{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 6.0pt\vbox{\hbox{$\sim$}}}}\ }mv^{2}.

iii) In spite of working with a Schrödinger-like equation, US gluons are still incorporated in a second-quantized, systematic and gauge-invariant fashion.

iv) The construction of pNRQCD follows a step-by-step procedure. It goes along the standard effective field theory idea of integrating out scale by scale , although properly speaking it should be said region by region. Since we are dealing with non-relativistic theories, the energy and the three-momentum are not related by the relativistic dispersion relation but by the non-relativistic one, producing asymmetric ultraviolet cut-off in the effective theories. The advantage of working region by region is that in each (perturbative) region the Feynman integrals become homogeneous111I thank M. Beneke for stressing this point to me., and much easier to compute in dimensional regularization, with just the one scale that we want to integrate out (on general grounds the matching calculation in each new effective theory is zero). On the other hand even in a non-perturbative situation it is expected that a better control in the dynamics would be obtained. In particular, perturbative regions can be disantangled from non-perturbative ones in a better controlled way.

Finally we would like to note that, once the UV cut-offs of pNRQCD are fixed, no ambiguity is left in the definition of pNRQCD. Still, for an easy comparison with the first line of research mentioned above, we mention that pNRQCD is obtained from NRQCD after integrating out soft quarks and gluons and potential gluons. We would like to stress here that, in particular, the so-called soft region was already included in Ref.  that, as a matter of fact, was the first place where a power counting including the hard, soft, potential and ultrasoft region was given within an effective field theory framework.

Leaving aside non-perturbative effects, one could be worried about smaller dynamical scales $`mv^3`$,… that could appear by fine-tuning the energy to be near some poles. We would like to note that these eventual scales will not be a problem since they would be encoded in pNRQCD, but then it may be convenient to integrate out some further degrees of freedom to accurately describe physics at these smaller scales.

The developments in effective field theories explained above (with ΛQCD<

mv2subscriptΛ𝑄𝐶𝐷<

𝑚superscript𝑣2\Lambda_{QCD}{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 6.0pt\vbox{\hbox{$\sim$}}}}\ }mv^{2}) can be applied to several physical situations such as the $`\mathrm{{\rm Y}}(1S)`$ , bottomonium sum rules  or $`t`$$`\overline{t}`$ production near threshold , where we believe that put on more solid theoretical grounds the formalism used.

There has been some recent attempts to obtain renormalization group equations for non-relativistic systems. For lack of space we do not discuss them here in detail (see ). We only briefly comment on the work of , from which it could be naively concluded that pNRQCD could not reproduce the anomalous dimension of the electromagnetic current. This is not true, which should be evident from the fact that the evaluation of the anomalous dimension in pNRQCD would go along similar lines to the one performed in . We believe that the approach in  may suffer from the drawback of having to treat the different scales entangled. This is in contradiction with the philosophy advocated here of working region by region and may jeopardize the huge simplification obtained in perturbative calculations. Moreover, since they treat the soft and US scale perturbatively, it is hard to imagine how to introduce non-perturbative effects in this approach.

So far we have restricted our considerations to the situation ΛQCD<

mv2subscriptΛ𝑄𝐶𝐷<

𝑚superscript𝑣2\Lambda_{QCD}{\ \lower-1.2pt\vbox{\hbox{\hbox to0.0pt{$<$\hss}\lower 6.0pt\vbox{\hbox{$\sim$}}}}\ }mv^{2}. It is doubtful whether we can consider most of the charmonium and bottomonium spectrum to be in this situation but rather in the (generic) non-perturbative case with $`mv\mathrm{\Lambda }_{QCD}`$. Then, it is not clear, a priori, what is the power counting that should be used for these systems222In fact, a different, non-standard, power counting of the matrix elements of NRQCD may eventually explain the apparent difficulties that NRQCD is facing to explain the polarization of prompt $`J/\psi `$ data, and to accurately determine the different matrix elements (see ).. In particular, it is less clear how to obtain a rigorous connection between NRQCD and potential models (if it exists), although naively one would expect that, to some extent, the same philosophy as used previously to obtain pNRQCD could also be followed here. In this line of thinking, we would like to end by briefly reporting on some new results in our understanding of NRQCD and pNRQCD in the non-perturbative regime obtained in Ref. ; there, for the first time, a method to obtain a controlled derivation of the potential from QCD in terms of Wilson loops at arbitrary orders in $`1/m`$ has been proposed. This has permitted to obtain new potentials for the heavy quarkonium (as well as for the heavy hybrids) at $`O(1/m)`$ and $`O(1/m^2)`$ previously missed in the literature, as well as to correct some of the already existing ones.

Acknowledgements. I thank J. Soto for his careful reading of the manuscript and M. Beneke for discussions.