text:
# Numerical evidence for ‘multi-scalar stars’

## I Introduction

The nature of dark matter in the universe is currently an open question in physics, with many models being proposed to fill this gap in our understanding, some of which resort to the use of exotic matter. Of interest to us is one class of models composed of massive scalar fields coupled to the general relativistic gravitational field, from which compact, star-like solutions can be formed, solutions which go by the names of “oscillating soliton stars” (or “oscillatons”) OSS ; Lopez  for real fields and “boson stars” Kaup ; RB ; Colpi ; Jetzer ; Mielke ; Bizon  for complex fields. These star-like solutions have received renewed attention recently, and a substantial body of evidence has been advanced in an effort to show that these fields may be key players on both galactic JWLee ; Francisco ; Miguel  and cosmological Matos  scales. Boson stars have been suggested as alternatives to primordial black holes Mielke2  as well as supermassive black holes in the centers of galaxies Torres1 , and their gravitational lensing properties have been explored Dabrowski , further developing the treatment of these solutions as objects of astrophysical interest.

Apart from the possible astrophysical relevance of these star-like objects, we find them interesting to study from a mathematical standpoint as well, for their properties as soliton-like solutions in the nonlinear dynamics of general relativity. The “solution space” of general relativity is still largely unexplored, and these scalar objects comprise simple systems with which to conduct investigations. It is from this viewpoint that we will proceed in this paper.

In 1991, Seidel and Suen OSS  considered the model of a real massive scalar field, minimally coupled to the general relativistic gravitational field, with the additional simplifying assumption of spherical symmetry. These authors were interested in the existence of “nontopological solitons” in the model: that is, whether the equations of motion admitted stable, localized, non-singular distributions of matter which could be interpreted as “scalar stars”. A theorem due to Rosen Rosen  suggested that, should such solutions exist, they could not be static. Thus, Seidel and Suen looked for periodic configurations by substituting a particular Fourier ansatz into the equations of motion, and solving the resulting hierarchy of ODEs via a generalized shooting technique. The authors found strong evidence that periodic star-like solutions did exist, and, via direct numerical simulation, demonstrated that their “oscillating soliton stars”, if not absolutely stable, had lifetimes many orders of magnitude longer than the stars’ intrinsic dynamical times.111More precisely, the oscillating stars constitute a one-parameter family which may be parametrized by the mean (period-averaged) central density, $`\rho _0`$. As with other relativistic stellar models, a plot of total (ADM) mass versus $`\rho _0`$ exhibits a maximum at $`\rho _0^{}`$, which seems to coincide with a transition from stable to unstable configurations. As expected, only stars with $`\rho _0<\rho _0^{}`$ could be stably evolved for long times.

These results were surprising to some researchers, particularly since the model has no conserved Noether current, which, it had been argued, was responsible for the existence of solitonic solutions in other non-linear field theories involving scalar fields Coleman ; Lee . However, at least heuristically, we can understand the existence of the oscillating stars as arising from a balance between the attractive gravitational interaction and the effective repulsive self-interaction generated by the mass of the scalar field (i.e. via the dispersion relation of the Klein-Gordon equation).

Recently it was shown by Ureña-López Lopez  that approximate solutions for boson stars and oscillating soliton stars, or “oscillatons” as he calls them, can both be derived from a single set of equations in a sort of “stationary limit”. The similarities seen between boson stars and oscillating soliton stars in terms of their curves relating mass, radius and central density can thus be related formally.

In this paper, we build on the works of Seidel and Suen and Ureña-López by considering a matter content consisting of multiple scalar fields, and we find some further ways in which boson stars and oscillating soliton stars are similar. For the specific case of two scalar fields, we find evidence for a new family of quasi-periodic, solitonic configurations.222Although the oscillations appear to be periodic, we do not have a proof of their periodicity and so we use the term “quasi-periodic” to describe their temporal behavior. Small departures from strict periodicity over long times scales are visible in the simulation results; however, these departures become smaller as we decrease the radial mesh spacing $`\mathrm{\Delta }r`$, and thus it is plausible that in the limit $`\mathrm{\Delta }r0`$, the solutions are truly periodic. Together with previous results, this suggests that solitonic solutions are generic to models which couple massive scalar fields through the Einstein gravitational field. We should note, however, that we make no attempt to address the question of whether these multiple scalar fields actually exist in nature; rather we are interested in their existence as valid mathematical solutions in the Einstein-Klein-Gordon system.

We begin by considering $`n`$ real, massive Klein-Gordon fields $`\varphi _i,i=1,2,\mathrm{}n`$, without additional self-interaction, minimally coupled to the general relativistic gravitational field. Specifically, choosing units such that $`c=1`$ and $`G=1`$, the Lagrangian density for the coupled system is

$$=\sqrt{g}\left(R\underset{i=1}{\overset{n}{}}\left(\varphi _{i}^{}{}_{}{}^{;a}\varphi _{i}^{}{}_{;a}{}^{}m_i^2\varphi _i^2\right)\right)$$
(1)
where $`gdetg_{\mu \nu }`$, $`R`$ is the Ricci scalar, and $`m_i`$ is the mass of the $`i`$-th scalar field.

We now restrict attention to spherical symmetry, and adopt the “polar/areal” coordinate system, so that the metric takes the form:

$$ds^2=\alpha ^2(t,r)dt^2+a^2(t,r)dr^2+r^2d\mathrm{\Omega }^2$$
(2)
The complete evolution of the scalar fields and spacetime can then be given in terms of a Klein-Gordon equation for each of the $`\varphi _i`$, and two constraints derived from the Einstein field equations and the coordinate conditions used to maintain the metric in the form (2). We solve these equations using the same scheme adopted for the critical phenomena study described in HawleyChop , and only briefly review that scheme here.

We define the following auxiliary scalar field variables:

$$\mathrm{\Phi }_i\varphi _i^{},\mathrm{\Pi }_i\frac{a}{\alpha }\dot{\varphi }_i$$
(3)
where all variables are functions of $`t`$ and $`r`$, $`\dot{}/t`$ and $`{}_{}{}^{}/r`$. The Klein-Gordon equation is written as the following system:

$`\dot{\mathrm{\Pi }}_i`$ $`=`$ $`{\displaystyle \frac{1}{r^2}}\left({\displaystyle \frac{r^2\alpha }{a}}\mathrm{\Pi }_i\right)^{}m_i^2\alpha a\varphi _i`$ (4)
$`\dot{\mathrm{\Phi }}_i`$ $`=`$ $`\left({\displaystyle \frac{\alpha }{a}}\mathrm{\Pi }_i\right)^{}`$ (5)
$`\varphi _i(t,r)`$ $`=`$ $`\varphi _i(t,r_{\mathrm{max}}){\displaystyle _{r_{\mathrm{max}}}^r}\mathrm{\Phi }(t,\stackrel{~}{r})𝑑\stackrel{~}{r}`$ (6)
where $`r=r_{\mathrm{max}}`$ is the outer boundary of the computational domain. The constraint equations are the “Hamiltonian constraint”

$$a^{}=a\frac{1a^2}{2r}+2\pi ra\underset{i=1}{\overset{n}{}}\left(\mathrm{\Pi }_{i}^{}{}_{}{}^{2}+\mathrm{\Phi }_{i}^{}{}_{}{}^{2}+a^2m_i^2\varphi _i^2\right),$$
(7)
and the “slicing condition”

$$\alpha ^{}=\alpha \left(\frac{a^21}{r}+\frac{a^{}}{a}4\pi ra^2\underset{i=1}{\overset{n}{}}m_i^2\varphi _i^2\right).$$
(8)
For diagnostic purposes, we have also found it useful to compute and monitor the quantities, $`M_i(t,r_{\mathrm{max}})`$, defined by

$$M_i(t,r_{\mathrm{max}})4\pi _0^{r_{\mathrm{max}}}\stackrel{~}{r}^2\rho _i(t,\stackrel{~}{r})𝑑\stackrel{~}{r},$$
(9)
where

$$\rho _i(t,r)=\frac{\mathrm{\Pi }_i{}_{}{}^{2}+\mathrm{\Phi }_i{}_{}{}^{2}+a^2\varphi _i^2}{2a^2}.$$
(10)
Loosely speaking, we can interpret $`M_i(t,r_{\mathrm{max}})`$ as the total contribution of field $`i`$ to the ADM mass of the spacetime. In particular, so long as no matter out-fluxes through $`r=r_{\mathrm{max}}`$, we have

$$\underset{i=1}{\overset{n}{}}M_i(t,r_{\mathrm{max}})=\text{const.}$$
(11)
We solve Eqs. (6)-(8) subject to the the boundary conditions $`a(t,0)=1`$ (local flatness at the origin) and $`\alpha (t,r_{\mathrm{max}})=1/a(t,r_{\mathrm{max}})`$ (so that $`t`$ measures proper time as $`r\mathrm{}`$). As in HawleyChop , we use the Sommerfeld condition for a massless field to set the values $`\varphi (t,r_{\mathrm{max}}),\mathrm{\Phi }(t,r_{\mathrm{max}})`$ and $`\mathrm{\Pi }(t,r_{\mathrm{max}})`$. Since the Sommerfeld condition is not ideal for a massive field, we ran our simulations with different values of $`r_{\mathrm{max}}`$, testing for any periodicity or other effect which might be a function of $`r_{\mathrm{max}}`$, and usually ran with an $`r_{\mathrm{max}}`$ which was large compared to the time for which we ran the simulation. Even with smaller $`r_{\mathrm{max}}`$, we found our results to be essentially independent of $`r_{\mathrm{max}}`$ and attribute this to the fact that there is very little scalar radiation emitted from the soliton-like objects considered here. Our results are also essentially independent of the resolution of the finite differencing algorithm and the Courant-Friedrichs-Levy factor, $`\mathrm{\Delta }t/\mathrm{\Delta }r`$, and we confirm that our results converge in a second-order-accurate manner using independent residual evaluations.

## II “Phase-Shifted Boson Stars”

We start by considering the case $`n=1`$, so that our matter content is a single scalar field, $`\varphi (t,r)`$. We note that the Hamiltonian constraint (7) and the slicing condition (8) are unchanged if we trivially decompose $`\varphi `$ into two identical fields (i.e. now choosing $`n=2`$), namely $`\varphi _1(t,r)`$ and $`\varphi _2(t,r)`$ = $`\varphi _1(t,r)`$ (with $`m_2=m_1=1`$), such that

$$\varphi =\frac{1}{\sqrt{2}}\left(\varphi _1+\varphi _2\right).$$
(12)
Further, we note that for fixed $`\alpha (t,r)`$ and $`a(t,r)`$, if $`\varphi (t,r)`$ is a solution of the Klein-Gordon equation (6), then so is $`\kappa \varphi (t,r)`$ where $`\kappa `$ is an arbitrary real constant. Since a soliton solution of the system (6)-(8) is the oscillating soliton star, we see that a trivial multi-scalar soliton solution can be obtained by constructing an oscillating soliton star with a single field $`\varphi `$ (as described in OSS ) and then reinterpreting it as a two-field solution in which $`\varphi _1\varphi _2\varphi /\sqrt{2}`$.

Moreover, if we wish to model a boson star Kaup ; RB ; Colpi ; Jetzer ; Mielke  with no self-interaction potential (often called a “mini-boson star” as in Bizon ), then we have one massive complex scalar field $`\stackrel{~}{\varphi }(t,r)`$, for which the real and imaginary parts behave like two real-valued scalar fields: $`\stackrel{~}{\varphi }(t,r)=\varphi _1(t,r)+i\varphi _2(t,r)`$. The boson star ansatz is $`\stackrel{~}{\varphi }(t,r)=\widehat{\varphi }(r)\mathrm{exp}(\pm i\omega t),`$ where $`\widehat{\varphi }(r)`$ is real. This implies

$`\varphi _1(t,r)`$ $`=`$ $`\widehat{\varphi }(r)\mathrm{cos}(\omega t)`$ (13)
$`\varphi _2(t,r)`$ $`=`$ $`\widehat{\varphi }(r)\mathrm{cos}(\omega t+\delta ),`$ (14)
where $`\delta =\pi /2`$.

Thus we see that soliton stars and boson stars can both be obtained using two real scalar fields with a constant temporal phase difference. For soliton stars, the fields are identical for all $`r`$ and $`t`$; whereas for boson stars, the fields have identical $`r`$-dependence, and the $`t`$-dependence is the same to within a phase. (In each case, the central density $`\rho _0`$ is uniquely fixed by the value of the field at the origin, e.g. $`\widehat{\varphi }(0)`$.)

The work described in the remainder of this paper began in the midst of our numerical evolutions of boson stars HawleyChop . The question arose, “What happens if we solve for the boson star initial data, then ‘manually’ change the phase between the two fields, and then re-solve the constraints to obtain the metric variables?” For future reference, we term such a modified-boson-star configuration a “phase-shifted boson star.” This modification to the boson star data was motivated more by “practical” reasons that “physical” ones — we were interested in studying oscillating soliton stars but found them difficult to construct.

Taking the boson star initial data and setting $`\varphi _2(0,r)\varphi _1(0,r)`$ resulted in what might be termed a “poor man’s soliton star” (PMSS). In Figures 1, 2 and 3, we show that the resulting solution is very similar to the true soliton star solution, and can perhaps best be regarded as a soliton star with a small perturbation added.

We then took the boson star initial data $`\widehat{\varphi }(r)`$ and distributed it to $`\varphi _1`$ and $`\varphi _2`$ via Eq. (14) using some different value of $`\delta `$, such as $`\delta =\pi /6`$. (Note that we only apply Eq. (14) for the initial data, i.e. at $`t=0`$; one cannot expect Eq. (14) to describe the fields for all $`t`$ if $`\delta \pm \pi /2`$, as the $`\delta =0`$ case demonstrates.) One aspect of the evolution for such a system can be seen in Figure 4.

For each of the many values of $`\delta `$ we tried, we found an apparently stable solution which oscillated in some essentially periodic manner for very long times. (The phase was preserved throughout the evolution, i.e. it is not the case that the system reverted to a simple “perturbed boson star” over time.) These results led us to conjecture that there may exist a continuous family of periodic soliton-like solutions (parameterized by the phase $`\delta `$) of which our “phase-shifted boson stars” are perturbations. We hope in the future to construct such a family directly via a periodic ansatz of the form used by Seidel and Suen for their oscillating soliton stars (with additional terms incorporated to account for the nonlinear coupling between the two scalar fields).

We would like to mention that, for a given $`\widehat{\varphi }(0)`$, varying $`\delta `$ has very little effect on the metric variable $`a(0,r)`$. Consequently the so-called “stability curves” relating total mass, radius and central density are essentially the same as the curves for boson stars ($`\delta =\pm \pi /2`$); consequently, we do not consider it relevant to plot them here.

## III Conclusions

We have strong numerical evidence for the existence of a two-parameter family of soliton-like solutions to the Einstein-Klein-Gordon system (parameterized by central density $`\rho _0`$ and phase $`\delta `$) for the case of two scalar fields. We speculate that extending the system to more scalar fields will yield similar results (especially if one uses a trivial extension like Eq. (12) ). The solutions we refer to as “phase-shifted boson stars” consist of boson star initial data for which the phase difference $`\delta `$ between the real and imaginary components of the field have been altered. These solutions oscillate in a seemingly periodic manner for very long times; thus they appear to be stable. For the case of $`\delta =0`$, we obtain close approximations to the oscillating soliton stars of Seidel and Suen. For other values of $`\delta `$, we find solutions which also appear to be stable and periodic; furthermore we can see mass-energy being exchanged between the two fields.

###### Acknowledgements.

We wish to thank Wai-Mo Suen and Edward Seidel for their valuable input and for sharing the computer code they used to generate initial data for oscillating soliton stars. S.H. Hawley would also like to thank Alan Rendall, Piotr Bizon and Francisco Guzman for helpful discussions. This research was supported by NSERC and CIAR. To perform the simulations referred to in this paper, we used the Silicon Graphics workstations of the Center for Relativity at the University of Texas at Austin, the CFI-funded vn Linux cluster of the University of British Columbia, and a laptop computer provided by the Albert Einstein Institute.


format: text

text:
# A study of QCD coupling constant and power corrections in the fixed target deep inelastic measurements Talk presented at X Int. Workshop on Deep Inelastic Scattering (DIS2002), Cracow, 30 April - 4 May 2002. A shorter version will appear in the Proceedings.

## 1  Introduction

The deep inelastic scattering (DIS) leptons on hadrons is the basical process to study the values of the parton distribution functions (PDF) which are universal (after choosing of factorization and renormalization schemes) and can be used in other processes. The accuracy of the present data for deep inelastic structure functions (SF) reached the level at which the $`Q^2`$-dependence of logarithmic QCD-motivated terms and power-like ones may be studied separately (for a review, see  and references therein).

In the present paper we sketch the results of our analysis  at the next-to-leading order (NLO) of perturbative QCD for the most known DIS SF $`F_2(x,Q^2)`$ taking into account experimental data - of SLAC, NM, BCDMS and BFP Collaborations. We stress the power-like effects, so-called twist-4 (i.e. $`1/Q^2`$) contributions. To our purposes we represent the SF $`F_2(x,Q^2)`$ as the contribution of the leading twist part $`F_2^{pQCD}(x,Q^2)`$ described by perturbative QCD, when the target mass corrections are taken into account, and the nonperturbative part (“dynamical” twist-four terms):

$$F_2(x,Q^2)F_2^{full}(x,Q^2)=F_2^{pQCD}(x,Q^2)\left(1+\frac{\stackrel{~}{h}_4(x)}{Q^2}\right),$$
(1)
where $`\stackrel{~}{h}_4(x)`$ is magnitude of twist-four terms.

Contrary to standard fits (see, for example, \- ) when the direct numerical calculations based on Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation  are used to evaluate structure functions, we use the exact solution of DGLAP equation for the Mellin moments $`M_n^{tw2}(Q^2)`$ of SF $`F_2^{tw2}(x,Q^2)`$ and the subsequent reproduction of $`F_2^{full}(x,Q^2)`$ at every needed $`Q^2`$-value with help of the Jacobi Polynomial expansion method 111We note here that there is similar method , based on Bernstein polynomials. The method has been used in the analyses at the NLO level in  and at the NNLO level in .  (see similar analyses at the NLO level  and at the next-next-to-leading order (NNLO) level and above .

In this paper we do not present exact formulae of $`Q^2`$-dependence of SF $`F_2`$ which are given in . We note only that the PDF at some $`Q_0^2`$ is theoretical input of our analysis and the twist-four term $`\stackrel{~}{h}_4(x)`$ is considered as a set of free parameters (one constant $`\stackrel{~}{h}_4(x_i)`$ per $`x_i`$-bin): $`\stackrel{~}{h}_4^{free}(x)=_{i=1}^I\stackrel{~}{h}_4(x_i)`$, where $`I`$ is the number of bins.

## 2  Fits of $`F_2`$

First of all, we choose the cut $`Q^21`$ GeV2 in all our studies. For $`Q^2<1`$ GeV2, the applicability of twist expansion is very questionable. Secondly, we choose $`Q_0^2`$ = 90 GeV2 ($`Q_0^2`$ = 20 GeV2) for the nonsinglet (combine nonsinglet and singlet) evolution, i.e. quite large values of the normalization point $`Q_0^2`$: our perturbative formulae should be applicable at the value of $`Q_0^2`$. Moreover, the higher order corrections $`\alpha _s^k(Q_0^2)`$ and $`(\alpha _s(Q^2)\alpha _s(Q_0^2))^k`$ ($`k2`$) should be less important at these $`Q_0^2`$ values.

We use MINUIT program  for minimization of $`\chi ^2(F_2)`$. We consider free normalizations of data for different experiments. For the reference, we use the most stable deuterium BCDMS data at the value of energy $`E_0=200`$ GeV ($`E_0`$ is the initial energy lepton beam). Using other types of data as reference gives negligible changes in our results. The usage of fixed normalization for all data leads to fits with a bit worser $`\chi ^2`$.

### 2.1  BCDMS $`{}_{}{}^{12}C+H_2+D_2`$ data

We start our analysis with the most precise experimental data  obtained by BCDMS muon scattering experiment at the high $`Q^2`$ values. The full set of data is 762 points.

It is well known that the original analyses given by BCDMS Collaboration itself (see also Ref. ) lead to quite small values $`\alpha _s(M_Z^2)=0.113`$. Although in some recent papers (see, for example, ) more higher values of the coupling constant $`\alpha _s(M_Z^2)`$ have been observed, we think that an additional reanalysis of BCDMS data should be very useful.

Based on study  we proposed in  that the reason for small values of $`\alpha _s(M_Z^2)`$ coming from BCDMS data was the existence of the subset of the data having large systematic errors. We studied this subject by introducing several so-called $`Y`$-cuts 222Hereafter we use the kinematical variable $`Y=(E_0E)/E_0`$, where $`E`$ is scattering energies of lepton. (see ):

$`y0.14\text{ when }0.3<x0.4,y0.16\text{ when }0.4<x0.5`$
$`yY_{cut3}\text{ when }0.5<x0.6,yY_{cut4}\text{ when }0.6<x0.7`$
$`yY_{cut5}\text{ when }0.7<x0.8`$ (2)
and several $`N`$ sets for the cuts at $`0.5<x0.8`$:

From the Figs. 1 and 2 we can see that the $`\alpha _s`$ values are obtained for $`N=1÷6`$ of $`Y_{cut3}`$, $`Y_{cut4}`$ and $`Y_{cut5}`$ are very stable and statistically consistent.

### 2.2  SLAC, BCDMS, NM and BFP data

After these $`Y`$-cuts have been incorporated (with $`N=6`$) for BCDMS data, the full set of combine data is 1309 points. The results of the fits are compiled in Summary.

## 3  Summary

We have demonstrated several steps of our study  of the $`Q^2`$-evolution of DIS structure function $`F_2`$ fitting all fixed target experimental data.

From the fits we have obtained the value of the normalization $`\alpha _s(M_Z^2)`$ of QCD coupling constant. First of all, we have reanalyzed the BCDMS data cutting the range with large systematic errors. As it is possible to see in the Figs. 1 and 2, the value of $`\alpha _s(M_Z^2)`$ rises strongly when the cuts of systematics were incorporated. In another side, the value of $`\alpha _s(M_Z^2)`$ does not dependent on the concrete type of the cut within modern statistical errors.

Fitting SLAC, BCDMS, NM and BFP data, we have found in  that at $`Q^210÷15`$ GeV2 the formulae of pure perturbative QCD (i.e. twist-two approximation together with target mass corrections) are in good agreement with all data 333We note that at small $`x`$ values, the perturbative QCD works well starting with $`Q^2=1.5÷2`$ GeV2 and higher twist corrections are important only at very low $`Q^2`$: $`Q^20.5`$ GeV2 (see  and references therein). As it is was observed in  (see also discussions in ) the good agreement between perturbative QCD and experiment seems connect with large effective argument of coupling constant at low $`x`$ range.. When we have added twist-four corrections, we have very good agreement between QCD (i.e. first two coefficients of Wilson expansion) and the data starting already with $`Q^2=1`$ GeV2, where the Wilson expansion should begin to be applicable. The results for $`\alpha _s(M_Z^2)`$ are very similar (see ) for the both types of analyses and have the following form:

$`\alpha _s(M_Z^2)`$ $`=`$ $`0.1177\pm 0.0007\text{(stat)}\pm 0.0021\text{(syst)}\pm 0.0009\text{(norm)},`$ (3)
where the symbols “stat”, “syst” and “norm” mark the statistical error, systematic one and the error of normalization of experimental data.

We would like to note that we have good agreement also with the analysis  of combined H1 and BCDMS data, which has been given by H1 Collaboration very recently. Our results for $`\alpha _s(M_Z^2)`$ are in good agreement also with the average value for coupling constant, presented in the recent studies (see  and references therein) and in famous reviews .

At the end of our paper we would like to discuss the contributions of higher twist corrections. In our study  we have reproduced well-known $`x`$-shape of the twist-four corrections at the large and intermediate values of Bjorken variable $`x`$ (see the Fig. 3 and ).

Note that there is a small-$`x`$ rise of the magnitude of twist-four corrections, when we use flat parton distributions at $`x0`$. As we have discussed in Ref , there is a strong correlation 444 This correlation comes because of very limited numbers of experimental data used here lie at the low $`x`$ region. Indeed, only the NMC experimental data contribute there. We hope to incorporate the HERA data  in our future investigations. between the small-$`x`$ behavior of twist-four corrections and the type of the corresponding asymptotics of the leading-twist parton distributions. The possibility to have a singular type of the asymptotics leads (in our fits) to the appearance of the rise of sea quark and gluon distributions as $`x^{0.18}`$ at low $`x`$ values, that is in full agreement with low $`x`$ HERA data and with theoretical studies . At this case the rise of the magnitude of twist-four corrections is completely canceled. This cancellation is in full agreement with theoretical and phenomenological studies (see ).

Acknowledgments. The study is supported in part by the Heisenberg-Landau program. A.V.K. would like to express his sincerely thanks to the Organizing Committee of the X Int. Workshop on Deep Inelastic Scattering (DIS2002) for the kind invitation and for fruitful discussions. He was supported in part by Alexander von Humboldt fellowship and INTAS grant N366.