text:
# COMPARISON OF EISCAT RADAR DATA ON SPACE DEBRIS WITH MODEL PREDICTIONS BY THE MASTER MODEL OF ESA

## INTRODUCTION

Currently ESA monitors the space debris population in low Earth orbit (LEO) in irregular intervals (Mehrholz et al., 2002). This is done in the frame of 24 hour radar beam-park experiments. The data from these measurements are taken in order to constrain models of the population evolution of space debris in orbit (Klinkrad et al., 2001). The most prominent of the European space debris models is MASTER/PROOF’99, which allows spacecraft operators to assess the risk of a debris impact for their mission. As the 24 hour experiments only provide a snapshot of the environment, the temporal behaviour of the space debris population is difficult to determine. Routine monitoring of LEO is required in order to better constrain the debris models. The most suitable technique for monitoring LEO is radar detection by high-power, large aperture facilities.

In order to obtain inexpensive access to such observations, ESA started in 2000 a study with the Sodankylä Geophysical Observatory (SGO) to use the radar facilities of the European Incoherent SCATter (EISCAT) scientific association. These facilities operate up to 2000 hours per year for ionospheric research purposes. In the study, the SGO team was able to demonstrate the exploitation of ionospheric measurements for space debris observations. A comparison of the retrieved, still preliminary, data with the ESA MASTER/PROOF’99 model is presented here.

## THE EISCAT FACILITY AT TROMSØ

The space debris measurements were performed in the frame of a demonstration campaign in February 2001. A radar facility located at $`69^{}35^{}\mathrm{N}`$ and $`19^{}14^{}\mathrm{E}`$ near Tromsø, Norway, was used. The UHF radar at Tromsø operates at a frequency of $`931\mathrm{MHz}`$ with a pulse repetition frequency of $`200\mathrm{kHz}`$. The transmitter can transmit peak powers up to $`1.5\mathrm{MW}`$ through the main antenna, but it was restricted to just below $`1\mathrm{MW}`$ during the measurement campaign. The $`35\mathrm{m}`$ antenna was pointed parallel to the field lines of the local geo-magnetic field, azimuth $`183^{}18^{}`$ and elevation $`77^{}6^{}`$. During the test campaign two different phase codes called TAU2 and CP1LT were transmitted. Of the TAU2 data $`1.2`$ hours have been analysed as well as $`2.8`$ hours of the CP1LT data. For both experiments the pulse length was about $`0.5\mathrm{ms}`$. The range window of the CP1LT experiment covered ranges from $`400`$ to $`1,400\mathrm{km}`$ with a gap between $`700`$ and $`800\mathrm{km}`$, and TAU covered ranges from $`500`$ to $`1,800\mathrm{km}`$ with a gap between $`800`$ and $`1,100\mathrm{km}`$.

## CATALOGUED OBJECT

The strongest of the high-altitude TAU2 events during the test campaign began at 22:19:06 UTC, 20 February 2001. We identified the target as a large catalogued object, a Tsyklon upper stage with the COSPAR designator 1994-11G. According to the catalogue the object has a total mass of $`1,390\mathrm{kg}`$ and is cylindrical in shape, with a diameter of $`2.7\mathrm{m}`$, a height of $`2.2\mathrm{m}`$, and a radar cross section of $`8.3\mathrm{m}^2`$. The orbital information in the TLE catalogue, with TLE epoch one day before the experiment epoch, lets us predict that it should have passed the radar beam at the off-axis angular distance of $`1^{}16^{}`$. Figure 1 shows the analysis summary plot of the event. The analysis was done using $`0.27\mathrm{s}`$ coherent integration.

The top panel of Figure 1 shows a quantity proportional to the signal to noise power ratio during the beam transit. The antenna side-lobe structure is clearly visible. Markers indicate those scans where the ratio was larger than the detection threshold of $`4.5`$. The off-centre angular distance of the pass was fit to the data (Markkanen et al., 2002). The best fit was achieved by assuming that the transit occurred $`31^{}`$ off-axis. It seems difficult to reconcile the predicted offset $`1^{}16^{}`$ with the inferred value of $`31^{}`$.

The middle and bottom panels of Figure 1 show the measured slant range and Doppler velocity. The solid dark curves represent quadratic and linear fits to the data points. The large dots represent the predicted values. The measured range is about $`7\mathrm{km}`$ larger than predicted. The slope of the velocity curve is as predicted, but there is a discrepancy of about $`0.1\mathrm{km}\mathrm{s}^1`$ in the actual velocity values. The circular-orbit, vertical-beam acceleration is $`27.5\mathrm{m}\mathrm{s}^2`$, which is consistent with the value $`27.1\mathrm{m}\mathrm{s}^2`$ from the velocity fit. The timing accuracy for the catalogued objects is of the order of $`10\mathrm{s}`$, while we believe the measured timing to be accurate to within about $`0.1\mathrm{s}`$. However, the range and velocity discrepancies cannot be removed simply by adjusting the relative timing, because the required correction would be about six seconds for the range data, but only about three seconds for the Doppler data. A small misalignment of the antenna can, however, explain the observed deviation.

This example of the detection of a large catalogued object by the EISCAT radar demonstrates the usefulness of routine space debris observations for the validation of large-object catalogues.

## SMALL-SIZE DEBRIS

While collisions of operational spacecraft with large catalogued objects like the upper stage discussed above can be avoided by manoeuvring the spacecraft, objects smaller than $`10\mathrm{cm}`$, that normally are not found in the catalogue, pose a serious threat to manned and unmanned spacecraft. Complex spacecraft like the International Space Station Alpha are shielded from impacts of objects smaller than about $`1\mathrm{cm}`$. Thus, the population in the size range between $`1`$ and $`10\mathrm{cm}`$ is particularly important to observe.

In the demonstration campaign performed by SGO, $`56`$ objects in the size range above $`1\mathrm{cm}`$ were detected within the $`3\mathrm{hours}`$ of analysed data. The measured distribution in size and altitude distributions discussed below are compared to the prediction by the ESA MASTER/PROOF’99 model. Due to the pointing of the radar beam along the local geo-magnetic field line, orbital parameters could not be determined from the measurements. The predicted distribution of orbital elements is discussed below anyhow.

### Size Distribution

Figure 2 shows the distribution of object diameter, as determined by the analysis of the reflected signal amplitude, as well as the prediction by the MASTER/PROOF’99 model. Only objects detected in the range window from $`800`$ to $`1,400\mathrm{km}`$ have been considered.

Under the assumption of incoherent integration, objects larger than about $`7\mathrm{cm}`$ are detectable reliably in the range window between $`800`$ and $`1,400\mathrm{km}`$. While the majority of the objects predicted to cross the beam are smaller than $`2\mathrm{cm}`$, sizes of $`2`$ to $`3\mathrm{cm}`$ have been derived for the objects actually detected in the experiment. From the figure it is evident that while no object with a derived size larger than $`6\mathrm{cm}`$ has been detected, more than $`20`$ of those large objects should have crossed the beam according to the prediction by MASTER/PROOF’99. Also, the number of $`3\mathrm{cm}`$ objects is three times higher than predicted, while the total number of detected objects is in agreement. Apparently, the sizes derived for the objects larger than $`10\mathrm{cm}`$ have been underestimated. Because it is unknown where the objects cross the beam pattern, the EISCAT radars can only give a lower bound for the radar cross section. Consequently, the translation of the measured radar cross section to object diameter systematically underestimates the object size.

### Altitude Distribution

The distribution of orbital altitudes of detected objects is shown in Figure 3. At the altitudes of maximum debris density between $`850`$ and $`950\mathrm{km}`$, $`12`$ objects were detected. According to the MASTER/PROOF’99 model, the predicted number of objects detectable in this altitude bin by incoherent integration is below $`8`$, significantly less than actually detected. This is also true for all but the altitude bin around $`1,200\mathrm{km}`$. The sensitivity of the radar appears to be higher than for the standard incoherent detection method. This can be explained by the fact that the detection algorithm exploits the phase information by coherently integrating the detected signal over $`0.1\mathrm{s}`$. MASTER/PROOF’99 can only simulate radar detection with incoherent integration. Since the detection algorithm used by the EISCAT team is based on a coherent integration technique, the predicted number of objects is lower than the EISCAT data.

### Inclination Distribution

Assuming circular orbits, the inclination of an object’s orbit can in theory be derived from the Doppler velocity measured by the radar. However, the functional dependence of the Doppler velocity on the orbit inclination varies with the orientation of the radar beam. As can be seen in Figures 4 and 5 the mapping from Doppler velocity to orbit inclination is ambiguous for an antenna azimuth around $`180^{}`$ (South). An easterly pointing (antenna azimuth $`90^{}`$) will provide an unambiguous relationship between the measured Doppler velocity and orbit inclination.

As the EISCAT transmitter at Tromsø was pointed parallel to the geo-magnetic field lines, which follow the north-south direction, orbit inclinations can not be derived from the data obtained during the measurement campaign. Using the PROOF application of the MASTER/PROOF’99 model, we can however predict the inclination distribution expected for the measurement campaign. The prediction is shown in Figure 6. It can be seen that inclinations from $`60^{}`$ to $`110^{}`$ are covered by beam-park experiments with the transmitter at Tromsø. This covers debris created in classical Sun-synchronous orbits as well as debris created by most launches from Plesetsk.

## CONCLUSION AND OUTLOOK

The test campaign performed by SGO resulted in raw data on one large space debris object as well as in the size and altitude distribution of $`55`$ small space debris objects. The comparison with the ESA MASTER/PROOF’99 model shows that the conversion from signal amplitudes to object sizes has to be reviewed. The more straight forward measurement of the objects’ altitude provides a distribution that is well in accord with the prediction by the model. As the model is regularly checked against results from other beam-park experiments, we can conclude that the demonstration measurements with the EISCAT transmitter at Tromsø provides reliable measurements. The absolute number of detections is significantly larger than the predicted number when assuming incoherent integration as a detection technique. Thus, the coherent integration of the received signal indeed increases the detection sensitivity.

As expected, we find that, when linked to ionospheric measurements, we will not obtain orbital information using the EISCAT radars. This restriction can, however, easily be removed if dedicated observation time is acquired. In this case an optimum antenna pointing will allow us to determine at least Doppler inclinations, which are identical to true orbit inclinations for circular orbits. As the PROOF tool of the new MASTER/PROOF-2001 model allows to predict the distribution of Doppler inclinations, this would allow us to constrain the statistical orbit properties of the small-size LEO space debris population.

In summary we find that the measurements performed during the demonstration campaign in February 2001 prove the value of space debris data obtained by exploiting ionospheric measurements. The mere amount of data, that will be available in case simultaneous operations with ionospheric measurements in the order of $`1000\mathrm{hours}`$ per year can be performed, provides an important contribution to the understanding of the evolution of LEO space debris.

In order to realise measurements that operate simultaneously with the ionospheric investigations, a real-time detection method has to be developed, because the volume of data that would be needed to record the raw signal from the radar for $`1000\mathrm{hours}`$ per year is prohibitive. With real-time detection only the object parameters as well as the signal of the passage will be recorded. A study to develop real-time measurements will be started by the end of 2002.


format: text

text:
# Centrality Scaling of the 𝑝_𝑇 Distribution of Pions

## Abstract

From the preliminary data of PHENIX on the centrality dependence of the $`\pi ^0`$ spectrum in $`p_T`$ at midrapidity in heavy-ion collisions, we show that a scaling behavior exists that is independent of the centrality. It is then shown that $`p_T`$ degrades with increasing $`N_{\mathrm{part}}`$ exponentially with a decay constant that can be quantified. A scaling distribution in terms of an intuitive scaling variable is derived that is analogous to the KNO scaling. No theoretical models are used in any part of this phenomenological analysis.

In a recent paper hy  we reported on the finding of a scaling property of the $`p_T`$ distribution of pions produced in heavy-ion collisions that is independent of the collision energy. Here we present an extension of that scaling property to include centrality variations and show that a KNO-type kno  scaling behavior exists over the entire range of $`p_T`$ measured. The investigation is primarily a phenomenological analysis with no assumptions about the hard and soft collisions, nor about the parton energy losses.

Recently, a scaling behavior of the transverse-mass spectrum has been reported in jsb . That work was motivated by color glass condensate and the saturation of the gluon density in nuclear collisions. Our investigation has no theoretical motivation other than the search for the simplest form that can represent the data. The dynamical origin of the $`p_T`$ distribution is complicated. At low $`p_T`$ the statistical model seems to work well, as does the hydrodynamical description up to $`p_T=`$3 GeV/c ph . At high $`p_T`$ hard parton scattering will create jets, which can lose energy due to multiple scatterings of partons in the dense medium xnw . A universal description of the hadron distribution over all $`p_T`$ is nonexistent, if not meaningless from the point of view of the sectarian nature of the dynamical theories that claim validities in different domains. However, if a universal scaling behavior can be found phenomenologically, it can serve as a common goal for different dynamical approaches to aim at.

From the preliminary PHENIX data of $`\pi ^0`$ produced in Au+Au collisions at the relativistic heavy-ion collider (RHIC) dde , we have the pion distribution, $`(2\pi p_T)^1dN_\pi /dp_T`$, at midrapidity for $`\sqrt{s}=200`$ GeV and for a wide range of centrality that has 9 bins from 0-10% to 80-92%. To unify the 9 distributions, it is necessary to define a scaling variable $`z`$. First, we use the number of participants, $`N_{\mathrm{part}}`$, to quantify centrality; those numbers for different bins are taken from kn , which agree well with those given by PHENIX ad . Next, we define, for fixed $`\sqrt{s}`$ (at 200 GeV),

$$z=p_T/K(N),$$
(1)
where $`K`$ depends on $`N_{\mathrm{part}}`$, for which we use the abbreviated notation $`N=N_{\mathrm{part}}`$ hereafter. For every centrality bin we vary $`K`$ by plotting the data of $`(2\pi p_T)^1dN_\pi /dp_T`$ in terms of $`z`$ and adjusting the normalization so that all data points lie on a universal curve. That is, we define

$$\mathrm{\Phi }(z)=A(N)K^2(N)\frac{1}{2\pi p_T}\frac{dN_\pi }{dp_T},$$
(2)
and find $`A(N)`$ and $`K(N)`$ such that $`\mathrm{\Phi }(z)`$ has no explicit dependence on $`N`$. That turns out to be possible, as evidenced by Fig. 1. For clarity we show only 5 bins of centrality in that figure. It is a remarkable property of the centrality dependence of the pion spectra that such a universal scaling distribution exists.

The values of $`K(N)`$ used to obtain the scaling behavior are shown in Fig. 2(a) in units of GeV/c. The dependence of $`K(N)`$ on $`N`$ can be well fitted by

$$K(N)=1.2266.36\times 10^4N,$$
(3)
such that $`K(N_{\mathrm{max}})=1`$ at $`N=N_{\mathrm{max}}=350`$. The effects of the degradation of parton momenta are hidden in this formula. Any change of the overall scale of $`K(N)`$ is trivial and does not affect the scaling behavior that we have found. Although the normalization factor $`A(N)`$ does not have a simple dependence on $`N`$, it turns out to depend simply on the number of binary collisions $`N_c`$. The values of $`A(N_c)`$ needed to achieve the scaling $`\mathrm{\Phi }(z)`$ are shown in Fig. 2(b) in a log-log plot. They can be fitted by

$$A(N_c)=530N_c^{0.9}.$$
(4)
From the tables listed in Refs. kn ; ad , $`N_c`$ and $`N`$ can be related by $`N_c=0.44N^{1.33}`$. Note that the normalization of $`\mathrm{\Phi }(z)`$ is set by the most central bin by choosing $`A(N)=1`$ at $`N=N_{\mathrm{max}}`$. If $`A(N_c)`$ were to behave as $`N_c^1`$, it would suggest that the average multiplicity of pions at midrapidity is proportional to $`N_c`$, which is a variable that measures the number of hard collisions. Thus the factor $`N_c^{0.9}`$ in Eq.(4) is an indication that the centrality dependence of the midrapidity multiplicity scales as $`N_c^{0.9}`$ from the $`pp`$ collisions, revealing the effect of suppression of $`p_T`$ in the nuclear medium.

To fit the scaling curve the $`\pi ^0`$ data are insufficient to give us guidance in the small $`z`$ region, since they do not extend below $`p_T`$=1 GeV/c. For $`0<p_T<1`$ GeV/c, we use the $`\pi ^+`$ data of PHENIX for 0-5% centrality ch  shown in Fig. 1. The combined $`\pi ^0`$ and $`\pi ^+`$ data can be well fitted by

$$\mathrm{\Phi }(z)=1200(z^2+2)^{4.8}(1+25e^{4.5z}),$$
(5)
which is shown by the solid line in Fig. 1. We can check its normalization by evaluating the integral

$$I=_0^{10}𝑑zz\mathrm{\Phi }(z)=46.2=\frac{A(N)}{2\pi }\frac{dN_{\pi ^0}}{d\eta }.$$
(6)
For $`N=200`$, say, this gives $`dN_{\pi ^0}/d\eta =149`$, which compares satisfactorily to $`dN_{\mathrm{ch}}/d\eta /(0.5N)=3.2`$ at the same $`N`$ ad . Since the $`\pi ^\pm `$ data do not extend into the $`p_T>2`$ GeV/c region, we do not consider them for centrality analysis here.

The exponential term in Eq. (5) is mainly to fit the low-$`z`$ data that contain thermodynamical effects. At high $`z`$, $`\mathrm{\Phi }(z)`$ behaves as a power law that represents the effects of hard collisions and jet quenching. For all $`z`$, $`\mathrm{\Phi }(z)`$ is a succinct summary of all dynamical effects for all centralities.

In terms of $`\mathrm{\Phi }(z)`$ it is now possible to have an analytic expression of the inclusive distribution of the pions in $`p_T`$ at midrapidity. For convenience, we shall write it in terms of the momentum fraction $`x`$

$$x=p_T/K_0,$$
(7)
where $`K_0`$ is a fixed scale, beyond which no physics of interest need be of concern here. We set $`K_0=10`$ GeV/c for now, although increasing it later, if necessary, is a simple matter. In view of Eq. (1) we thus have

$$z=x\mathrm{\Lambda }(N),\mathrm{\Lambda }(N)=K_0/K(N).$$
(8)
Converting $`(2\pi z)^1dN_\pi /dz`$ to the $`x`$ variable, we define the corresponding pion distribution to be

$$H(x,N)=A^1(N)\mathrm{\Phi }(x,N),$$
(9)
where $`A(N)=A(N_c(N))`$. To see the evolution of the pion distribution with increasing $`N`$, it is more enlightening to study the normalized distribution, defined by

$$P(x,N)=H(x,N)/_0^1𝑑xxH(x,N),$$
(10)
where the upper limit of integration is set to 1 on the assumption that the contribution from $`p_T>K_0`$ is insignificant. Thus $`P(x,N)`$ is the probability distribution of producing a $`\pi ^0`$ at $`x`$, for which the differential phase space is $`xdx`$ due to the 2D nature of $`\stackrel{}{p}_T`$.

In Fig. 3 we show $`P(x,N)`$ for 4 values of $`N`$. Note how $`P(x,N)`$ decreases at high $`x`$ but increases at low $`x`$, when $`N`$ is increased. That is the behavior we expect when high-$`p_T`$ partons are suppressed, giving rise to low-$`p_T`$ partons. The crossover occurs at around $`x=0.06`$, corresponding to $`p_T=0.6`$ GeV/c.

Such an evolution of the $`x`$-distribution is reminiscent of the evolution of the parton distribution in ln$`Q^2`$ in perturbative QCD. Although no precise relationship between the two has been established, it is known that in the latter case the analytical description is simpler in terms of the moments. Thus let us define the moments

$$P_n(N)=_0^1𝑑xx^{n+1}P(x,N).$$
(11)
From Eqs. (5), (9) and (10) we can calculate the $`N`$ dependencies of $`P_n(N)`$, which are shown in Fig. 4 for $`n=1`$ to 5. Evidently, ln$`P_n(N)`$ can be well approximated by linear dependence on $`N`$, i.e.,

$$\mathrm{ln}P_n(N)=a_nb_nN.$$
(12)
The slope parameters $`b_n`$ are shown in the inset of the same figure. The dependence of $`b_n`$ on $`n`$ is also linear. Thus we may rewrite Eq. (12) as

$$\frac{d}{dN}\mathrm{ln}P_n(N)=\lambda n,\lambda =5.542\times 10^4.$$
(13)
This is a very economical way of describing the degradation property of the pion distribution in terms of one basic parameter $`\lambda `$.

A physical interpretation can readily be given for $`\lambda `$ when we consider $`n=1`$, for which $`P_1(N)=x_N`$, the average $`x`$ at $`N`$. From Eq. (13) we obtain

$$x_N=x_{N_0}\mathrm{exp}[\lambda (NN_0)],$$
(14)
which exhibits explicitly the exponential decrease of $`x_N`$ with increasing $`N`$, a behavior that solidifies our physical notion of what the dense medium does to $`p_T`$. For $`N_0=2`$ and $`N=N_{\mathrm{max}}`$, we get

$$x_{N_{\mathrm{max}}}/x_2=0.825,$$
(15)
which gives a quantitative measure of the degree of degradation. From Eq. (13) it is easy also to show that

$$\frac{d}{dN}\frac{x^n_N}{x_N^n}=0,$$
(16)
where $`x^n_N=P_n(N)`$. Hence, the normalized moments of $`P(x,N)`$ are invariant in $`N`$. That is a clue to another invariant form of the distribution.

Before we examine the implications of that clue, we note that the properties of $`P_n(N)`$ displayed in Fig. 4 and described by Eq. (12) cannot be expected to be valid for arbitrarily large n, since the definition of $`P_n(N)`$ in Eq. (11) puts more weight on the high end of $`x`$ when $`n`$ is large. Our cutoff at $`x=1`$, corresponding to $`p_T=K_0=10`$ GeV/c, is based partly on the lack of data at higher $`p_T`$ and partly on the recognition that the contribution from $`p_T>K_0`$ is unimportant when $`n`$ is not too large. To test the validity of our procedure, we have carried out the analysis for $`K_0=20`$ GeV/c, using the same $`\mathrm{\Phi }(z)`$, and found that Eq. (13) remains to be an excellent approximation of the $`n`$ dependence shown in Fig. 4, and that the value of $`\lambda `$ is larger by just 2%, which is less than the experimental errors. Thus we claim that our analysis is stable under variations of $`K_0`$ so long as we consider $`K_010`$ GeV/c and $`n5`$.

The invariance of the normalized moments in Eq. (16) suggests that we should consider yet another scaling variable

$$u=x/x_N=p_T/p_T_N$$
(17)
for any fixed $`N`$. Let us now define

$$\mathrm{\Psi }(u,N)=x_N^2P(x,N),$$
(18)
whose moments are defined by

$$\mathrm{\Psi }_n(N)=_0^{x_N^1}𝑑uu^{n+1}\mathrm{\Psi }(u,N).$$
(19)
Transforming this integral to an integration over $`x`$, we find that

$$\mathrm{\Psi }_n(N)=x_N^nP_n(N).$$
(20)
It then follows from Eq. (16), that

$$d\mathrm{\Psi }_n(N)/dN=0.$$
(21)
Hence, $`\mathrm{\Psi }_n(N)`$ is independent of $`N`$ and we have a scaling function $`\mathrm{\Psi }_n`$, which in turn implies that $`\mathrm{\Psi }(u)`$ is also independent of $`N`$. Indeed, from Eqs. (9) and (10) we see that (18) can be reexpressed as

$$\mathrm{\Psi }(u)=\mathrm{\Phi }(z(u))/𝑑uu\mathrm{\Phi }(z(u)),$$
(22)
where, by virtue of Eqs. (8) and (17),

$$z(u)=x_N\mathrm{\Lambda }(N)u.$$
(23)
Although $`x_N\mathrm{\Lambda }(N)`$ may appear to depend on $`N`$, it actually is a constant

$$\gamma =x_N\mathrm{\Lambda }(N)=z=\frac{𝑑zz^2\mathrm{\Phi }(z)}{𝑑zz\mathrm{\Phi }(z)}=0.414.$$
(24)
Thus using $`z=\gamma u`$ in Eq. (22), we obtain the scaling function $`\mathrm{\Psi }(u)`$

$$\mathrm{\Psi }(u)=2.1\times 10^4(u^2+11.65)^{4.8}(1+25e^{1.864u}).$$
(25)
In Fig. 5 we show $`\mathrm{\Psi }(u)`$ whose shape differs from that of $`\mathrm{\Phi }(z)`$ at low $`u`$ because of the difference in the power-law violating constants. $`\mathrm{\Psi }(u)`$ is a universal form of all the $`P(x,N)`$ shown in Fig. 3.

The scaling property of $`\mathrm{\Psi }(u)`$ is analogous to the KNO scaling of the multiplicity distributions $`P_m(s)`$ in hadronic collisions for $`\sqrt{s}<200`$ GeV kno . It was found that in terms of the scaling variable $`z=m/m`$, where $`m`$ is the multiplicity, the KNO function $`\psi (z)=mP_m(s)`$ is independent of $`s`$. Here, we find that $`\mathrm{\Psi }(u)`$, defined in Eq. (18), is independent of centrality when the scaling variable, $`u=p_T/p_T`$, is used. As it is with KNO scaling, we have

$$u^n=𝑑uu^{n+1}\mathrm{\Psi }(u)=1,\text{for}n=0,1.$$
(26)
The higher moments are what characterize the scaling function, and perhaps scaling violation at some point.

Recall now that the energy scaling distribution found in hy  is the same as the one in Fig. 1 here, although the scaling factor $`K(s)`$ there is different. Since Eq. (22) is independent of $`K(s)`$ or $`K(N)`$, $`\mathrm{\Psi }(u)`$ is thus also the scaling distribution for any energy.

It should be emphasized that no theoretical models have been used in any part of this investigation. The discovery of a scaling behavior over the whole $`p_T`$ range that has been measured offers a simple form of the $`p_T`$ distribution for dynamical models to describe at any centrality and energy. The scaling distribution provides us with not only a simple picture of the complex $`p_T`$ problem, but also a way of quantifying the degree of degradation of the transverse momentum in the dense medium. More importantly, the mere existence of the scaling behavior presents a phenomenological obstacle to the realization of the theoretical expectation that deconfinement results in an anomalous dependence of the $`p_T`$ distribution on centrality. The distribution, $`\mathrm{\Phi }(z)`$ or $`\mathrm{\Psi }(u)`$, indicates that there is no irregularity suggestive of scaling violation, as $`N`$ is varied over the whole range allowed by the Au+Au collisions.

We are grateful to W.A. Zajc and D. d’Enterria for very helpful comments. This work was supported, in part, by the U. S. Department of Energy under Grant No. DE-FG03-96ER40972.


text:
# Field-Induced Thermal Metal-to-Insulator Transition in Underdoped LSCO

## Abstract

The transport of heat and charge in cuprates was measured in single crystals of La2-xSrxCuO4+δ (LSCO) across the doping phase diagram at low temperatures. In underdoped LSCO, the thermal conductivity is found to decrease with increasing magnetic field in the $`T0`$ limit, in striking contrast to the increase observed in all superconductors, including cuprates at higher doping. In heavily-underdoped LSCO, where superconductivity can be entirely suppressed with an applied magnetic field, we show that a novel thermal metal-to-insulator transition takes place upon going from the superconducting state to the field-induced normal state.

In underdoped La2-xSrxCuO4+δ (LSCO), resistivity measurements have revealed the field-induced normal state to be a charge insulator Ando et al. (1995). On the other hand, the superconducting state of underdoped LSCO is a thermal metal, in the sense that there is a clear $`T`$-linear contribution to the thermal conductivity at $`T0`$ Sutherland et al. ; Takeya et al. (2002). Given that in all superconductors investigated to date (including cuprates) heat transport at $`T0`$ is always seen to increase as one goes from the superconducting state to the field-induced normal state, these two observations point to a violation of the Wiedemann-Franz law in underdoped cuprates. Note that this universal law is violated in the electron-doped cuprate Pr2-xCexCuO4+δ (PCCO) at optimal doping Hill et al. (2001), in that low-temperature heat conduction was found to exceed the expected charge conduction by a factor of approximately two. However, the law is recovered in the strongly overdoped regime Proust et al. (2002).

In this Letter, we show the natural assumption that heat conduction will increase upon going from the superconducting state to the field-induced normal state to be incorrect in underdoped LSCO. Indeed, in the $`T0`$ limit the thermal conductivity decreases in the vortex state and the residual linear term drops to a value below our resolution limit in the field-induced normal state. This argues strongly for a thermally insulating normal state and reveals a novel thermal metal-to-insulator transition. These findings shed new light on the nature of the intriguing state of underdoped cuprates above $`H_{c2}`$, for which several proposals have been put forward recently, including stripe order Carlson et al. , d-density wave order Chakravarty et al. (2001) and a Wigner crystal of d-wave hole pairs Chen et al. (2002). As for the Wiedemann-Franz law, a strict test will prove elusive as transport in both charge and heat channels show insulating behavior (in the normal state).

Measurements of the thermal conductivity ($`\kappa `$) were performed down to 40 mK in magnetic fields up to 13 T on single crystals of La2-xSrxCuO4+δ with Sr doping $`x=0`$, 0.05, 0.06, 0.09, 0.17 and 0.20. All samples were grown in an image furnace using the travelling solvent floating zone technique. Sample $`x=0`$ was annealed at 700 oC in flowing argon gas overnight, in order to set the oxygen content at the stoichiometric O4. Other samples were not annealed after growth. The samples were oriented using Laue diffraction and cut into rectangular samples of typical dimensions of length equal to 1-2 mm, thickness equal to 0.1-0.3 mm and width equal to 0.4-1 mm, with the current flowing along the $`a`$ axis of the tetragonal unit cell. Electrical contacts to the samples were made using silver epoxy cured at 500 oC in flowing O2. (In the undoped sample, argon was used instead of oxygen.) Silver wires were then attached to the contact pads using silver epoxy. The error in the geometric factor resulting from the finite width of the contact pads is less than 15%. The resistively determined transition temperature $`T_c`$ is 5.5, 16, 34 and 33.5 K for the $`x=0.06`$ 0.09, 0.17 and 0.20 samples, respectively, where $`T_c`$ is defined as the temperature where $`\rho =0`$. (Samples with $`x=0`$ and 0.05 are not superconducting.) Thermal conductivity measurements on the same samples in zero field are discussed in Ref. Sutherland et al. , along with details of the measurement technique and data analysis. Resistivity was also measured on these samples using the same contacts.

In Fig. 1 the thermal conductivity is plotted as $`\kappa /T`$ vs. $`T^{\alpha 1}`$, where $`\alpha `$ is a free fitting parameter. This type of plot is used to separate the electronic ($`\kappa _{el}`$) and lattice ($`\kappa _{ph}`$) contributions to $`\kappa `$ by making use of their different power-law temperature dependences in the $`T0`$ limit. In the limit $`T0`$ $`\kappa _{el}`$ is linear in $`T`$ for a $`d`$-wave superconductor on account of nodal quasiparticle excitations Durst and Lee (2000); Graf et al. (1996). Quite generally, a linear contribution to $`\kappa `$ at $`T0`$ is direct evidence for fermionic excitations. The phonon contribution can be modeled as $`\kappa _{ph}T^\alpha `$ for phonons limited to scattering from the boundaries of the sample. Note that this in general differs from the simple $`T^3`$ dependence assumed in previous work, a departure which arises from specular reflection of phonons off smooth crystal faces. Thus, $`\kappa _{el}`$ and $`\kappa _{ph}`$ can be separated by fitting the data at low-temperatures to

$$\frac{\kappa }{T}=\frac{\kappa _o}{T}+BT^{\alpha 1}.$$
(1)
The two distinct contributions are identified in Fig. 1 as the intercept and slope of the curves, respectively, when plotting the data as $`\kappa /T`$ vs. $`T^{\alpha 1}`$. Eq. 1 provides an excellent fit to the data for the underdoped samples. Quite generally for almost any material and sample we have investigated, this fitting procedure is a definite improvement upon the standard approach of fitting to $`\kappa /T=\kappa _o/T+BT^2`$, which requires one to typically only make use of data below 150 mK and often results in an overestimate of $`\kappa _o/T`$ (see Ref. Sutherland et al. ).

Although Eq. 1 works well for the underdoped samples, it does not provide a good description of the $`x=0.20`$ and 0.17 data. For this reason, the data for these two more highly doped samples is plotted simply vs. $`T`$ in Fig. 1. The inapplicability of the fit in these two cases is due to a downturn in the data below about 150 mK, as observable in Figs. 1a and 1b. These downturns, however, only impact upon the $`T`$ dependence and not the $`H`$ dependence of $`\kappa `$. Thus, although the origin of the downturns remains uncertain the field dependence can be understood independently of the $`T`$ dependence, which is left as a topic of future investigation.

Having described our analysis of the data, several observations can be made. Firstly, in zero field (solid circles) the data reproduces the results of Refs. Takeya et al. (2002) and Sutherland et al.  whereby a finite residual linear term in $`\kappa (T)`$ is resolved for $`x0.06`$. This proves the existence of delocalized zero-energy quasiparticles throughout the superconducting region. (Among other implications, this essentially rules out the possibility of a $`d+ix`$ order parameter, where $`x=s`$ or $`d`$; see Sutherland et al. ; Proust et al. (2002).) In other words, the $`d`$-wave superconducting state is a thermal metal at all dopings (see also Sutherland et al. ). Outside the superconducting region, i.e. below $`x=0.055`$, the residual linear term becomes extremely small. The power-law fitting procedure described above to extrapolate to $`T=0`$ yields a value of $`\kappa _0/T=3\mu `$W K-2 cm-1 for both the $`x=0`$ and $`x=0.05`$ samples. Now, even though the power-law fit does provide a good description of the data (see Figs. 1e and 1f), all the way up to 0.4 K, the fact that $`\kappa _0/T`$ is 5 times smaller than the value of $`\kappa /T`$ at the lowest data point (40 mK) means that one has to view the extrapolated value with caution. The conservative position is to assume that the parent compound $`x=0`$ is a heat insulator as well as a charge insulator, to regard this minute linear term of $`3\mu `$W K-2 cm-1 as the resolution limit of our technique (for this series of samples) end (a), and to treat the $`x=0`$ data as our reference (for an insulating state in LSCO crystals). We emphasize that the $`x=0.05`$ sample is no more conductive than the parent compound (see Fig. 2), and hence is also taken to be a thermal insulator. By contrast, the linear term in the $`x=0.06`$ sample (at zero field), of magnitude $`12\mu `$W K-2 cm-1, is clearly above the reference limit (by a factor 4) and is thus unambiguously a thermal metal.

This brings us to the second, and principal observation that $`\kappa `$ decreases with increasing field for the underdoped samples ($`x=0.06`$ and 0.09). In Fig. 2, the field evolution of $`\kappa _{el}/T`$ is shown at fields intermediate between 0 and 13 T, as $`\kappa _o/T`$ for LSCO $`x=0.06`$ and 0.09, and as isotherms of $`\kappa /T`$ at 60 mK for $`x=0.17`$ and 0.20 (because of the previously stated difficulties in extrapolating to $`T`$=0 the data for these two samples). By contrast to the underdoped samples, the electronic heat conductivity in the more highly doped samples ($`x=0.17`$ and 0.20) increases with field, as it does in all known superconductors at $`T0`$ end (b). This increase in $`\kappa `$ at 60 mK is qualitatively consistent with the $`T=0`$ field dependence observed in optimally doped YBCO Chiao et al. (1999), which is roughly $`\sqrt{H}`$ and described by semi-classical models Kübert and Hirschfeld (1998); Vekhter and Houghton (1999). Note that $`\kappa `$ is totally independent of magnetic field in our reference sample ($`x=0`$), as in the $`x=0.05`$ sample. This shows that field dependence is a property of the superconducting state. We can therefore use this criterion to establish that the non-superconducting normal state is reached in the bulk by 11 T in sample $`x=0.06`$. Indeed, as seen in Fig. 1d, a further increase of the field to 13 T causes no further change in $`\kappa `$. This claim is supported by resistivity measurements, shown in Fig. 3a, where the resistive onset of superconductivity is entirely absent for fields of 12 T and above (down to 40 mK). We take this as an additional indication that the field-induced (non-superconducting) normal state has been reached by 13 T at $`x=0.06`$ (in the bulk). A zoom on the $`x=0.06`$ data is shown in Fig. 3b, where we can see that $`\kappa _0/T`$ drops by a factor 4 from $`H=0`$ to $`H=13`$ T, where it reaches a value equal to that of the reference sample, namely $`\kappa _0/T=3\mu `$W K-2 cm-1. We conclude that the field-induced normal state in underdoped LSCO is a thermal insulator. This implies the existence of an unprecedented kind of thermal metal-to-insulator transition. The superconducting state is a thermal metal by virtue of its delocalized nodal quasiparticles, while the field-induced normal state in the same sample is a thermal insulator, with either no fermionic excitations or localized fermionic excitations. The well-known crossover from charge metal to charge insulator, identified as the change from a positive $`d\rho /dT`$ at low $`T`$ near and above optimal doping to a negative $`d\rho /dT`$ in underdoped LSCO Ando et al. (1995), now finds a parallel in the heat sector where the metal-insulator crossover is identified as the change from a positive $`d\kappa /dH`$ at low $`T`$ to a negative $`d\kappa /dH`$ (see Fig. 2).

The fundamental question is: what does this imply for the nature of the field-induced normal state in underdoped cuprates? The answer to this question depends on the role of disorder. Is the “normal” state intrinsically metallic, with fermionic excitations which would be delocalized in the absence of disorder? If so, an explanation must be found for why quasiparticles escape localization in the superconducting state. Alternatively, is the normal state intrinsically insulating, with no fermionic excitations (at low energy)?

A tantalizing possibility is that the thermal metal-to-insulator transition is a signature of a competing order that co-exists with superconductivity (see for example Chakravarty et al. (2001); Demler et al. (2001); Kivelson et al. (2002)). In this picture the application of a field suppresses superconductivity and thereby allows the competing phase to define the ground state excitations. Some candidates for the competing phase are thermal insulators (e.g. antiferromagnetic order, Wigner crystal), while others are thermal metals (staggered flux phase, spin-charge separation). The doping evolution of the field dependence of $`\kappa `$ with increasing doping follows naturally from the decreasing importance of the competing order. For instance by $`x=0.17`$ the ground state may have already undergone a quantum critical transition to a state where the only ordered state is that of the d-wave superconductor, thus accounting for the recovery of the conventional increase in $`\kappa `$ with field. One way to possibly distinguish between the two scenarios for the origin of insulating behavior is to reduce the level of disorder and see whether a thermal metal is uncovered.

In closing we mention that because of the insulating character of both heat and charge conduction in $`x=0.06`$ at 13 T, it is technically not possible to perform a real test of the Wiedemann-Franz law. Consequently, the possibility that the two conductivities cease to be equal as they are in overdoped Tl-2201 Proust et al. (2002) but rather diverge as they do in PCCO Hill et al. (2001) is still an open question.

In summary, we have observed in underdoped LSCO a decrease in thermal conductivity with magnetic field upon going from the superconducting state to the field-induced normal state. We show that this result is due to a novel thermal metal-to-insulator transition. By contrast, for optimally and overdoped LSCO the thermal conductivity increases with field, as it does in all other known superconductors. These results impose clear constraints on models of the field-induced normal state of underdoped cuprates, and are expected to help distinguish between various proposed phases that are postulated to compete with superconductivity in this part of the phase diagram. To this end, it will be important to understand the role of disorder.

###### Acknowledgements.

This work was supported by the Canadian Institute for Advanced Research and NSERC of Canada. We would like to thank H. Dabkowska, B.D. Gaulin, G. Luke and K. Hirota for use of their image furnaces and P. Fournier and W.A. MacFarlane for help with the measurements.

$``$ Present address: Laboratoire National des Champs Magnétiques Pulsés, 31432 Toulouse, France.

$``$ Present address: Department of Physics, Cornell University, Ithaca, NY 14853, USA.

$``$ Present address: Départment de Physique, Université de Sherbrooke, Sherbrooke, Québec, Canada.


text:
# Pareto Law in a Kinetic Model of Market with Random Saving Propensity

## Abstract

We have numerically simulated the ideal-gas models of trading markets, where each agent is identified with a gas molecule and each trading as an elastic or money-conserving two-body collision. Unlike in the ideal gas, we introduce (quenched) saving propensity of the agents, distributed widely between the agents ($`0\lambda <1`$). The system remarkably self-organizes to a critical Pareto distribution of money $`P(m)m^{(\nu +1)}`$ with $`\nu 1`$. We analyse the robustness (universality) of the distribution in the model. We also argue that although the fractional saving ingredient is a bit unnatural one in the context of gas models, our model is the simplest so far, showing self-organized criticality, and combines two century-old distributions: Gibbs (1901) and Pareto (1897) distributions.

Considerable investigations have already been made to study the nature of income or wealth distributions in various economic communities, in particular, in different countries. For more than a hundred years, it is known that the probability distribution $`P(m)`$ for income or wealth of the individuals in the market decreases with the wealth $`m`$ following a power law, known as Pareto law pareto :

$$P(m)m^{(1+\nu )},$$
(1)
where the value of the exponent $`\nu `$ is found to lie between 1 and 2 olids ; fuji ; ls . It is also known that typically less than 10% of the population in any country possesses about 40% of the wealth and follow the above power law. The rest of the low-income group population, in fact the majority, clearly follows a different law, identified very recently to be the Gibbs distribution marjit ; dy2 ; anbkcepjb . Studies on real data show that the high income group indeed follow Pareto law, with $`\nu `$ varying from $`1.6`$ for USA dy2 , to $`1.82.2`$ in Japan fuji . The value of $`\nu `$ thus seem to vary a little from economy to economy.

We have studied here numerically a gas model of a trading market. We have considered the effect of saving propensity of the traders. The saving propensity is assumed to have a randomness. Our observations indicate that Gibbs and Pareto distributions fall in the same category and can appear naturally in the century-old and well-established kinetic theory of gas land : Gibbs distribution for no saving and Pareto distribution for agents with quenched random saving propensity. Our model study also indicates the appearance of self-organized criticality bak  in the simplest model so far, namely in the kinetic theory of gas models, when the stability effect of savings sam  is incorporated.

We consider an ideal-gas model of a closed economic system where total money $`M`$ and total number of agents $`N`$ is fixed. No production or migration occurs and the only economic activity is confined to trading. Each agent $`i`$, individual or corporate, possess money $`m_i(t)`$ at time $`t`$. In any trading, a pair of traders $`i`$ and $`j`$ randomly exchange their money marjit ; anbkcepjb ; dy , such that their total money is (locally) conserved and none end up with negative money ($`m_i(t)0`$, i.e, debt not allowed):

$$m_i(t)+m_j(t)=m_i(t+1)+m_j(t+1);$$
(2)
time ($`t`$) changes by one unit after each trading. The steady-state ($`t\mathrm{}`$) distribution of money is Gibbs one:

$$P(m)=(1/T)\mathrm{exp}(m/T);T=M/N.$$
(3)
Hence, no matter how uniform or justified the initial distribution is, the eventual steady state corresponds to Gibbs distribution where most of the people have got very little money. This follows from the conservation of money and additivity of entropy:

$$P(m_1)P(m_2)=P(m_1+m_2).$$
(4)
This steady state result is quite robust and realistic too! In fact, several variations of the trading, and of the ‘lattice’ (on which the agents can be put and each agent trade with its ‘lattice neighbors’ only), whether compact, fractal or small-world like olids , leaves the distribution unchanged. Some other variations like random sharing of an amount $`2m_2`$ only (not of $`m_1+m_2`$) when $`m_1>m_2`$ (trading at the level of lower economic class in the trade), lead to even drastic situation: all the money in the market drifts to one agent and the rest become truely pauper anijmpc ; hayes .

In any trading, savings come naturally sam . A saving propensity factor $`\lambda `$ is therefore introduced in the same model anbkcepjb  (see dy  for model without savings), where each trader at time $`t`$ saves a fraction $`\lambda `$ of its money $`m_i(t)`$ and trades randomly with the rest:

$$m_i(t+1)=m_i(t)+\mathrm{\Delta }m;m_j(t+1)=m_j(t)\mathrm{\Delta }m$$
(5)
where

$$\mathrm{\Delta }m=(1\lambda )[ϵ\{m_i(t)+m_j(t)\}m_i(t)],$$
(6)
$`ϵ`$ being a random fraction, coming from the stochastic nature of the trading.

The market (non-interacting at $`\lambda =0`$ and $`1`$) becomes ‘interacting’ for any non-vanishing $`\lambda (<1)`$: For fixed $`\lambda `$ (same for all agents), the steady state distribution $`P_f(m)`$ of money is exponentially decaying on both sides with the most-probable money per agent shifting away from $`m=0`$ (for $`\lambda =0`$) to $`M/N`$ as $`\lambda 1`$ (Fig. 1(a)). This self-organizing feature of the market, induced by sheer self-interest of saving by each agent without any global perspective, is quite significant as the fraction of paupers decrease with saving fraction $`\lambda `$ and most people end up with some fraction of the average money in the market (for $`\lambda 1`$, the socialists’ dream is achieved with just people’s self-interest of saving!). Interestingly, self-organisation also occurs in such market models when there is restriction in the commodity market acspbkc . Although this fixed saving propensity does not give yet the Pareto-like power-law distribution, the Markovian nature of the scattering or trading processes (eqn. (4)) is lost and the system becomes co-operative. Indirectly through $`\lambda `$, the agents get to know (start interacting with) each other and the system co-operatively self-organises towards a most-probable distribution ($`m_p0`$).

FIG.1: Steady state money distribution (a) $`P_f(m)`$ for the fixed (same for all agents) $`\lambda `$ model, and (b) $`\stackrel{~}{P_f}(m)`$ for some typical values of $`\lambda `$ in the distributed $`\lambda `$ model. The data is collected from the ensembles with $`N=200`$ agents. The inset in (b) shows the scaling behavior of $`\stackrel{~}{P_f}(m)`$. For all cases, agents start with average money per agent $`M/N=1`$.

In a real society or economy, $`\lambda `$ is a very inhomogeneous parameter: the interest of saving varies from person to person. We move a step closer to the real situation where saving factor $`\lambda `$ is widely distributed within the population. One again follows the same trading rules as before, except that

$$\mathrm{\Delta }m=ϵ(1\lambda _j)m_j(t)(1\lambda _i)(1ϵ)m_i(t)$$
(7)
here; $`\lambda _i`$ and $`\lambda _j`$ being the saving propensities of agents $`i`$ and $`j`$. The agents have fixed (over time) saving propensities, distributed independently, randomly and uniformly (white) within an interval $`0`$ to $`1`$ (see bkcsv  for preliminary results): agent $`i`$ saves a random fraction $`\lambda _i`$ ($`0\lambda _i<1`$) and this $`\lambda _i`$ value is quenched for each agent ($`\lambda _i`$ are independent of trading or $`t`$). Starting with an arbitrary initial (uniform or random) distribution of money among the agents, the market evolves with the tradings. At each time, two agents are randomly selected and the money exchange among them occurs, following the above mentioned scheme. We check for the steady state, by looking at the stability of the money distribution in successive Monte Carlo steps $`t`$. Eventually, after a typical relaxation time ($`10^5`$ for $`N=200`$ and uniformly distributed $`\lambda `$) dependent on $`N`$ and the distribution of $`\lambda `$, the money distribution becomes stationary. After this, we average the money distribution over $`10^3`$ time steps. Finally we take configurational average over $`10^5`$ realizations of the $`\lambda `$ distribution to get the money distribution $`P(m)`$. It is found to follow a strict power-law decay. This decay fits to Pareto law (1) with $`\nu =1.02\pm 0.02`$ (Fig. 2). Note, for finite size $`N`$ of the market, the distribution has a narrow initial growth upto a most-probable value $`m_p`$ after which it falls off with a power-law tail for several decades. As can be seen from the inset of Fig. 2, this Pareto law (with $`\nu 1`$) covers the entire range in $`m`$ of the distribution $`P(m)`$ in the limit $`N\mathrm{}`$. We checked that this power law is extremely robust: apart from the uniform $`\lambda `$ distribution used in the simulations in Fig. 2, we also checked the results for a distribution

$$\rho (\lambda )|\lambda _0\lambda |^\alpha ,\lambda _01,0<\lambda <1,$$
(8)
of quenched $`\lambda `$ values among the agents. The Pareto law with $`\nu =1`$ is universal for all $`\alpha `$. The data in Fig. 2 corresponds to $`\lambda _0=0`$, $`\alpha =0`$. For negative $`\alpha `$ values, however, we get an initial (small $`m`$) Gibbs-like decay in $`P(m)`$ (see Fig. 3).

FIG.2: Steady state money distribution $`P(m)`$ in the model for distributed $`\lambda `$ ($`0\lambda <1`$) for $`N=1000`$ agents. Inset shows that the most probable peak $`m_p`$ shifts towards $`0`$ (indicating the same power law for the entire range of $`m`$) as $`N\mathrm{}`$: results for four typical system sizes $`N=100,200,500,1000`$ are shown. For all cases, agents play with average money per agent $`M/N=1`$.

FIG.3: Steady state money distribution $`P(m)`$ in the model for $`N=100`$ agents with $`\lambda `$ distributed as $`\rho (\lambda )\lambda ^\alpha `$ with different values of $`\alpha `$. The inset shows the region of validity of Pareto law with $`1+\nu =2`$. For all cases, agents play with average money per agent $`M/N=1`$.

In case of uniform distribution of saving propensity $`\lambda `$ ($`0\lambda <1`$), the individual money distribution $`\stackrel{~}{P_f}(m)`$ for agents with any particular $`\lambda `$ value, although differs considerably, remains non-monotonic: similar to that for fixed $`\lambda `$ market with $`m_p(\lambda )`$ shifting with $`\lambda `$ (see Fig. 1). Few subtle points may be noted though: while for fixed $`\lambda `$ the $`m_p(\lambda )`$ were all less than of the order of unity (Fig. 1(a)), for distributed $`\lambda `$ case $`m_p(\lambda )`$ can be considerably larger and can approach to the order of $`N`$ for large $`\lambda `$ (see Fig. 1(b)). The other important difference is in the scaling behavior of $`\stackrel{~}{P_f}(m)`$, as shown in the inset of Fig. 1(b). In the distributed $`\lambda `$ ensemble, $`\stackrel{~}{P_f}(m)`$ appears to have a very simple scaling:

$$\stackrel{~}{P}_f(m)(1\lambda )(m(1\lambda )),$$
(9)
for $`\lambda 1`$, where the scaling function $`(x)`$ has non-monotonic variation in $`x`$. The fixed (same for all agents) $`\lambda `$ income distribution $`P_f(m)`$ do not have any such comparative scaling property. It may be noted that a small difference exists between the ensembles considered in Fig 1(a) and 1(b): while $`mP_f(m)𝑑m=M`$ (independent of $`\lambda `$), $`m\stackrel{~}{P_f}(m)𝑑m`$ is not a constant and infact approaches to order of $`M`$ as $`\lambda 1`$. There is also a marked qualitative difference in fluctuations (see Fig. 4): while for fixed $`\lambda `$, the fluctuations in time (around the most-probable value) in the individuals’ money $`m_i(t)`$ gradually decreases with increasing $`\lambda `$, for quenched distribution of $`\lambda `$, the trend gets reversed (see Fig. 4).

FIG.4: Money of the $`i`$-th trader with time. For fixed $`\lambda `$ for all agents in the market: (a) with $`\lambda =0`$, (b) $`\lambda =0.5`$, (c) $`\lambda =0.9`$. For distributed $`\lambda `$ ($`0\lambda <1`$): for agents with (d) $`\lambda =0.1`$, (e) $`\lambda =0.5`$ and (f) $`\lambda =0.9`$ in the market. All data are for $`N=200`$ ($`M/N=1`$).

FIG.5: Money distribution in cases where saving propensity $`\lambda `$ is distributed uniformly within a range of values: (a) When width of $`\lambda `$ distribution is $`0.5`$, money distribution shows power-law for $`0.5<\lambda <1.0`$; (b) When width of $`\lambda `$ distribution is $`0.2`$, money distribution starts showing power-law when $`0.7<\lambda <0.9`$. Note, the exponent $`\nu 1`$ in all these cases when the power law is seen. All data are for $`N=100`$ ($`M/N=1`$).

We now investigate on the range of distribution of the saving propensities in a certain interval $`a<\lambda _i<b`$, where, $`0<a<b<1`$. For uniform distribution within the range, we observe the appearance of the same power law in the distribution but for a narrower region. As may be seen from Fig. 5, as $`ab`$, the power-law behavior is seen for values $`a`$ or $`b`$ approaching more and more towards unity: For the same width of the interval $`|ba|`$, one gets power-law (with same $`\nu `$) when $`b1`$. This indicates, for fixed $`\lambda `$, $`\lambda =0`$ corresponds to Gibbs distribution, and one gets Pareto law when $`\lambda `$ has got non-zero width of its distribution extending upto $`\lambda =1`$. This of course indicates a crucial role of these high saving propensity agents: the power law behavior is truely valid upto the asymptotic limit if $`\lambda =1`$ is included. Indeed, had we assumed $`\lambda _0=1`$ in (8), the Pareto exponent $`\nu `$ immediately switches over to $`\nu =1+\alpha `$. Of course, $`\lambda _01`$ in (8) leads to the universality of the Pareto distribution with $`\nu =1`$ (independent of $`\lambda _0`$ and $`\alpha `$). Indeed this can be easily rationalised from the scaling behavior (9): $`P(m)_0^1\stackrel{~}{P_f}(m)\rho (\lambda )𝑑\lambda `$ $``$ $`m^2`$ for $`\rho (\lambda )`$ given by (8) and $`m^{(2+\alpha )}`$ if $`\lambda _0=1`$ in (8) (for large $`m`$ values).

FIG.6: Cumulative distribution $`Q(m)=_m^{\mathrm{}}P(m)𝑑m`$ of wealth $`m`$ in USA dy2  in 1997 and Japan fuji  in 2000. Low-income group follow Gibbs law (shaded region) and the rest (about 5%) of the rich population follow Pareto law. The inset shows the cumulative distribution for a model market where the saving propensity of the agents is distributed following (8) with $`\lambda _0=0`$ and $`\alpha =0.7`$. The dotted line (for large $`m`$ values) corresponds to $`\nu =1.0`$.

These model income distributions $`P(m)`$ compare very well with the wealth distributions of various countries: Data suggests Gibbs like distribution in the low-income range dy2  (more than 90% of the population) and Pareto-like in the high-income range fuji  (less than 10% of the population) of various countries (Fig. 6). In fact, we have compared one model simulation of the market with saving propensity of the agents distributed following (8), with $`\lambda _0=0`$ and $`\alpha =0.7`$. This model result is shown in the inset of Fig. 6. The qualitative resemblance of the model income distribution with the real data for Japan and USA in recent years is quite intriguing. In fact, for negative $`\alpha `$ values in (8), the density of traders with low saving propensity is higher and since $`\lambda =0`$ ensemble yields Gibbs-like income distribution (3), we see an initial Gibbs-like distribution which crosses over to Pareto distribution (1) with $`\nu =1.0`$ for large $`m`$ values. The position of the crossover point depends on the magnitude of $`\alpha `$. The important point to note is that any distribution of $`\lambda `$ near $`\lambda =1`$, of finite width, eventually gives Pareto law for large $`m`$ limit. The same kind of crossover behavior (from Gibbs to Pareto) can also be reproduced in a model market of mixed agents where $`\lambda =0`$ for a finite fraction of population and $`\lambda `$ is distributed uniformly over a finite range near $`\lambda =1`$ for the rest of the population.

We also considered annealed randomness in the saving propensity $`\lambda `$: here $`\lambda _i`$ for any agent $`i`$ changes from one value to another within the range $`0\lambda _i<1`$, after each trading. Numerical studies for this annealed model did not show any power law behavior for $`P(m)`$; rather it again becomes exponentially decaying on both sides of a most-probable value.

We have numerically simulated here ideal-gas like models of trading markets, where each agent is identified with a gas molecule and each trading as an elastic or money-conserving two-body collision. Unlike in the ideal gas, we introduce (quenched) saving propensity of the agents, distributed widely between the agents ($`0\lambda <1`$). For quenched random variation of $`\lambda `$ among the agents the system remarkably self-organizes to a critical Pareto distribution (1) of money with $`\nu 1.0`$ (Fig. 2). The exponent is quite robust: for savings distribution $`\rho (\lambda )|\lambda _0\lambda |^\alpha `$, $`\lambda _01`$, one gets the same Pareto law with $`\nu =1`$ (independent of $`\lambda _0`$ or $`\alpha `$). It may be noted that the trading market model we have talked about here has got some apparent limitations. The stochastic nature of trading assumed here in the trading market, through the random fraction $`ϵ`$ in (6), is of course not very straightforward as agents apparently go for trading with some definite purpose (utility maximization of both money and commodity). We are however, looking only at the money transactions between the traders. In this sense, the income distribution we study here essentially corresponds to ‘paper money’, and not the ‘real wealth’. However, even taking money and commodity together, one can argue (see hayes ) for the same stochastic nature of the tradings, due to the absence of ‘just pricing’ and the effects of bargains in the market.

Apart from the intriguing observation that Gibbs (1901) and Pareto (1897) distributions fall in the same category and can appear naturally in the century-old and well-established kinetic theory of gas, that this model study indicates the appearance of self-organized criticality in the simplest (gas) model so far, when the stability effect of savings incorporated, is remarkable.


text:
# References

Anomalous dimensions of Wilson operators
in N=4 SYM Theory.

A. V. Kotikov

Bogoliubov Laboratory of Theoretical Physics
Joint Institute for Nuclear Research
141980 Dubna, Russia

L. N. Lipatov and V. N. Velizhanin 222velizh@thd.pnpi.spb.ru
Theoretical Physics Department
Petersburg Nuclear Physics Institute
Orlova Roscha, Gatchina
188300, St. Petersburg, Russia

Abstract

We present the results of two-loop calculations of the anomalous dimension matrix for the Wilson twist-2 operators in the N=4 Supersymmetric Yang-Mills theory for polarized and unpolarized cases. This matrix can be transformed to a triangle form by the same similarity transformation as in the leading order. The eigenvalues of the anomalous dimension matrix are expressed in terms of an universal function with its argument shifted by integer numbers. In the conclusion we discuss relations between the weak and strong coupling regimes in the framework of the AdS/CFT correspondence.

Parton distributions in QCD satisfy the Balitsky-Fadin-Kuraev-Lipatov (BFKL)  and Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP)  equations. Next-to-leading corrections to the BFKL equation were calculated only recently . It is natural to generalize these equations to the supersymmetric case (see Refs.  and references therein). Indeed, the supersymmetric field theories have a number of amazing properties, such as a cancellation of quadratic divergencies and non-renormalization theorems for interaction terms in the lagrangian. Moreover, the supersymmetry is an excellent technical playground for QCD. For example, the empirically established Dokshitzer relation  among elements of the leading order anomalous dimension matrix in the N=1 supersymmetric limit provides a non-trivial check of results of higher order calculations. Another interesting example is the relation between the BFKL and DGLAP equations in the N=4 Supersymmetric Yang-Mills (SYM) theory . In this model one can obtain the anomalous dimensions of the multiplicatively renormalizable twist-2 operators from the eigenvalues of the BFKL kernel . These operators are certain linear combinations of the Wilson operators appearing in the theoretical description of the deep-inelastic $`ep`$ scattering  (note, that in the N=4 SYM the beta function is zero and the Bjorken scaling for structure functions is strongly violated). Using some assumptions the authors of Ref.  derived also an expression for the universal anomalous dimension for the N=4 model in the two-loop approximation. Moreover, the eigenvalues of the anomalous dimension matrices in the polarized and unpolarized cases were obtained from this universal anomalous dimension by an appropriate integer shift of its argument. In this paper we present the results of direct two-loop calculations of these matrices in the N=4 SYM theory.

Now the anomalous dimensions of twist-2 operators in QCD are known up to two loops both for the unpolarized  and polarized  cases. In the N=4 SYM theory  there are one gluon $`g`$, four Majorana fermions $`q`$, three scalars and three pseudoscalars which can be unified in three complex scalars $`\phi `$. All particles belong to the adjoint representation of the gauge group $`SU(N_c)`$. The transition from QCD to the N=4 SYM theory can be performed if one puts in the final expressions $`C_A=C_F=N_c`$, $`T_f=2N_c`$ (the last substitution follows from the fact, that each gluino $`q_i`$ from four Majorana particles gives a half of the Dirac spinor contribution). Furthermore, one should take into account the diagrams with virtual scalars in the polarized structure functions and the graphs with external scalars in the non-polarized distributions. In the last case the anomalous dimension matrix extends to $`3\times 3`$. Below we calculate the anomalous dimensions of the following gauge-invariant twist-2 operators:

$`𝒪_{\mu _1,\mathrm{},\mu _j}^g`$ $`=`$ $`\widehat{S}G_{\rho \mu _1}D_{\mu _2}D_{\mu _3}\mathrm{}D_{\mu _{j1}}G_{\rho \mu _j},`$ (1)
$`\stackrel{~}{𝒪}_{\mu _1,\mathrm{},\mu _j}^g`$ $`=`$ $`\widehat{S}G_{\rho \mu _1}D_{\mu _2}D_{\mu _3}\mathrm{}D_{\mu _{j1}}\stackrel{~}{G}_{\rho \mu _j},`$ (2)
$`𝒪_{\mu _1,\mathrm{},\mu _j}^q`$ $`=`$ $`\widehat{S}\overline{\mathrm{\Psi }}\gamma _{\mu _1}D_{\mu _2}\mathrm{}D_{\mu _j}\mathrm{\Psi },`$ (3)
$`\stackrel{~}{𝒪}_{\mu _1,\mathrm{},\mu _j}^q`$ $`=`$ $`\widehat{S}\overline{\mathrm{\Psi }}\gamma _5\gamma _{\mu _1}D_{\mu _2}\mathrm{}D_{\mu _j}\mathrm{\Psi },`$ (4)
$`𝒪_{\mu _1,\mathrm{},\mu _j}^\phi `$ $`=`$ $`\widehat{S}\overline{\mathrm{\Phi }}D_{\mu _1}D_{\mu _2}\mathrm{}D_{\mu _j}\mathrm{\Phi },`$ (5)
where $`D_\mu `$ are covariant derivatives; the spinor $`\mathrm{\Psi }`$ and field tensor $`G_{\rho \mu }`$ describe gluinos and gluons, respectively, and $`\mathrm{\Phi }`$ is the complex scalar field appearing in the N=4 supersymmetric model. The symbol $`\widehat{S}`$ implies a symmetrization of the tensor in the Lorenz indices $`\mu _1,\mathrm{},\mu _j`$ and a subtraction of its traces. The anomalous dimension matrices can be written as follows for the unpolarized

$$\gamma _{\mathrm{𝐮𝐧𝐩𝐨𝐥}}=\begin{array}{ccc}\gamma _{gg}& \gamma _{gq}& \gamma _{g\phi }\\ \gamma _{qg}& \gamma _{qq}& \gamma _{q\phi }\\ \gamma _{\phi g}& \gamma _{\phi q}& \gamma _{\phi \phi }\end{array}$$
(6)
and polarized cases

$$\gamma _{\mathrm{𝐩𝐨𝐥}}=\begin{array}{cc}\stackrel{~}{\gamma }_{gg}& \stackrel{~}{\gamma }_{gq}\\ \stackrel{~}{\gamma }_{qg}& \stackrel{~}{\gamma }_{qq}\end{array}.$$
(7)
Note, that in the super-multiplet of twist-2 operators there are also operators with fermion quantum numbers and operators anti-symmetric in two Lorentz indices .

Our approach is similar to that of Refs. . In particular we calculated unrenormalized matrix elements of the partonic operators sandwiched between the scalars, fermion and gluon states 444 For the calculations we used the program DIANA , which calls QGRAF  for the generation of Feynman diagrams, and the package MINCER  for FORM  for the evaluation of two-loop diagrams. and the anomalous dimensions were extracted from the expansion of the matrix elements through the renormalization group coefficients, with the condition, that the renormalized matrix elements satisfy the Callan-Symanzik equations. In our calculations we used the modified minimal subtraction scheme ($`\overline{\mathrm{MS}}`$). Because this scheme violates the supersymmetry, the results were transformed to the dimensional reduction scheme ($`\overline{\mathrm{DR}}`$) , explicitly preserving supersymmetry at least in the two-loop level. For this purpose we used the same procedure as in Ref. . Namely, the difference of two-loop results in $`\overline{\mathrm{MS}}`$ and $`\overline{\mathrm{DR}}`$-schemes was related to the difference of the finite contributions of the corresponding one-loop results.

In the polarized case one needs an appropriate choice for the $`\gamma _5`$-prescription. Our procedure is analogous to that of Ref. , which based on ”reading point” method . To begin with, in each trace of the $`\gamma `$-matrix product, we pushed $`\gamma _5`$ to the right hand side using the property of the trace cyclicity. After that we simplified in a straightforward way the product of $`\gamma `$-matrices leaving the $`\gamma _5`$-matrix untouched and used the relation $`\mathrm{Tr}\gamma _\mu \gamma _\nu \gamma _\rho \gamma _\sigma \gamma _5=4iϵ_{\mu \nu \rho \sigma }`$. Then the integration over the loop momenta in the space-time dimension $`D=42\epsilon `$ was performed and the contraction between two Levi-Civita tensors (second one appearing from the projector) in four dimensions was done. In the end we introduced an additional renormalization constant to restore the anticommutativity of $`\gamma _5`$ with other $`\gamma `$-matrices in an accordance with Ref. .

The final two-loop result for the elements of the anomalous dimension matrix in N=4 SYM theory in the $`\overline{\mathrm{DR}}`$-scheme has the following form (multiplied by $`\alpha _s^2N_c^2/(4\pi )^2`$) in the unpolarized case (for even $`j`$)

$`\gamma _{gg}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{500}{9(j1)}}{\displaystyle \frac{16}{j^3}}+{\displaystyle \frac{72}{j^2}}+{\displaystyle \frac{140}{3j}}+{\displaystyle \frac{24}{(j+1)^2}}{\displaystyle \frac{236}{3(j+1)}}{\displaystyle \frac{16}{(j+2)^3}}+{\displaystyle \frac{176}{3(j+2)^2}}`$
$`+`$ $`{\displaystyle \frac{788}{9(j+2)}}16K(j1)+16K(j)16K(j+1)+16K(j+2)+\widehat{Q}(j),`$
$`\gamma _{gq}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{500}{9(j1)}}{\displaystyle \frac{16}{j^3}}+{\displaystyle \frac{72}{j^2}}+{\displaystyle \frac{140}{3j}}{\displaystyle \frac{22}{3(j+1)}}+{\displaystyle \frac{32}{3(j+2)^2}}+{\displaystyle \frac{152}{9(j+2)}}`$
$``$ $`16K(j1)+16K(j)8K(j+1),`$
$`\gamma _{g\phi }^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{500}{9(j1)}}{\displaystyle \frac{16}{j^3}}+{\displaystyle \frac{72}{j^2}}+{\displaystyle \frac{140}{3j}}{\displaystyle \frac{8}{(j+1)^2}}+{\displaystyle \frac{16}{j+1}}{\displaystyle \frac{16}{3(j+2)^2}}{\displaystyle \frac{64}{9(j+2)}}`$
$``$ $`16K(j1)+16K(j),`$
$`\gamma _{qg}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{320}{9(j1)}}+{\displaystyle \frac{32}{j^3}}{\displaystyle \frac{96}{j^2}}+{\displaystyle \frac{8}{3j}}{\displaystyle \frac{96}{(j+1)^2}}+{\displaystyle \frac{944}{3(j+1)}}+{\displaystyle \frac{64}{(j+2)^3}}{\displaystyle \frac{704}{3(j+2)^2}}`$
$``$ $`{\displaystyle \frac{3152}{9(j+2)}}32K(j)+64K(j+1)64K(j+2),`$
$`\gamma _{qq}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{320}{9(j1)}}+{\displaystyle \frac{32}{j^3}}{\displaystyle \frac{96}{j^2}}+{\displaystyle \frac{8}{3j}}+{\displaystyle \frac{88}{3(j+1)}}{\displaystyle \frac{128}{3(j+2)^2}}{\displaystyle \frac{608}{9(j+2)}}32K(j)`$
$`+`$ $`32K(j+1)+\widehat{Q}(j),`$
$`\gamma _{q\phi }^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{320}{9(j1)}}+{\displaystyle \frac{32}{j^3}}{\displaystyle \frac{96}{j^2}}+{\displaystyle \frac{8}{3j}}+{\displaystyle \frac{32}{(j+1)^2}}{\displaystyle \frac{64}{j+1}}+{\displaystyle \frac{64}{3(j+2)^2}}+{\displaystyle \frac{256}{9(j+2)}}`$
$``$ $`32K(j),`$
$`\gamma _{\phi g}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{64}{3(j1)}}+{\displaystyle \frac{24}{j^2}}{\displaystyle \frac{48}{j}}+{\displaystyle \frac{72}{(j+1)^2}}{\displaystyle \frac{236}{j+1}}{\displaystyle \frac{48}{(j+2)^3}}+{\displaystyle \frac{176}{(j+2)^2}}+{\displaystyle \frac{788}{3(j+2)}}`$
$``$ $`48K(j+1)+48K(j+2),`$
$`\gamma _{\phi q}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{64}{3(j1)}}+{\displaystyle \frac{24}{j^2}}{\displaystyle \frac{48}{j}}{\displaystyle \frac{22}{j+1}}+{\displaystyle \frac{32}{(j+2)^2}}+{\displaystyle \frac{152}{3(j+2)}}24K(j+1),`$
$`\gamma _{\phi \phi }^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{64}{3(j1)}}+{\displaystyle \frac{24}{j^2}}{\displaystyle \frac{48}{j}}{\displaystyle \frac{24}{(j+1)^2}}+{\displaystyle \frac{48}{j+1}}{\displaystyle \frac{16}{(j+2)^2}}{\displaystyle \frac{64}{3(j+2)}}+\widehat{Q}(j)`$
and in the polarized case (for odd $`j`$)

$`\stackrel{~}{\gamma }_{gg}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{32}{j^2}}{\displaystyle \frac{280}{3j}}{\displaystyle \frac{32}{(j+1)^3}}+{\displaystyle \frac{64}{(j+1)^2}}+{\displaystyle \frac{280}{3(j+1)}}32K(j)+32K(j+1)+\widehat{Q}(j),`$
$`\stackrel{~}{\gamma }_{gq}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{16}{j^2}}{\displaystyle \frac{140}{3j}}{\displaystyle \frac{8}{(j+1)^3}}+{\displaystyle \frac{32}{(j+1)^2}}+{\displaystyle \frac{142}{3(j+1)}}16K(j)+8K(j+1),`$
$`\stackrel{~}{\gamma }_{qg}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{64}{j^2}}+{\displaystyle \frac{568}{3j}}+{\displaystyle \frac{64}{(j+1)^3}}{\displaystyle \frac{128}{(j+1)^2}}{\displaystyle \frac{560}{3(j+1)}}+32K(j)64K(j+1),`$
$`\stackrel{~}{\gamma }_{qq}^{(1)}(j)`$ $`=`$ $`{\displaystyle \frac{32}{j^2}}+{\displaystyle \frac{284}{3j}}+{\displaystyle \frac{16}{(j+1)^3}}{\displaystyle \frac{64}{(j+1)^2}}{\displaystyle \frac{284}{3(j+1)}}+16K(j)16K(j+1)+\widehat{Q}(j),`$
where

$`\widehat{Q}(j)={\displaystyle \frac{4}{3}}S_1(j)+16S_1(j)S_2(j)+8S_3(j)8\stackrel{~}{S}_3(j)+16\stackrel{~}{S}_{1,2}(j),`$ (8)
$`K(j)={\displaystyle \frac{1}{j}}\left({\displaystyle \frac{S_1(j)}{j}}+S_2(j)+\stackrel{~}{S}_2(j)\right),`$ (9)
$`S_k(j)={\displaystyle \underset{i=1}{\overset{j}{}}}{\displaystyle \frac{1}{i^k}},`$ (10)
$`\stackrel{~}{S}_k(j)={\displaystyle \underset{i=1}{\overset{j}{}}}{\displaystyle \frac{(1)^i}{i^k}},`$ (11)
$`\stackrel{~}{S}_{k,l}(j)={\displaystyle \underset{i=1}{\overset{j}{}}}{\displaystyle \frac{1}{i^k}}\stackrel{~}{S}_l(i).`$ (12)
The analytical continuation of functions $`\gamma _{ab}^{(1)}(j)`$ ($`a,b=g,q,\phi `$) and $`\stackrel{~}{\gamma }_{ab}^{(1)}(j)`$ ($`a,b=g,q`$) to the complex values of $`j`$ can be done analogously to Refs. . The procedure of the analytic continuation together with a detailed description of our method of calculations will be presented elsewhere.

The eigenvalues of the anomalous dimension matrices are given below

$`\gamma _I^{(1)}(j)=\gamma _+^{(1)}(j)=\widehat{Q}(j2),`$ (13)
$`\gamma _{II}^{(1)}(j)=\gamma _0^{(1)}(j)=\widehat{Q}(j),`$ (14)
$`\gamma _{III}^{(1)}(j)=\gamma _{}^{(1)}(j)=\widehat{Q}(j+2),`$ (15)
$`\gamma _{IV}^{(1)}(j)=\stackrel{~}{\gamma }_+^{(1)}(j)=\widehat{Q}(j1),`$ (16)
$`\gamma _V^{(1)}(j)=\stackrel{~}{\gamma }_{}^{(1)}(j)=\widehat{Q}(j+1).`$ (17)
In fact they coincide with the expressions predicted in Ref. . Indeed, using the two-loop result

$`\gamma _+(j)=\stackrel{~}{\gamma }_+(j1)=\gamma _0(j2)=\stackrel{~}{\gamma }_{}(j3)=\gamma _{}(j4)`$
$`=\gamma (j)={\displaystyle \frac{\alpha _sN_c}{\pi }}S_1(j2)+({\displaystyle \frac{\alpha _sN_c}{4\pi }})^2\widehat{Q}(j2)`$ (18)
for the universal anomalous dimension $`\gamma (j)`$ we can redefine $`\alpha _s\alpha _s(1\alpha _sN_c/(12\pi ))`$ to remove in $`\widehat{Q}(j)`$ the term proportional to $`S_1(j)`$ (note, however, that it is an additional redefinition of $`\alpha _s`$ in comparison with the transition from $`\overline{\mathrm{MS}}`$ to $`\overline{\mathrm{DR}}`$-scheme). After this substitution the above universal function $`\widehat{Q}(j)`$ in two loops coincides with $`16Q(j)`$ from Ref. .

For the polarized case the Dokshitzer relation is similar to original one (below $`\gamma _{ab}^{(1)}(j)=\gamma _{ab}`$ and $`\stackrel{~}{\gamma }_{ab}^{(1)}(j)=\stackrel{~}{\gamma }_{ab}`$)

$$\stackrel{~}{\gamma }_{gg}+\frac{1}{2}\stackrel{~}{\gamma }_{qg}=\stackrel{~}{\gamma }_{qq}+2\stackrel{~}{\gamma }_{gq}$$
(19)
and we can find that

$`\stackrel{~}{\gamma }_{gg}+{\displaystyle \frac{1}{2}}\stackrel{~}{\gamma }_{qg}`$ $`=`$ $`\widehat{Q}(j1),`$ (20)
$`\stackrel{~}{\gamma }_{gg}2\stackrel{~}{\gamma }_{gq}`$ $`=`$ $`\widehat{Q}(j+1).`$ (21)
There are three relations for the unpolarized case

$`\gamma _{gg}+\gamma _{qg}+\gamma _{sg}=\gamma _{gq}+\gamma _{qq}+\gamma _{sq}=\gamma _{gs}+\gamma _{qs}+\gamma _{ss},`$ (22)
$`\gamma _{gg}4\gamma _{gq}+3\gamma _{gs}={\displaystyle \frac{\gamma _{qg}}{4}}+\gamma _{qq}{\displaystyle \frac{3\gamma _{qs}}{4}}={\displaystyle \frac{\gamma _{sg}}{3}}{\displaystyle \frac{4\gamma _{sq}}{3}}+\gamma _{ss},`$ (23)
$`12\gamma _{gq}12\gamma _{gs}+3\gamma _{qg}3\gamma _{qs}+4\gamma _{sg}4\gamma _{sq}=0`$ (24)
and one can verify that

$`\gamma _{gg}+\gamma _{qg}+\gamma _{sg}`$ $`=`$ $`\widehat{Q}(j2),`$ (25)
$`\gamma _{gg}4\gamma _{gq}+3\gamma _{gs}`$ $`=`$ $`\widehat{Q}(j),`$ (26)
$`\gamma _{gg}\gamma _{gq}{\displaystyle \frac{\gamma _{sg}}{3}}+{\displaystyle \frac{\gamma _{sq}}{3}}`$ $`=`$ $`\widehat{Q}(j+2).`$ (27)
A complete diagonalization of the above anomalous dimension matrices corresponds to the use of a slightly modified basis of the multiplicatively renormalized twist-2 operators in comparison with the leading order  due to the breakdown of the superconformal invariance (cf. ). But Eqs. (19)-(27) are correct in all orders with the replacement of $`\widehat{Q}(j2)`$ by the exact universal dimension $`\gamma (j)`$.

Following the analysis in Ref.  the $`\alpha _s^3`$ correction to the universal anomalous dimensions $`\gamma (j)`$ will be constructed as soon as the corresponding QCD anomalous dimensions will be calculated (see the recent papers  and references therein).

Recently there was a great progress in the investigation of the N=4 SYM theory in a framework of the AdS/CFT correspondence  where the strong-coupling limit $`\alpha _sN_c\mathrm{}`$ is described by a classical supergravity in the anti-de Sitter space $`AdS_5\times S^5`$. In particular, a very interesting prediction  (see also ) was obtained for the large-$`j`$ behavior of the anomalous dimension for twist-2 operators

$$\gamma (j)=a(z)\mathrm{ln}j,z=\frac{\alpha _sN_c}{\pi }$$
(28)
in the strong coupling regime (see Ref.  for asymptotic corrections):

$$\underset{z\mathrm{}}{lim}a=z^{1/2}+\frac{3\mathrm{ln}2}{8\pi }+𝒪\left(z^{1/2}\right).$$
(29)
Here we took into account, that in our normalization $`\gamma (j)`$ contains the extra factor $`1/2`$ in comparison with that in Ref. .

On the other hand, all anomalous dimensions $`\gamma _i(j)`$ and $`\stackrel{~}{\gamma }_i(j)`$ ($`i=+,0,`$) coincide at large $`j`$ and our results for $`\gamma (j)`$ allow one to find two first terms of the small-$`z`$ expansion of the coefficient $`a(z)`$

$$\underset{z0}{lim}\stackrel{~}{a}=z+\frac{\pi ^21}{12}z^2+\mathrm{}.$$
(30)
To go from this expansion to the strong coupling regime we perform a resummation of the perturbative result using a method similar to the Pade approximation and taking into account, that for large $`N_c`$ the perturbation series has a finite radius of convergency. Namely, we present $`\stackrel{~}{a}`$ as a solution of the simple algebraic equation

$$z=\stackrel{~}{a}+\frac{\pi ^21}{12}\stackrel{~}{a}^2.$$
(31)
From this equation the following large-$`z`$ behaviour of $`\stackrel{~}{a}`$ is obtained:

$$\stackrel{~}{a}1.1632z^{1/2}+0.67647+𝒪\left(z^{1/2}\right)$$
(32)
in a rather good agreement with Eq. (29) based on the AdS/CFT correspondence. Note, that if we write for $`\stackrel{~}{a}`$ the more general equation

$$z^n=\underset{r=n}{\overset{2n}{}}C_r\stackrel{~}{a}^r,$$
(33)
the coefficients $`C_r`$ for $`n3`$ can be chosen in such way to include all known information about $`a`$. For $`n=2`$ we can impose on $`\stackrel{~}{a}`$ apart from a correct small-$`z`$ expansion also the condition $`a=z^{1/2}`$ for large-$`z`$ and as a result the asymptotic correction $`𝒪\left(1\right)`$ will be in a better agreement with Eq. (29).

Further, for $`j2`$ due to the energy-momentum conservation

$$\gamma (j)=(j2)\gamma ^{}(2)+\mathrm{},$$
(34)
where the coefficient $`\gamma ^{}(2)`$ can be calculated from our results in two first orders of the perturbation theory:

$$\gamma ^{}(2)=\frac{\pi ^2}{6}z+1.2158z^2+\mathrm{}.$$
(35)
Using the same method of resummation as we used above for $`\stackrel{~}{a}`$, we obtain for large $`z`$

$$\gamma ^{}(2)=0.9071z^{1/2}+0.6768+\mathrm{}.$$
(36)
Let us take into account, that in this limit $`\gamma =1/2+i\nu +(j1)/21+(j2)/2`$ for the principal series of unitary representations of the Möbius group appearing in the BFKL equation . Therefore we obtain for large $`z`$

$$j=21.1024z^{1/2}0.2148z^1.$$
(37)
in an agreement with the result, that the Pomeron in the strong coupling regime coincides with the graviton . The correction $`z^{1/2}`$ to the graviton spin $`j=2`$ coincides in form with that obtained in Ref.  from the AdS/CFT correspondence and the coefficient in front of $`z^{1/2}`$ was not calculated yet. Note, that for the soft Pomeron the correction is $`z^1`$ .

One can attempt to calculate the intercept of the BFKL Pomeron also using its perturbative expansion in Ref.

$$j1=2.7726z5.0238z^2.$$
(38)
After the Pade resummation we obtain in the strong coupling regime $`j2.53010.8444z^1`$ in an reasonable agreement with the AdS/CFT estimate (see ). Note, however, that in the upper orders of the perturbation theory the BFKL equation should be modified by including the contributions from multi-gluon components of the Pomeron wave function.

In the conclusion we want to stress again, that the AdS/CFT correspondence unified with a resummation procedure gives a possibility to relate weak and strong coupling results.
Acknowledgments.
We thank Yu. Makeenko for the helpful discussions. This work was supported by the Alexander von Humboldt Foundation fellowship (A.V.K.), the RFBR grants 00-15-96610, 02-02-17513 and INTAS grant 00-366.


text:
# Primes in the denominators of Igusa class polynomials

## 1. Introduction

The purpose of this note is to suggest an analogue for genus $`2`$ curves of part of Gross and Zagier’s work on elliptic curves \[GZ84\]. Experimentally, for genus $`2`$ curves with CM by a quartic CM field $`K`$, it appears that primes dividing the denominators of the discriminants of the Igusa class polynomials all have the property 1) that they are bounded by $`d`$, the absolute value of the discriminant of $`K`$, and 2) that they divide $`dx^2`$, for some integer $`x`$ whose square is less than $`d`$. A slightly stronger condition is given in Section 3. Such primes are primes of bad reduction for the genus $`2`$ curve and primes of supersingular reduction for the Jacobian of the genus $`2`$ curve.

## 2. Algorithm to generate Igusa class polynomials of quartic CM fields

The steps of this algorithm are as follows:

1. Choose a CM field $`K`$ of degree $`4`$ which is either not Galois or Galois with Galois group $`/4`$. Let $`K_0`$ be the totally real quadratic subfield. Choose a CM-type, $`\mathrm{\Phi }`$ (a choice of two complex embeddings of $`K`$, such that one is not the complex conjugate of the other).
2. For each element of the ideal class group of $`𝒪_K`$, choose a representative ideal $`𝔞_i`$, and find an integral basis $`\{1,\tau _i\}`$ for it over the ring of integers of $`K_0`$.
3. Each ideal class corresponds to an isomorphism class of an abelian variety via Shimura’s theory \[Sh94, p.126\]. The principal polarization(s) on the abelian variety and the corresponding $`2\times 2`$ period matrix (matrices) $`\mathrm{\Omega }_i`$ are given in \[Sp94, Section 4.2 and p.62\]. The entries of $`\mathrm{\Omega }_i`$ are given in terms of $`\tau _i`$, the CM-type $`\mathrm{\Phi }`$, and $`\omega `$, the generator of the ring of integers of $`K_0`$.
4. For each period matrix, evaluate the ten even theta constants up to some amount of precision.
5. The three absolute Igusa invariants, $`j_1`$, $`j_2`$, $`j_3`$, associated to each period matrix are defined as combinations of the ten even theta constants by formulas given in the appendix.
6. Take the product over all possible period matrices to form the three Igusa class polynomials:

$$h_1(x)=\underset{\mathrm{\Omega }_i}{}(xj_1(\mathrm{\Omega }_i)),$$

$$h_2(x)=\underset{\mathrm{\Omega }_i}{}(xj_2(\mathrm{\Omega }_i)),$$

$$h_3(x)=\underset{\mathrm{\Omega }_i}{}(xj_3(\mathrm{\Omega }_i)).$$

This algorithm was implemented using the software packages Pari and MAPLE in \[CL01\]. Implementations can be found in the literature in \[Sp94\], \[vW99\], \[We03\].

###### Remark 1.

The above algorithm is analogous to computing the Hilbert class polynomial associated to an imaginary quadratic field $`K`$. The $`j`$-invariant for each ideal class is evaluated up to some amount of precision, then $`H(x)`$ is formed by taking the product over all ideal classes of $`(xj(\tau _i))`$, where $`\tau _i`$ is a particular algebraic integer associated to each ideal class. Then the Hilbert class polynomial, $`H(x)`$, has integer coefficients.

The case of a quartic CM field is different because there are three class polynomials, their coefficients are rational numbers, and the amount of precision required to compute them is not known in advance. The next section proposes a constraint on the primes appearing in the denominators which, together with a bound on the power to which each prime appears, would give a bound on the amount of precision required to compute the coefficients of the class polynomials.

## 3. Conjectural formula for the primes in the denominators of the Igusa class polynomials

Let $`h_1`$, $`h_2`$, and $`h_3`$ be the three Igusa class polynomials with rational coefficients obtained from the algorithm in Section 2. Let $`\{q_i\}`$ be the primes appearing in the denominators of the discriminants of all three $`h_i`$. Let $`d`$ be the absolute value of the discriminant of the number field $`K`$ and let $`d_0`$ be the absolute value of the norm of the discriminant of $`K/K_0`$.

Experimentally I have observed that the primes $`q_i`$ satisfy the following property:

Property (1) Each $`q_i`$ divides $`dx^2`$, for some $`x`$, an integer such that $`x^2d`$. In fact in every case so far, each $`q_i`$ divides $`d_0x^2`$, for some $`x`$, an integer such that $`x^2d_0`$.

###### Remark 2.

I have only tested this property on a few handfuls of examples, but Annegret Weng has also tested it on at least that many others. I am currently working on a proof of this property jointly with Annegret Weng, Farshid Hajir, Fernando Rodriguez-Villegas, and Tonghai Yang.

## 4. Appendix: Igusa invariants

Let $`\{\theta _i\}_{i=1,10}`$ be the ten even theta constants associated to a given period matrix. Define the functions $`f`$ and $`g`$ as follows:

$$f(k_1,k_2,k_3,k_4,k_5,k_6,k_7,k_8)=(\underset{i=1,\mathrm{},8}{}\theta _{k_i})^4,$$

$$g(k_1,k_2,k_3,k_4,k_5,k_6)=(\underset{i=1,\mathrm{},6}{}\theta _{k_i})^4.$$

Now $`h_4`$, $`h_{10}`$, $`h_{12}`$, and $`h_{16}`$ are values of modular forms of weights $`4`$, $`10`$, $`12`$, and $`16`$:

$`h_4=_{i=1,\mathrm{},10}\theta _i^8`$,

$`h_{10}=_{i=1,\mathrm{},10}\theta _i^2`$,

$`h_{12}=g(1,5,2,9,6,10)+g(1,2,9,6,8,3)+g(5,9,6,8,10,7)+g(5,2,6,8,3,7)+g(1,5,2,10,3,7)+g(1,9,8,10,3,7)+g(1,5,2,8,10,4)+g(1,5,9,8,3,4)+g(5,9,6,10,3,4)+g(2,6,8,10,3,4)+g(1,2,9,6,7,4)+g(1,5,6,8,7,4)+g(2,9,8,10,7,4)+g(5,2,9,3,7,4)+g(1,6,10,3,7,4),`$

$`h_{16}=f(8,1,5,2,9,6,8,10)+f(5,1,5,2,9,6,8,3)+f(10,1,2,9,6,8,10,3)+f(3,1,5,2,9,6,10,3)+f(1,1,5,9,6,8,10,7)+f(2,5,2,9,6,8,10,7)+f(1,1,5,2,6,8,3,7)+f(9,5,2,9,6,8,3,7)+f(9,1,5,2,9,10,3,7)+f(6,1,5,2,6,10,3,7)+f(5,1,5,9,8,10,3,7)+f(2,1,2,9,8,10,3,7)+f(6,1,9,6,8,10,3,7)+f(8,1,5,2,8,10,3,7)+f(10,5,2,6,8,10,3,7)+f(3,5,9,6,8,10,3,7)+f(7,1,5,2,9,6,10,7)+f(7,1,2,9,6,8,3,7)+f(9,1,5,2,9,8,10,4)+f(6,1,5,2,6,8,10,4)+f(2,1,5,2,9,8,3,4)+f(6,1,5,9,6,8,3,4)+f(1,1,5,9,6,10,3,4)+f(2,5,2,9,6,10,3,4)+f(1,1,2,6,8,10,3,4)+f(5,5,2,6,8,10,3,4)+f(9,2,9,6,8,10,3,4)+f(8,5,9,6,8,10,3,4)+f(10,1,5,9,8,10,3,4)+f(3,1,5,2,8,10,3,4)+f(5,1,5,2,9,6,7,4)+f(2,1,5,2,6,8,7,4)+f(9,1,5,9,6,8,7,4)+f(8,1,2,9,6,8,7,4)+f(1,1,2,9,8,10,7,4)+f(5,5,2,9,8,10,7,4)+f(6,2,9,6,8,10,7,4)+f(10,1,2,9,6,10,7,4)+f(10,1,5,6,8,10,7,4)+f(1,1,5,2,9,3,7,4)+f(6,5,2,9,6,3,7,4)+f(8,5,2,9,8,3,7,4)+f(5,1,5,6,10,3,7,4)+f(2,1,2,6,10,3,7,4)+f(9,1,9,6,10,3,7,4)+f(8,1,6,8,10,3,7,4)+f(10,5,2,9,10,3,7,4)+f(3,1,2,9,6,3,7,4)+f(3,1,5,6,8,3,7,4)+f(3,2,9,8,10,3,7,4)+f(7,1,5,2,8,10,7,4)+f(7,1,5,9,8,3,7,4)+f(7,5,9,6,10,3,7,4)+f(7,2,6,8,10,3,7,4)+f(4,1,5,2,9,6,10,4)+f(4,1,2,9,6,8,3,4)+f(4,5,9,6,8,10,7,4)+f(4,5,2,6,8,3,7,4)+f(4,1,5,2,10,3,7,4)+f(4,1,9,8,10,3,7,4).`$

The Igusa invariants are

$$I_2=\frac{h_{12}}{h_{10}},$$

$$I_4=h_4,$$

$$I_6=\frac{h_{16}}{h_{10}},$$

$$I_{10}=h_{10},$$

and the absolute Igusa invariants are

$$j_1=\frac{I_2^5}{I_{10}},$$

$$j_2=\frac{I_4I_2^3}{I_{10}},$$

$$j_3=\frac{I_6I_2^2}{I_{10}}.$$