text:
# Chapter 1 VHE Observations of BL Lacertae Objects: 1995-2000

## Chapter 1 VHE Observations of BL Lacertae Objects: 1995-2000

### Abstract

The results of observations of 29 BL Lacertae objects taken with the Whipple Observatory 10 m gamma-ray Telescope between 1995 and 2000 are presented.

### 1. Introduction

Among blazars, BL Lacertae objects (BL Lacs) are believed to be the best candidates for VHE emission. Consistent with these expectations, the four confirmed sources of extragalactic gamma rays (E $`>`$ 250 GeV) are BL Lacs, Mrk 421 , Mrk 501 , H1426  and 1ES1959 . To improve our understanding of the emission mechanisms in BL Lacs, more need to be detected at very high energies.

### 2. Observations and Analysis

The observations presented here were taken as part of three BL Lac campaigns. The first was a survey of all known (circa 1995) BL Lacs with redshift $`<0.1`$. The second was a search for TeV emission from high frequency peaked BL Lacs in the redshift range from 0.1 to 0.2. More recent observations were taken as part of a “Snapshot Survey” in which many BL Lacs were observed for 10 minutes each on a regular basis in the hope of catching one of them in a flaring state .

Throughout the course of these observations, the imaging camera on the Whipple telescope was upgraded a number of times, the triggering criteria were changed and different light concentrators were installed in front of the camera. Therefore, the peak response energy of the instrument varied during this period. These observations were all taken in the “TRACKING mode” described in .

### 3. Results And Discussion

Tables 1 & 2 summarise the results of the BL Lac observations. No significant excesses were detected from any of the objects on time scales of days, months or years. The mean exposure on each object was 5.5 hours for each season that it was observed. Typically, in order to detect a signal at the 4$`\sigma `$ level during a 5 hour exposure, the object would need to have had a flux of at least 0.4 times that of the Crab above the peak response energy for that season.

Flux upper limits were calculated for each object for each season. Costamante & Ghisellini  have made predictions for the TeV flux from fourteen of the BL Lacs included in this paper using two different methods. The upper limits presented here were compared with these predictions; those of six BL Lacs, highlighted in Tables 1 & 2, were found, during a number of observing seasons, to be lower than the predicted fluxes according to the Fossati approach  adapted in . It should be noted however, that the upper limits quoted here pertain only to the specific period during which the observations were made. Indeed, 1ES2344, H1426 and 1ES1959 were initially observed as part of this BL Lac campaign and, like the objects listed here, were not detected. In subsequent years, continuous monitoring, similar to that described here, revealed these objects to be TeV emitters when in more active states.

The analysis of these results is ongoing. Spectral energy distributions are being constructed and the flux predictions are being corrected to account for pair-production of the gamma rays with the infra red background radiation so that our upper limits can be compared to these predictions in a more meaningful way.

### 4. References

1. Catanese M. A. et al. 1998, ApJ 501, 616

2. Costamante L. & Ghisellini G. 2002, A&A 384, 56

3. D’Vali M. et al. 1999, in Proc of 26th ICRC, Salt Lake City, 3, 422

4. Fossati G. et al. 1998, MNRAS 299, 433

5. Horan D. et al. 2001, AIP Conf. Proc. 587, 324

6. Nishiyama T. et al. 2000, AIP Conf. Proc. 516, 26, 369

7. Punch M. et al. 1992, Nature 358, 477

8. Quinn J. et al. 1998, ApJ 456, L83

9. Vernon-Cetty M. P. & Vernon P. 1993, A&AS 100, 521


format: text

text:
# Chapter 1 Measurement of Cosmic-Ray Proton, Antiproton and Muon Spectra at Mountain Altitude

## Chapter 1 Measurement of Cosmic-Ray Proton, Antiproton and Muon Spectra at Mountain Altitude

### Abstract

Measurement of cosmic-ray proton, antiproton and muon spectra was carried out at mountain altitude. We observed $`2\times 10^5`$ protons and $`10^2`$ antiprotons in a kinetic energy region of 0.25 – 3.3 GeV. Zenith-angle dependence of proton fluxes was obtained. Atmospheric muon spectra were measured simultaneously. The observed antiproton spectrum showed some deviation from theoretical predictions particularly in a low energy region.

### 1. Introduction

Primary cosmic rays hit the Earth’s atmosphere and produce baryons and mesons via hadronic interactions. Absolute fluxes of these “secondary cosmic rays” can be calculated by using the primary cosmic-ray intensity and interaction cross sections. Observation of the secondary cosmic rays is very important to verify, or to improve, theoretical calculations. It is essentially important to understand propagation process of the secondary particles inside the atmosphere. We report new measurement of secondary cosmic ray spectra at mountain altitude.

### 2. Observations

We performed cosmic-ray observation at Norikura Observatory, ICRR, University of Tokyo, Japan, in September 1999, with the BESS detector , which was the same apparatus as we utilized to measure the primary protons and antiprotons , as well as atmospheric muons at sea level . The observatory is located at 2,770m above sea level. The vertical cutoff rigidities is 11.2 GV . During the observation, the mean atmospheric depth was 742 $`\mathrm{g}/\mathrm{cm}^2`$.

### 3. Results and Discussion

The measured energy spectra of protons and antiprotons are shown in Fig. 1, together with the previous measurements at mountain altitude . The antiproton spectrum is compared with theoretical predictions . In this analysis, the zenith angle $`(\theta _z)`$ was limited as $`\mathrm{cos}\theta _z0.95`$ for protons, and $`\mathrm{cos}\theta _z0.84`$ for antiprotons, thus the obtained fluxes are “near-vertical” fluxes.

There is some disagreement among the proton fluxes shown in Fig. 1. According to simple Monte Carlo simulations, the deviations can be explained by the different altitudes and cutoff rigidities at their observation sites.

Fig. 2 shows zenith angle dependence of the observed proton flux in two kinetic energy regions.

The zenith angle dependence can be expected in a simple one-dimensional approximation, $`F(\mathrm{cos}\theta _z)=F_0\mathrm{exp}(X/\lambda (11/\mathrm{cos}\theta _z))`$, where $`X`$ is the atmospheric depth and $`\lambda `$ is the absorption mean free path of protons inside the atmosphere. Dotted lines in Fig. 2 show the expectation, in which $`X/\lambda =6`$ was assumed. Relatively good agreement was found between the observed data and the calculation in the higher energy region. In the lower energy region, however, significant discrepancy was found between them. These facts are most likely due to the effect of angular spread of secondary protons produced via nuclear interactions. The solid lines shown in Fig. 2 give the results of the analytic calculation, in which the angular spread was taken into account. They reproduced the observed data better than one-dimensional calculation in the whole energy range.

The observed antiproton spectrum agrees with a theoretical calculation by Stephens  above 1 GeV as shown in Fig. 1. On the other hand, the flux below 1 GeV shows significant disagreement. In Ref. , production spectra of antiprotons are calculated, and the result shows a sharp peak around 2 GeV. This means that most antiprotons observed below 1 GeV are tertiary antiprotons, those which have been produced inside the atmosphere and then lost their energies during the propagations in the atmosphere. In this case, cross sections in $`\overline{p}+A(nuclei)`$ processes are to be precisely treated for an accurate evaluation of antiproton spectrum at mountain altitude. A recent work made with Monte Carlo simulation reported a preliminary result which shows better agreement with the observed spectrum .

The atmospheric muon spectra in a momentum range of 0.6 – 106 GeV/$`c`$ have been measured simultaneously. The results are reported in Ref. . The observed muon spectra showed much better agreement with theoretical prediction than that of antiprotons. It shows that secondary meson productions are treated rather properly in the theoretical calculations. However, secondary baryon interaction cross sections have to be modified so as to reproduce the observed antiproton spectrum.

### 4. Summary

We have measured proton, antiproton and muon spectra at Mt. Norikura, Japan, where the atmospheric depth was 742 $`\mathrm{g}/\mathrm{cm}^2`$. The zenith angle dependence in the proton flux was observed. It suggests an importance of three-dimensional effect of angular spread in secondary baryon productions. The calculated antiproton flux agrees with our antiproton flux above 1 GeV. In a lower energy region, however, our measurement gives much lower flux than that of the prediction.

### Acknowledgment

This study was supported by the Joint Research Program of ICRR, the University of Tokyo and Grants-in-Aid, KAKENHI(11694104, 11440085, 09304033), from MEXT and JSPS. We would like to thank NASA, ISAS, KEK and ICEPP, the University of Tokyo for their continuous support.

a Present address: Kanagawa University, Yokohama, Kanagawa 221-8686, Japan

b Present address: CERN, CH-1211 Geneva 23, Switzerland

c Present address: University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan

d Present address: NASA/GSFC, Greenbelt, MD 20771, USA

deceased.

1. Ajima Y. et al. 2000, Nucl. Instr. and Meth. A 443, 71

2. Asaoka Y. et al. 1998, Nucl. Instr. and Meth. A 416, 236

3. Asaoka Y. et al. 2002, Phys. Rev. Lett. 88, 051101

4. Barber H. B. et al. 1980, Phys. Rev. D 22, 2667

5. Bowen T. and Moats A. 1986, Phys. Rev. D 22, 2667

6. Buénerd M. 2002, Int. J. Mod. Phys. A 17, 1665

7. Kocharian N. M. 1954, J. Exper. Theoret. Phys. (USSR) 28, 160

8. Kocharian N. M. et al. 1958, J. Exper. Theoret. Phys. (USSR) 35, 1335

9. Maeno T. et al. 2001, Astropart. Phys. 16, 121

10. Motoki M. et al. 2003, Astropart. Phys. 19, 113

11. Orito S. 1987, Proc. ASTROMAG Workshop, KEK Report KEK87-19, 111

12. Orito S. et al. 2000, Phys. Rev. Lett. 84, 1078

13. Sanuki T. et al. 2000, ApJ 545, 1135

14. Sanuki T. et al. 2002, Phys. Lett. B 541, 234

15. Sembroski G. H. et al. 1986, Phys. Rev. D 33, 639

16. Shea M. A. and Smart D. F. 2001, Proc. 27th ICRC(Hamburg), 4063

17. Shikaze Y. et al. 2000, Nucl. Instr. and Meth. A 455, 596

18. Stephens S. A. 1997, Astropart. Phys. 6, 229


text:
# Cold scalar-tensor black holes coupled to a massless scalar field

## Abstract

New four-dimensional black hole solutions of Brans-Dicke equations with a negative constant $`\omega `$, coupled to a massless scalar field, are presented. The temperature of these black holes is zero and the horizon area is infinite. An astrophysical application is also discussed.

Scalar-tensor theories of gravity are considered as the most natural generalizations of general relativity. From theoretical point of view one should mention that specific scalar-tensor theories arise naturally as a low energy limit of string theory. One of the most interesting aspects of scalar-tensor theories is the possible existence of black holes different form those in general relativity. Campanelli and Lousto  pointed out that, in Brans-Dicke theory with a negative coupling constant $`\omega `$, there exist solutions possessing all black holes properties but with horizons of infinite area (type B black holes). Let us note that that Brans-Dicke theory with $`\omega <0`$ is also worth to study for more than its intrinsic theoretical interest. For sufficiently negative $`\omega `$, the theory passes all the experimental constrains as does for possitive $`\omega `$ .

Later on the work of Campanelli and Lousto was generalized in  where it was shown, in the framework of a general class of scalar-tensor theories, that nontrivial black hole solutions can exist for the coupling function $`2\omega (\mathrm{\Phi })+3<0`$. These black holes solutions are divided into two classes: class $`B1`$, where horizons are attained by infalling particles in a finite proper time, and class $`B2`$, for which this proper time is infinite. The structure and stability of the scalar-tensor black holes were investigated in . Charged scalar-tensor black holes were studied in  where it was shown that they exist for anomalous versions of the scalar-tensor theories with a negative kinetic term in the lagrangian. The temperature of these black holes is zero and the horizon area is (in most cases) infinite.

Thermodynamics of the black holes with infinite horizon area was considered in .

Stationary, axisymmetric black holes in Bans-Dicke theory with $`5/4\omega <3/2`$ were discussed in . These black holes are characterize with zero Hawking temperature and finite horizon area.

The aim of the present work is to present new type scalar-tensor black hole solutions of Brans-Dicke equations with $`2\omega +3<0`$, coupled to a real massless scalar field.

The action for the scalar-tensor gravity in the presence of a real massless scalar filed $`\sigma `$ is:

$$S=\frac{1}{16\pi }d^4x\sqrt{g}\left(\mathrm{\Phi }R\frac{\omega }{\mathrm{\Phi }}g^{\mu \nu }_\mu \mathrm{\Phi }_\nu \mathrm{\Phi }2g^{\mu \nu }_\mu \sigma _\nu \sigma \right)$$
(1)
The field equations are given by:

$`R_{\mu \nu }{\displaystyle \frac{1}{2}}g_{\mu \nu }R`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Phi }}}T_{\mu \nu }+{\displaystyle \frac{\omega }{\mathrm{\Phi }^2}}\left(_\mu \mathrm{\Phi }_\nu \mathrm{\Phi }{\displaystyle \frac{1}{2}}g_{\mu \nu }_\alpha \mathrm{\Phi }^\alpha \mathrm{\Phi }\right)`$ (2)
$`+`$ $`{\displaystyle \frac{1}{\mathrm{\Phi }}}\left(_\mu _\nu \mathrm{\Phi }g_{\mu \nu }_\alpha ^\alpha \mathrm{\Phi }\right),`$ (3)
$`_\alpha ^\alpha \mathrm{\Phi }`$ $`=`$ $`{\displaystyle \frac{T}{3+2\omega }},`$ (4)
$`_\alpha ^\alpha \sigma `$ $`=`$ $`0,`$ (5)
where

$$T_{\mu \nu }=2_\mu \sigma _\nu \sigma g_{\mu \nu }_\alpha \sigma ^\alpha \sigma $$
(6)
and $`T=g^{\mu \nu }T_{\mu \nu }`$.

We were able to find111This class of exact solutions can be obtained by using solution generating techniques developed in . the following static, spherically symmetric class of exact solutions of Eq.(2)-Eq.(6) for $`2\omega +3<0`$:

$`ds^2`$ $`=`$ $`\left({\displaystyle \frac{1\beta ^2f^b(r)}{1\beta ^2}}\right)^2\left(f^{ab}(r)dt^2+f^{ab}(r)dr^2+f^{1ab}(r)r^2d\mathrm{\Omega }^2\right),`$ (7)
$`\mathrm{\Phi }^1(r)`$ $`=`$ $`\left({\displaystyle \frac{1\beta ^2f^b(r)}{1\beta ^2}}\right)^2f^b(r),`$ (8)
$`\sigma (r)`$ $`=`$ $`\pm \beta 3+2\omega ^{1/2}{\displaystyle \frac{1f^b(r)}{1\beta ^2f^b(r)}}`$ (9)
where

$$f(r)=1\frac{2\lambda }{r}.$$
(10)
The found class depends on three essential parameters $`\lambda `$, $`\beta `$ ($`\beta ^2<1`$) and $`a`$. The parameters $`b`$ and $`a`$ are related by

$$(2\omega +3)b^2=1a^2$$
(11)
and therefore $`a^21`$. Here we shall consider the solutions with $`\lambda >0`$.

Let us note that the usual vacuum Brans-Dicke solution (with $`2\omega +3<0`$) is recovered for $`\beta =0`$.

Black hole solutions with a metric

$$ds^2=g_{00}(r)dt^2+g_{11}(r)dr^2+g_{22}(r)d\mathrm{\Omega }^2$$
(12)
are selected by the following criteria :

1. There exists a Killing horizon: at some $`r=r_{}`$, $`g_{00}(r_{})=0`$;
2. The integral $`\left(\frac{g_{11}}{g_{00}}\right)^{1/2}𝑑r`$ $``$ $`\mathrm{}`$ as $`rr_{}`$ (invisibility of the horizon for an observer at rest);
3. The Hawking temperaure $`T_H`$ is finite;
4. The invariants $`𝒦_1=R^2`$, $`𝒦_2=R_{\mu \nu }R^{\mu \nu }`$, $`𝒦_3=R_{\mu \nu \alpha \beta }R^{\mu \nu \alpha \beta }`$ are finite at $`r=r_{}`$.

It can be checked that all criteria are satisfied by the solutions with

$$1<a,\mathrm{\hspace{0.33em}\hspace{0.33em}0}<b,\mathrm{\hspace{0.33em}\hspace{0.33em}2}ab<a.$$
(13)
Therefore, these solutions describe scalar-tensor black holes with a regular horizon at $`r=r_{}=2\lambda `$. These black holes are cold since $`T_H=0`$ and have a horizon with an infinite area (i.e. type B black holes). It can be shown that the infalling particles reach the horizon for finite proper time when $`0<b<1`$ (type B1 black holes) and this time is infinite when $`b1`$ (type B2 black holes). It is worth also noting that the scalar field $`\sigma `$ is regular on the horizon. Following the approach presented in  it can be shown that our black holes are characterized with zero entropy ($`𝒮=0`$).

In this work we presented new type black hole solutions in Brans-Dicke theory with a negative coupling parameter $`\omega `$ : black holes coupled to a massless scalar field. These black hole solutions are interesting not only from a pure theoretical point of view. They, as well as the vacuum scalar-tensor black holes, could be of astrophysical relevance and some astrophysical applications were discussed in . Here we would like to point out another astrophysical application. Gravitational lensing of the scalar-tensor black holes in the strong field regime could be quite different from that of the general relativistic black holes and could exhibit some new features. In fact, the gravitational lensing of objects (naked singularities) with a space-time metric similar to that of the scalar-tensor black holes was investigated by Virbhadra and Ellis in . They found that the gravitational lensing of the strongly naked singularity is qualitatively different from that of a Schwarzschild black hole . However, while the metric studied in  is characterized with a parameter $`a<1`$ (in our notations, and $`\nu <1`$ in the notations of ), the metrics of scalar-tensor black holes have $`a>1`$. This could lead to qualitatively new effects. So, the strong field gravitational lensing may be used to distinguish between the scalar-tensor and general relativistic black holes. Of course, this question needs more careful investigation .

Acknowledgments

I would like to thank V. Rizov for reading the manuscript. This work was partially supported by Sofia University Grant No3429.


text:
# Modeling peculiar velocities of dark matter halos This paper was published in Mon. Not. R. Astron. Soc. 343, 1312-1318 (2003). Owing to an error in numerical computations, some incorrect results were presented there. Erratum is to be published in Mon. Not. R. Astron. Soc. Conclusions of the original version are unaffected by the correction. This version supersedes the original version.

## 1 Introduction

In the standard structure formation scenario, dark matter halos are thought to be formed via the gravitational amplification of initial small density fluctuations, and their motion results from gravitational forces acting over large scales. As a consequence, the root-mean-square (rms) peculiar velocities of halos at the present time depend on the total amount of matter in the universe, $`\mathrm{\Omega }_0`$, and on the normalization of the matter power spectrum, $`\sigma _8`$, as $`\sigma _8\mathrm{\Omega }_0^{0.6}`$. In addition, linear theory predicts that the halo peculiar velocity evolves as $`a(t)dD(t)/dt`$ where $`a(t)`$ is the scale factor normalized unity at the present time and $`D(t)`$ is the linear growth rate (see Section 2). Thus, a measurement of the rms peculiar velocities of halos can be used as a test of structure formation scenarios in principle.

A measurement of the peculiar velocity of a dark matter halo is possible through its baryon contents. Because a cluster of galaxies is formed in the gravitational potential of a massive dark matter halo, its peculiar motion traces that of the hosting halo. Therefore, the rms peculiar velocity of clusters of galaxies provides a reasonable measure of that of halos. So far, empirical distance indicators such as the Tully–Fisher and the $`D_n`$$`\sigma `$ relations have been used to measure the peculiar velocity of relatively nearby systems (e.g. Bahcall, Gramann & Cen 1994; Lauer & Postman 1994; Bahcall & Oh 1996; Moscardini et al. 1996; Borgani et al. 1997; 2000; Watkins 1997; Dale et al. 1999; Hudson et al. 1999; Colless et al. 2001). Future observations of the kinematic Sunyaev–Zel’dovich effect will provide measurements of line-of-sight peculiar velocity for distant systems (Sunyaev & Zel’dovich 1980; Rephaeli & Lahav 1991; Haehnelt & Tegmark 1996; Kashlinsky & Atrio-Barandela 2000; Aghanim, Górski & Puget 2001).

Theoretical predictions of the rms peculiar velocity of halos have been made by applying the peak theory (Bardeen et al. 1986), in which halo motions are modeled by motions of density peaks of a suitably smoothed version of the initial Gaussian density field evolved according to linear theory (Peebles 1980). $`N`$-body simulations have been used to test those predictions (e.g. Bahcall, Gramann, Cen 1994; Croft & Efstathiou 1994; Suhhonenko & Gramann 1999, Colberg et al. 2000; Sheth & Diaferio 2001; Yoshida, Sheth & Diaferio 2001). Colberg et al. (2000) found that, at $`z>4`$, the mean peculiar velocities of materials that end up in massive halos agree well with the predictions of linear theory for the corresponding peaks, but that linear theory systematically underestimates their subsequent growth; the rms peculiar velocities at $`z=0`$ evaluated from $`N`$-body simulations are larger than the linear theory prediction by up to 40 percent. This discrepancy may be ascribed to the following facts: that evolution of peculiar velocities of halos in high density regions are accelerated relative to the average (Tormen & Bertshinger 1996; Cole 1997; Kepner, Summers & Strauss 1997) and that halos tend to be formed in such higher-density regions because of the biased structure formation mechanism. Sheth & Diaferio (2001) have shown using $`N`$-body simulations that the evolution of the halo peculiar velocity indeed depends on the local matter density, rather than on the global density. They presented a simple model for the shape of the distribution function of the peculiar velocities that takes into account their dependence on the local density in an empirical manner.

The purpose of this paper is to present a simple and accurate model for the rms of peculiar velocities of halos. Hereafter, we denote the rms peculiar velocity of halos with mass $`M`$ as $`\sigma _{\mathrm{halo}}(M)`$.

The major improvements that we achieve here are to take into account the following two points that play an important role in understanding statistical properties of halo peculiar velocities. (i) The evolution of the halo peculiar velocity depends on the local matter density, rather than on the global density. We develop a model of the local density-dependent growth of the halo peculiar velocity in an empirical manner following Sheth & Diaferio (2001). (ii) Dark matter halos are biased tracers of the dark matter distribution, and thus halos tend to be located at high-density regions. In order to implement this, we use the halo biasing model developed by Mo & White (1986) and Sheth & Tormen (1999; 2002) combined with an empirical log-normal probability distribution function (PDF) of the mass density distribution (Kayo, Taruya & Suto; 2001).

The outline of this paper is as follows. In Section 2, we summarize linear theory predictions for the peculiar velocities of peaks. In Section 3, we describe $`N`$-body simulations that are used to calibrate and to test the model. We construct the model in Section 4. The model predictions are compared with $`N`$-body simulations in section 5. A summary and discussion are presented in Section 6.

## 2 Linear predictions for the peculiar velocities of peaks

In this section, we summarize basic equations of the linear theory predictions for peculiar velocities of density peaks.

According to linear theory of gravitational instability in a dust universe (Peebles 1980), the peculiar velocity field of mass grows as

$`v`$ $``$ $`a\dot{D}=\dot{a}D{\displaystyle \frac{d\mathrm{ln}D}{d\mathrm{ln}a}}`$ (1)
$`=`$ $`\dot{a}D(a)\left[{\displaystyle \frac{5\mathrm{\Omega }_0}{2D(a)X(a)^2}}+{\displaystyle \frac{\mathrm{\Lambda }_0a^2\mathrm{\Omega }_0/2a}{X(a)^2}}1\right],`$
where

$$X(a)=\frac{\dot{a}}{H_0}=\left[\frac{\mathrm{\Omega }_0}{a}+(1\mathrm{\Omega }_0\mathrm{\Lambda }_0)+\mathrm{\Lambda }_0a^2\right]^{\frac{1}{2}}.$$
(2)
The linear growth factor $`D(a)`$ is given by

$$D(a)=\frac{5\mathrm{\Omega }_0}{2}\frac{X(a)}{a}_0^a\frac{da^{}}{X(a^{})^3}.$$
(3)
Note that the frequently used approximation $`f(a)=d\mathrm{ln}D/d\mathrm{ln}a\mathrm{\Omega }^{0.6}`$ underestimates the growth rate by $`5`$ percent in the flat $`\mathrm{\Omega }_0=0.3`$ cosmology that we consider here. Thus, we will not adopt the approximation but we will use the numerical integration throughout the present paper.

We assume that the primordial (linear) density fluctuations obey the Gaussian statistics. Then the peculiar velocity field is isotropic and Gaussian, and its three-dimensional dispersion smoothed over a smoothing scale of $`R`$ is given by

$$\sigma _v(R,a)=\dot{a}f(a)\sigma _1(R,a).$$
(4)
Following Bardeen et al. (1986), $`\sigma _j`$ is defined for any integer $`j`$ as

$$\sigma _j^2(R)=\frac{1}{2\pi ^2}𝑑kk^{2+2j}P(k)W^2(kR),$$
(5)
where $`W(x)`$ is the Fourier transform of the smoothing window, and $`P(k)`$ is the linear matter power spectrum. Throughout this paper, we use the real space top-hat window function which corresponds to $`W(x)=(3/x^3)[\mathrm{sin}(x)x\mathrm{cos}(x)]`$. We adopt the fitting function of a cold dark matter (CDM) power spectrum given by Bardeen et al. (1986), who show that the rms peculiar velocity at the peaks of the smoothed density field differs systematically from $`\sigma _v`$ by

$$\sigma _p(R)=\sigma _v(R)\sqrt{1\sigma _0^4/\sigma _1^2\sigma _1^2}.$$
(6)
Note that this expression does not depend on the height of peaks.

The conventional assumptions in modeling the peculiar velocities of halos are as follows. (i) The velocities of halos are identical to those of peaks (equation 6). The smoothing scale is chosen such that the filtered mass is equal to that of the halo considered, i.e., $`M=4\pi \overline{\rho }R^3/3`$ in the case of the real space top-hat window function. (ii) The halo velocities evolve according to linear theory (equation 1). However, the prediction based on these assumptions seems to underestimate the rms peculiar velocities of halos (Colberg et al. 2000).

## 3 $`N`$-body simulation

In constructing an improved model of the peculiar velocity of halos, we first use $`N`$-body simulations of a $`\mathrm{\Lambda }`$CDM cosmology ($`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Lambda }_0=0.7`$, $`h=0.7`$, and $`\sigma _8=1.0`$). The details of the $`N`$-body simulations are described in Jing & Suto (1998); see also Kayo et al. (2001). Briefly, the simulations employ $`256^3`$ CDM particles in a cubic box of $`300h^1`$Mpc on a side. The simulations are performed using a P3M code with the gravitational softening length of $`ϵ120h^1\mathrm{kpc}`$. The initial matter power spectrum is computed using the CDM fitting function given by Bardeen et al. (1986). We have three realizations which differ only in the phases of the initial density fluctuations. The particle mass ($`m_{\mathrm{part}}=1.34\times 10^{11}h^1M_{}`$) of the simulations is sufficiently small to guarantee that there are practically no discreteness effects on dark matter clustering on scales down to the softening length in the redshift range of interest for our purposes (Hamana, Yoshida & Suto 2002).

We identify dark matter halos using the standard friends-of-friends algorithm with a linking parameter of $`b=0.164`$ (in units of the mean particle separation). We set the minimum mass of the halos as $`1.34\times 10^{12}h^1M_{}`$, which corresponds to the mass of 10 simulation particles. We define the peculiar velocity of each halo to be the center-of-mass peculiar velocity of member particles.

The smoothed mass density fields are computed on $`64^3`$ regular grids in the simulation box using the top-hat window function. For each halo, the smoothed mass density of the nearest grid is assigned as its local background density.

Fig. 1 shows the rms peculiar velocities of halos in the $`N`$-body simulation (filled circles) together with the linear theory prediction of peaks, (equation 6; dashed curves) as well as our predictions described in the next section (solid lines and long-dashed lines). Note that, in order to account for the finite size of the simulation box ($`L_{\mathrm{box}}=300h^1`$Mpc), we set a lower limit of the integration in equation (5) by $`k_{min}=2\pi /L_{\mathrm{box}}`$. This large-scale cut-off reduces $`\sigma _p`$ by a fractional amount of about 8 percent, and is applied to all model predictions in this paper. The plot confirms the previous finding that the linear theory predictions underestimate the rms peculiar velocities by about 10–40 percent for cases $`0<z<2.2`$ in a $`\mathrm{\Lambda }`$CDM cosmology. In addition, it is seen from this Figure that the linear theory predictions have a slightly steeper slope than that of $`N`$-body simulations, especially for the case of $`z=0`$. This discrepancy may be explained by the following: (i) massive halos tend to be formed in denser regions as a result of the biasing mechanism, and (ii) halos in denser regions move faster than those in less dense regions, as will be proved in the next section (Fig. 2 and 4). As a consequence, on average, deviations from linear theory appear larger for massive halos, and thus the slope of $`N`$-body data becomes flatter.

## 4 Model

In this section, we describe our model of the rms peculiar velocity of halos. We basically follow the framework developed by Sheth & Diaferio (2001); see also Diaferio & Geller (1996) and Sheth (1996). They originally developed the model for the distribution function of peculiar velocities of dark matter particles and halos, but their framework can be straightforwardly applied to the model for the rms of halo peculiar velocities.

Sheth & Diaferio (2001) pointed out that the evolution of the halo peculiar velocity depends on the local matter density (defined with an appropriate smoothing length which we discuss later). We thus first construct a model for the dependence of the halo peculiar velocity dispersion, $`\sigma _{\mathrm{halo}}^2(M,\delta )`$, on the local background density $`\delta `$. Following Sheth & Diaferio (2001), we consider a parametric model

$$\sigma _{\mathrm{halo}}^2(M,\delta )=[1+\delta (R_{\mathrm{local}})]^{2\mu (R_{\mathrm{local}})}\sigma _p^2(M),$$
(7)
where $`R_{\mathrm{local}}`$ is the smoothing scale with which the local background density $`\delta `$ is defined. The key question is how to define the appropriate smoothing length-scale. Evidently, $`R_{\mathrm{local}}`$ should enclose the gravitational coherence scale which is responsible for the local deviation of the peculiar motion of a halo from its global value (for instance, given by linear theory). Thus, the size of halo (e.g., the virial radius) sets a lower limit for $`R_{\mathrm{local}}`$. On the other hand, choosing values for $`R_{\mathrm{local}}`$ that are too large would result in an excessively smoothed density field, and local density fluctuations, by which the deviation from the global value is induced, are smoothed out. We may, thus, expect that $`R_{\mathrm{local}}`$ is in the range between the linear scale (such that $`\sigma _0(R)<1`$) and the non-linear scale ($`\sigma _0(R)>1`$). Concerning the power-law index $`\mu `$, the linear theory relation, $`f\mathrm{\Omega }^{0.6}`$, suggests $`\mu 0.6`$ (Sheth & Diaferio 2001). Note, however, that $`\mu `$ depends on $`R_{\mathrm{local}}`$ on which the local density is defined (a larger $`\mu `$ for a larger smoothing length, because the range of $`\delta `$ becomes narrower), thus $`\mu `$ does not need to be very close $`\mu 0.6`$, but it is expected that with an appropriate choice of the smoothing length-scale, $`\mu `$ would be close to 0.6.

To proceed further, we adopt an ansatz that $`R_{\mathrm{local}}`$ is given by a relation $`\sigma _0(R_{\mathrm{local}})=\sigma _{\mathrm{local}}`$, and we determine the model parameter $`\sigma _{\mathrm{local}}`$ empirically using $`N`$-body simulations. As expected from the discussion above, a suitable choice of $`\sigma _{\mathrm{local}}`$ would be within a range $`0.1<\sigma _{\mathrm{local}}<1`$. We find that choosing $`\mu =0.6`$ (constant) with $`R_{\mathrm{local}}`$ given by $`\sigma _{\mathrm{local}}=0.5(1+z)^{0.5}`$ (through $`\sigma _0(R_{\mathrm{local}})=\sigma _{\mathrm{local}}`$) provides reasonable fits to the results of the $`N`$-body simulations as shown in Fig. 2. For this choice, $`R_{\mathrm{local}}18`$ and $`16h^1`$Mpc for z=0 and 0.8, respectively. The value of $`\mu =0.6`$ is consistent with one suggested in linear theory ($`\mu 0.6`$). Thus, while the choice of the smoothing length may still be arbitrary, we decide to use the above sets of parameters for definiteness.

Fig. 3 shows the peculiar velocity dispersion of halo progenitors as a function of the local density at the initial redshift of the $`N`$-body simulation $`z=18`$. We mark the member particles of a halo identified at $`z=0`$, and then compute the position and velocity of the halo progenitor from the center-of-mass position and velocity of the marked particles at $`z=18`$. Fig. 3 clearly indicates that there is no clear correlation between the peculiar velocity and the local density at the initial epoch. Together with Fig. 2, this strongly supports the idea that most of the halo peculiar velocity is acquired through the subsequent evolution of the nearby density field. A closer look at Fig. 3 reveals that the initial peculiar velocity dispersions decline very slightly at both the lowest and highest bins. An interpretation of this is that, at the initial epoch, the halo progenitors locate close to positive (negative) density peaks tend to have relatively higher (lower) local densities. At the same time, a potential gradient at positions close to the density peaks tends to be small, and thus the progenitors at such regions move slower than others. As a result, the highest and lowest density bins have a slightly lower peculiar velocity than other bins.

As shown by Sheth & Diaferio (2001), if the range of halo masses and local background density is sufficiently small, the PDF of halo peculiar velocities is well approximated by a Maxwellian distribution, i.e., the PDF of one-dimensional velocity is a Gaussian. We denote by $`p(M|\delta )`$ the probability of finding a halo with mass $`M`$ at a region where the background density is $`\delta `$. Then the peculiar velocity dispersion of halos with mass $`M`$ is given by summing up the dispersion $`\sigma _{\mathrm{halo}}^2(M,\delta )`$ weighted over the probability of finding halos in regions of $`\delta `$:

$$\sigma _{\mathrm{halo}}^2(M)=\frac{𝑑\delta p(M|\delta )\sigma _{\mathrm{halo}}^2(M,\delta )}{𝑑\delta p(M|\delta )}.$$
(8)
Note that this expression is valid only if the background density is defined appropriately (i.e., with the smoothing scale $`R_{\mathrm{local}}`$ defined above). The expression for $`p(M|\delta )`$ is obtained as follows.

It is well established that dark matter halos are biased tracers of the underlying dark matter distribution. In this picture, it is useful to write the conditional probability as

$$p(M|\delta )=\frac{n(M|\delta )}{n(M)}p_{dm}(\delta ),$$
(9)
where $`n(M)`$ and $`n(M|\delta )`$ denote the halo mass function and the conditional mass function in regions with the local background density $`\delta `$, respectively, and $`p_{dm}(\delta )`$ is the PDF of the dark matter density field. To evaluate the first term of the right hand side of equation (9), i.e. the excess number density of halos, we adopt two approaches. The first is to directly use the mass functions derived from the extended model of Press-Schechter (1974); see, for example, Bond et al. (1991). Lacey & Cole (1991) and Mo & White (1996). In practice, we adopt the conditional mass function given by Sheth & Tormen (2002) that even takes into account the effect of the ellipsoidal collapse. The other approach is to use a simple approximation, $`n(M|\delta )[1+b(M)\delta ]n(M)`$. Note that this approximation is valid only if the smoothing scale defining the local density is much larger than the halo size, or equivalently if the total mass enclosed within the smoothing scale is much larger than the halo mass. Adopting this, equation (9) reduces to

$$p(M|\delta )=[1+b(M)\delta ]p_{dm}(\delta ).$$
(10)
Mo & White (1996) developed a Lagrangian halo bias model based on the extended Press-Schechter model. Jing (1998) and Sheth & Tormen (1999) discussed the correction for the mass dependence of the halo bias and proposed modified fitting functions that give better fits to numerical simulation results. In what follows, we use the linear bias model of Sheth & Tormen (1999). Although the former approach may lead to a better model than the latter, we examine the latter approximation to see if the simple approach could provide an acceptable model. Fig. 4 compares the above two theoretical predictions (dashed and dotted lines) for $`n(M|\delta )/n(M)`$ with the $`N`$-body data (solid histograms) for ranges of the local background density. The simple approximation $`n(M|\delta )[1+b(M)\delta ]n(M)`$ turns out to provide a good fit to the $`N`$-body data for smaller mass ranges $`M_{\mathrm{halo}}<10^{1314}h^1M_{}`$, but underestimates the conditional probability for larger mass ranges. The conditional mass function by Sheth & Tormen (2002) gives good fits for high local densities ($`\delta >1`$, upper two cases in the plots), while fits become degraded for low densities ($`\delta <1`$, lower two cases in the plots). Importantly, it reproduces the cut-off in the conditional mass function at a large mass nicely in marked contrast to the approximation which does not have the cut-off.

The dark matter PDF is known to be well fitted to the log-normal distribution (e.g., Coles & Jones 1991; Kofman et al. 1994; Kayo, Taruya & Suto 2001; see also Taruya, Hamana & Kayo 2003 for a possible explanation of the origin of the log-normal PDF)

$$p_{ln}(\delta )=\frac{1}{(2\pi \sigma _1^2)^{1/2}}\mathrm{exp}\left\{\frac{[\mathrm{ln}(1+\delta )+\sigma _1^2/2]^2}{2\sigma _1^2}\right\}\frac{1}{1+\delta }$$
(11)
where $`\sigma _1`$ is related to the variance of the nonlinear density field by

$$\sigma _1^2(R)=\mathrm{ln}[1+\sigma _{nl}^2(R)],$$
(12)
with

$$\sigma _{nl}^2(R)=\frac{1}{2\pi ^2}𝑑kk^2P_{nl}(k)W^2(kR).$$
(13)
In the above expression, $`P_{nl}(k)`$ denotes the nonlinear matter power spectrum for which we compute using the fitting function of Peacock & Dodds (1996).

## 5 Results

Fig. 1 compares our model predictions with $`N`$-body simulations (filled circles); long-dashed curves correspond to equation (9) with the conditional mass function by Sheth & Tormen (2002), while solid curves correspond to the approximation equation (10). It is clear that our models (especially the former one) substantially improve the accuracy of the prediction compared with that of linear theory (dotted lines) except for a very high mass range of $`M>10^{14}h^1M_{}`$. However, the predictions of the model adopting the approximation equation (10) rise very rapidly at a high mass range of $`M>10^{14}h^1M_{}`$ and significantly deviate from the simulation results. This is explained as follows: As shown in Fig. 4, the conditional mass function computed from $`N`$-body simulation declines at a high mass range. However, the approximation equation (10) does not reproduce this feature but keeps on rising toward a large halo mass. Consequently, the approximated model over-estimates the contribution from moderately over density regions. Here, it must be noted that with our choice of the smoothing scale ($`R_{\mathrm{local}}=18`$, 15 and 13$`h^1`$Mpc for $`z=1`$, 2 and 3, respectively), the most of the matter in the universe is in regions where the local density contrast being $`1<\delta <1`$. Therefore, a discrepancy in the conditional mass function at a high density range of $`\delta >1`$ between the model and the simulation does not have a large influence on the model prediction. Adopting the conditional mass function proposed by Sheth & Tormen (2002), which reasonably reproduces the decline as shown in Fig. 4, improves the accuracy of the prediction very well. In particular, this model predicts a flatter slope than that of the linear theory prediction, and it is very close to the simulation results. We, however, note that the model prediction adopting Sheth & Tormen conditional mass function rises at very large mass scales ($`M>10^{15}`$) and deviates from the simulation. The reason of this is that the conditional mass function becomes less accurate at such mass range. Because at that mass range, the number density of halos is too small to have a meaningful statistic, this disagreement is not serious.

Although there is still room for improvement in the model, our model adopting the conditional mass function of Sheth & Tormen (2002), provides a good fit to the $`N`$-body data over a wide range of the halo mass. To improve the model prediction, it may be necessary to adopt a further developed conditional mass function. Developing an accurate conditional mass function is beyond the scope of the present paper, and hence we do not attempt this here. In the following discussions, we adopt the model adopting Sheth & Tormen conditional mass function.

As Sheth & Diaferio (2001) originally showed, the model presented in this paper is applied to model the PDF of halo peculiar velocities

$$P(v|M)=\frac{𝑑\delta p(M|\delta )P(v|M,\delta )}{𝑑\delta p(M|\delta )},$$
(14)
where $`P(v|M)`$ is the PDF of peculiar velocities of halos with mass $`M`$, and $`P(v|M,\delta )`$ is the PDF for halos at regions with the local background density $`\delta `$. The PDF of one-dimensional peculiar velocities (which we denote by $`v_1`$), $`P(v_1|M,\delta )`$, is well approximated by the Gaussian distribution (Sheth & Diaferio 2001). We compute $`\sigma `$ of the Gaussian PDF using equation (7) with a reasonable isotropic assumption $`v_1^2=v^2/3`$. The results are plotted in Fig. 5 (solid lines), and are compared with $`N`$-body simulations (histograms) and Gaussian PDF (dotted lines). As Fig. 5 clearly shows, our model agrees with the $`N`$-body simulation very well. Also it is seen from the Figure that the PDF of one-dimensional peculiar velocity consists of two parts; a Gaussian core and an exponential tail. The exponential tail arises from the sum of Gaussian distributions with different dispersions (Sheth & Diaferio 2001; Kuwabara, Taruya & Suto 2002). Our model nicely reproduces the deviation from the Gaussian distribution at large peculiar velocities and agrees with the simulation very well.

The model presented in this paper can be also applied to the pairwise velocity dispersion of halos at large distances. The model of the pairwise velocity dispersion plays an important role in computing the redshift distortion of the two-point correlation function of halos. Because any velocity correlation between halos with large separations of $`r>20h^1`$Mpc may be safely ignored, the halo pairwise velocity dispersion is given by $`v_{12}^2=2/3v^2`$, where $`v_{12}`$ denotes the pairwise velocity of a pair of halos. Often $`\sigma _p^2`$ (equation 6) is used to compute $`v^2`$ (e.g., Hamana et al. 2001); see Sheth et al. 2001a,b for an alternative approach. However, as shown in the previous sections, the linear theory prediction underestimates the halo velocity and our model provides better fits to $`N`$-body data. We may, therefore, expect that our model improves predictions of the pairwise velocity dispersion. Fig. 6 compares the model prediction for the pairwise velocity dispersion (solid lines) with the result of our $`N`$-body simulation (filled circles). We also plot the linear theory prediction for peaks (dotted lines). As Fig. 6 clearly shows, our model indeed improves the accuracy of prediction at $`r>20h^1`$Mpc.

## 6 Summary and discussions

We present a model of peculiar velocities of dark matter halos which describes various statistical properties well and significantly improves the accuracy of the simple linear theory predictions, including the halo rms velocity dispersion, PDF, and the pairwise velocity dispersion at large separations. The model predictions, however, becomes less accurate for very large halo masses ($`M>10^{15}h^1M_{}`$). This is mainly due to the limitation of the adopted conditional halo mass function, and may be reconciled by using an improved model. Note that this shortcoming hardly causes any practical problem, because the number density of such very massive halos is too small to have a meaningful statistic.

The model presented in this paper is characterized by a single free parameter, $`\mu (R_{\mathrm{local}})`$, that we determine in an empirical manner. To do this, we need the appropriate smoothing scale for the current purpose. We make an anzats that $`R_{\mathrm{local}}`$ is in the range between the linear and the non-linear scale and it should be given by the relation $`\sigma _0(R_{\mathrm{local}})=\sigma _{\mathrm{local}}`$ for some value in a range $`0.1<\sigma _{\mathrm{local}}<1`$. This turns out to be reasonably successful, and we have found that a choice of $`\mu =0.6`$ with $`\sigma _{\mathrm{local}}=0.5(1+z)^{0.5}`$ works well for all redshifts we consider in this paper ($`z<2.2`$). This choice gives the smoothing scales of $`R_{\mathrm{local}}=18`$, 15, 11, 8.5 and 5.6$`h^1`$Mpc for $`z=0`$, 1, 3, 5 and 10, respectively. However, we calibrated the model with $`N`$-body data at $`0<z<2.2`$ only, thus the validity of this procedure is questionable for very higher redshifts (although in practice velocities of halos at very high redshifts are not observable currently). A more sophisticated approach to the choice of the smoothing scale is needed for a further improvement of the model. This may be closely connected to an understanding of the role of the local environment in the structure formation, and is a subject of work in progress.

## Acknowledgments

We thank the referee, Ravi Sheth, for detailed and constructive comments on the earlier manuscript which significantly improves various aspects of the paper. We thank Joerg M. Colberg for alerting us to the mistake in the original version of this paper. T.H. acknowledges supports from Japan Society for Promotion of Science (JSPS) Research Fellowships. Y.I.P. is supported in part by NKBRSF (G19990754) and by NSFC (No. 10125314). This work was supported in part by the Grant-in-Aid for Scientific Research of JSPS (12640231). Numerical computations presented in this paper were carried out at the Astronomical Data Analysis Center of the National Astronomical Observatory, Japan (project ID: mys02a).


text:
# Elliptic Flow of Identified Hadrons in Au+Au Collisions at √𝑠__NN = 200 GeV

## Abstract

The anisotropy parameter ($`v_2`$), the second harmonic of the azimuthal particles distribution, has been measured with the PHENIX detector in Au+Au collisions at $`\sqrt{s_{_{NN}}}=`$ 200 GeV for identified and inclusive charged particle production at central rapidities ($`|\eta |<0.35`$) with respect to the reaction plane defined at high rapidities ($`|\eta |=`$3–4). We observe that the $`v_2`$ of mesons falls below that of (anti)baryons for $`p_\mathrm{T}>`$ 2 GeV/c, in marked contrast to the predictions of a hydrodynamical model. A quark coalescence model is also investigated.

Event anisotropy is expected to be sensitive to the early stage of ultra-relativistic nuclear collisions at the Relativistic Heavy Ion Collider (RHIC). The possible formation of a quark-gluon plasma (QGP) could affect how the initial anisotropy in coordinate space is transferred into momentum space in the final state. The anisotropy parameter $`v_2`$ for a selection of produced particles is derived from the azimuthal distribution of those particles.

$$\frac{dN}{d\varphi }\mathrm{\hspace{0.33em}1}+\mathrm{\hspace{0.33em}2}v_2\mathrm{cos}2(\varphi \mathrm{\Phi }_{\mathrm{RP}})$$
(1)
where $`\varphi `$ is the azimuthal direction of the particle and $`\mathrm{\Phi }_{\mathrm{RP}}`$ is the direction of the nuclear impact parameter (“reaction plane”) in a given collision. Measurements of the parameter $`v_2`$ in RHIC collisions have been performed  for charged particles and for identified particles. The current work reports results for charged particles versus transverse momentum ($`p_\mathrm{T}`$) out to 5 GeV/c, and extends previous measurements for identified particles out to 3 GeV/c for $`\pi `$ and K, and to 4 GeV/c for protons. (Previous measurements of the $`v_2`$ for $`\pi `$, K, and p extended to 1 GeV/c at $`\sqrt{s_{_{NN}}}=`$ 130 GeV.) Detailed measurements of the azimuthal anisotropy are important to eventually discriminate among different possible scenarios for its physical origin. Such scenarios include: hydrodynamical flow of compressed hadronic matter, the production of multiple mini-jets, and an anisotropy developed during an early quark-matter phase of the collision. It has been observed that $`v_2`$ saturates at $`p_\mathrm{T}`$ $``$ 2 GeV/c and above . The cause of this saturation is not yet known; however, we note that at this momentum the particle composition is very different than at low momentum in that the proton yield is comparable to the pion yield . This makes the measurement of $`v_2`$ for separately identified particles especially interesting.

The measurements described here were carried out in the PHENIX experiment at RHIC . About 28 M minimum bias Au+Au collisions at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 200 GeV from the 2001-2002 run period (Run-2) are used in the analysis. Charged particles are measured in the central arm spectrometers ($`|\eta |<0.35`$)  where PHENIX has excellent particle identification capabilities. The drift chamber (DC) and the first pad chamber plane (PC1) together with the collision vertex define the charged particle tracks. In order to reduce background, the reconstructed tracks are confirmed by requiring matching hits in the outer detectors, i.e. the third pad chamber plane (PC3) and the electromagnetic calorimeter (EMCAL) or the time-of-flight detector (TOF). In this analysis, the TOF detector is used to identify charged particles up to 4 GeV/c in $`p_\mathrm{T}`$. Particle time-of-flight is measured using the TOF with respect to the collision time defined by beam counters (BBC), and is used to calculate mass squared using the particle momentum and the flight path length . The timing resolution of the system is $``$ 120 ps. A momentum dependent $`\pm 2`$ $`\sigma `$ cut on mass squared allows particle identification in the $`p_\mathrm{T}`$ range 0.2 $`<`$ $`p_\mathrm{T}`$ $`<`$ 3 GeV/c for pions, 0.3 $`<`$ $`p_\mathrm{T}`$ $`<`$ 3 GeV/c for kaons, and 0.5 $`<`$ $`p_\mathrm{T}`$ $`<`$ 4 GeV/c for protons. The contamination of mis-identified particles is less than 10$`\%`$. In addition to collision time, the BBC provide z-vertex position information. The two beam counters are located at $`|z|`$=1.5 m from the collision point, covering $`|\eta |=34`$. They consist of 64 photo-multiplier tubes (PMT) equipped with quartz Cherenkov radiators in front surrounding the beam pipe. The large charged multiplicity (a few hundred) in $`|\eta |=34`$ and the non-zero signal of event anisotropy in this $`\eta `$ range enables us to estimate the azimuthal angle of the reaction plane in each event using the BBC with full azimuthal angle coverage.

Since the $`v_2`$ parameter is in effect a quadrupole moment, the anisotropy which gives rise to a non-zero $`v_2`$ is often referred to as an “elliptic flow.” It is extracted by first determining the reaction plane angle $`\mathrm{\Phi }_{\mathrm{RP}}`$ for each event,

$$\mathrm{tan}2\mathrm{\Phi }_{\mathrm{RP}}=\frac{\mathrm{\Sigma }n_{\mathrm{ch}}\mathrm{sin}2\varphi _{\mathrm{PMT}}}{\mathrm{\Sigma }n_{\mathrm{ch}}\mathrm{cos}2\varphi _{\mathrm{PMT}}}$$
(2)
where $`n_{\mathrm{ch}}`$ is the number of charged particles per PMT (determined from the pulse height in each PMT) and $`\varphi _{\mathrm{PMT}}`$ is the azimuthal angle of each PMT. Then, it is calculated by the Fourier moment $`v_2=\mathrm{cos}2(\varphi \mathrm{\Phi }_{\mathrm{RP}})`$ over all particles, for all events in a given sample . Corrections  are applied to account for finite resolution in the reaction plane determination, and for possible azimuthal asymmetries in the reaction plane detector response. The bottom-left panel in Fig. 1 shows the average cosine of the difference between the two reaction planes defined by the beam counters at $`\eta `$ = $`34`$ and at $`\eta `$ = $`43`$ using the the elliptic (second) moment definition. In order to improve the reaction plane resolution, a combined reaction plane is defined by averaging the reaction plane angles obtained from each BBC, using the elliptic moment in each. The estimated resolution of the combined reaction plane, $`\mathrm{cos}2(\mathrm{\Phi }_{\mathrm{measured}}\mathrm{\Phi }_{\mathrm{true}})`$, has an average of 0.3 over centrality with a maximum of about 0.4. The estimated correction factor, which is the inverse of the resolution for the combined reaction plane, is shown in the top-left panel in Fig. 1.

The present technique is distinguished by defining the reaction plane angle using particles at high rapidity when measuring $`v_2`$ for particles at mid-rapidity. Other measurements of $`v_2`$ for mid-rapidity particles at RHIC have used reaction planes defined with mid-rapidity particles; or have employed a technique of measuring angular correlations between pairs of particles at mid-rapidity. While these different approaches generally seek to measure the same thing, they are not identical and a variety of physics effects can cause them to yield different results from the same collision sample . Because of the large rapidity gap between the reaction plane and the mid-rapidity acceptance of about 3 units, it is expected that this analysis is less affected by non-flow contributions. However, we do not observe any substantial difference between the $`v_2`$ results shown here and published results for the $`v_2`$ of charged particles at RHIC in the $`p_\mathrm{T}`$ range where they are available.

The centrality of each collision is defined using the simultaneous measurement of the total number of particles measured in the BBC and the total energy measured in the zero degree calorimeter . The middle panel in Fig. 1 shows the centrality dependence of $`v_2`$ for charged particles measured at mid-rapidity ($`|\eta |<0.35`$) with respect to the reaction plane defined above. The centrality is measured in percentile from the most central collision. The $`v_2`$ parameter decreases for both peripheral and central collisions with a maximum at about 50$`\%`$ of the geometric cross section. Beyond 70$`\%`$, the correction factor due to the reaction plane resolution is large, as shown the left-most panel in the Fig. 1. This limits the centrality range used in this analysis.

The right-most panel in Fig. 1 shows the transverse momentum dependence of $`v_2`$ for charged particles with respect to the reaction plane for minimum-bias events. The data above a $`p_\mathrm{T}`$ of 2 GeV/c clearly show a deviation from the monotonically increasing behavior seen at smaller $`p_\mathrm{T}`$. The systematic errors are shown as line bands, which are estimated by several reaction plane methods using the two single beam counters or combined beam counters and by several different ways to correct non-uniform reaction plane distribution: “inverse weighting,” “re-centering of sine and cosine summation,” “Fourier expansion” and combinations of those above . Those systematic errors are estimated to be about 10$`\%`$, depending on centrality, and are independent of $`p_\mathrm{T}`$. Above 3 GeV/c, background tracks result in an additional systematic error of about 10$`\%`$, depending on $`p_\mathrm{T}`$, which is included in the upper error band .

In Fig. 2, the transverse momentum dependence of $`v_2`$ for identified particles is shown. The top-left panel shows negatively charged particles, while the top-right panel shows positively charged particles as described in the figure caption. The statistical errors and the systematic errors are plotted independently. From the lambda particle spectra measured in the PHENIX central arm, it is determined that approximately 35$`\%`$ of the protons originate from lambda decays (“lambda feed-down”). The effect of the lambda feed-down on the measured $`v_2`$ of the proton is studied by varying the lambda $`v_2`$ with Monte Carlo simulation. Protons resulting from lambda feed-down increase the measured $`v_2`$ value. Using the value of the lambda $`v_2`$ measured at $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 130 GeV at RHIC , the effect on the proton $`v_2`$ would be less than 10$`\%`$. Less than 5$`\%`$ of protons originate from decays of particles not involving the lambda. Based on further simulations of their decays to protons, we estimate that the total systematic error due to feed-down is at most 11$`\%`$ depending on $`p_\mathrm{T}`$, which is included in the lower sytematic error band in Fig. 2.

The combined positive and negative particles are shown in the bottom-left panel. The lines in that panel represent a hydrodynamical calculation including a first-order phase transition with a freeze-out temperature of 120 MeV. The data show that at lower $`p_\mathrm{T}`$ ($`<2`$ GeV/c), the lighter mass particles have a larger $`v_2`$ at a given $`p_\mathrm{T}`$, which is reproduced by the model calculations. We note, however, that the difference between the charged kaons and charged pions is larger than the model predicts.

A striking feature observed at higher $`p_\mathrm{T}`$ is that the $`v_2`$ of $`p`$ and $`\overline{p}`$ are larger than for $`\pi `$ and $`K`$ at $`p_\mathrm{T}`$ $`>`$ 2 GeV/c. This is in sharp contrast to the hydrodynamical picture, which would predict the same mass-ordering for $`v_2`$ at all $`p_\mathrm{T}`$. In our data the mesons begin to show a departure from the hydrodynamical prediction at $`p_\mathrm{T}`$ of about 1.5 GeV/c, while the (anti)baryons agree with the prediction up until 3 GeV/c but may be deviating at higher $`p_\mathrm{T}`$. Such behavior is predicted by the quark coalescence mechanism , as shown in the bottom-right panel where both $`v_2`$ and $`p_\mathrm{T}`$ have been scaled by the number of quarks. This could be an indication that the $`v_2`$ of measured hadrons is already established in a quark-matter phase, although it does not explain why the quark $`v_2`$ would saturate with $`p_\mathrm{T}`$. There exist other scenarios that could be applicable at RHIC, but we have selected two simple models (hydrodynamical and quark coalesence) only to emphasize the experimental evidence of the crossing of $`v_2`$ for mesons and baryons.

As an additional illustration of the different behavior for mesons and baryons, the transverse momentum dependences of the $`v_2`$ parameter are shown in Fig. 3 for different particles and different centralities. Since the particle identification separation of $`K`$ and $`p`$ goes up to 4 GeV/c, the combined $`\pi `$ and $`K`$ can be compared with protons up to 4 GeV/c. The charged particle acceptance is larger than the TOF acceptance where the particle identification can be performed. Therefore, the statistical fluctuations for the charged particle $`v_2`$ are smaller than for the $`p`$, $`\overline{p}`$ and $`\pi +K`$. The trend exhibited in Fig. 2 for minimum bias spectra, in which the $`v_2`$ parameter for (anti)baryons exceed those for mesons at $`p_\mathrm{T}`$ $`>`$ 2 GeV/c, is shown here to occur for all centralities.

In summary, the value of the $`v_2`$ parameter for identified and inclusive charged particle production at mid-rapidity has been measured with respect to the reaction plane defined in the forward and backward rapidity regions in $`\sqrt{s_{_{\mathrm{NN}}}}`$ = 200 GeV Au+Au collisions, using the PHENIX experiment at RHIC. The value of $`v_2`$ for charged particles decreases for both peripheral and central collisions with a maximum at about the 50th percentile of the geometric cross section. We have observed that for charged particles $`v_2`$ increases with $`p_\mathrm{T}`$ up to about 2 GeV/c, then starts to saturate or decrease slightly. However, the detailed behavior is different for different particle species. The lighter particles have larger $`v_2`$ than the heavier particles for $`p_\mathrm{T}`$ below 2 GeV/c. This trend is partly reversed above 2 GeV/c where the proton and anti-proton have larger $`v_2`$ than mesons, a pattern which persists over all centralities. A hydrodynamical calculation can reproduce the mass-ordering and magnitude of $`v_2`$ for the different particles in the region up to 2 GeV/c, but fails to reproduce either in the $`p_\mathrm{T}`$ region above 2 GeV/c. As an alternative, we investigated the quark-coalescence scenario, in which the anisotropy of the final-state hadrons is largely inherited from the anisotropy of quarks in a preceding quark-matter phase. The quark-coalescence model makes a definite prediction for a simple scaling behavior between the $`v_2`$ for mesons and for (anti)baryons, and this scaling behavior is largely, though not perfectly, borne out in our data. Further measurements extending to higher $`p_\mathrm{T}`$ involving more identified species will be required to discriminate among alternative scenarios for the origin of elliptic flow at RHIC.

We thank the staff of the Collider-Accelerator and Physics Departments at BNL for their vital contributions. We acknowledge support from the Department of Energy and NSF (U.S.A.), MEXT and JSPS (Japan), CNPq and FAPESP (Brazil), NSFC (China), CNRS-IN2P3 and CEA (France), BMBF, DAAD, and AvH (Germany), OTKA (Hungary), DAE and DST (India), ISF (Israel), KRF and CHEP (Korea), RAS, RMAE, and RMS (Russia), VR and KAW (Sweden), U.S. CRDF for the FSU, US-Hungarian NSF-OTKA-MTA, and US-Israel BSF.


text:
# Scaling and universality in the micro-structure of urban space

## 1 Introduction

The morphology of urban settlements and its dynamics has captured the interest of physicists  as it may shed light on Zipf’s law for cities , challenge theoretical frameworks for cluster dynamics or improve predictions of urban growth .

The search for a unified theory of urban morphology has focused on the premise that cities can be conceptualized at several scales as fractals. At the regional scale, rank-order plots of city size follow a fractal distribution  and population scales with city area as a power-law . More recently, it has been observed that the area distribution of satellite cities, towns and villages around large urban centres also obeys a power-law with exponent $`2`$ . At the scale of transportation networks, railway networks appear to have a fractal structure . At the scale of the neighbourhood, it has been suggested that urban space resembles a Sierpinsky gasket . These scales are inter-related as summed up in \[1, pp. 241\] (author’s translation from French): ’*Polycentric growth, which is connected to the non-homogeneous distribution of pre-urban cores and the birth of a hierarchy of sub-centres, influences the morphology of the transport network, which plays in itself an important role in axial growth and therefore for the future spatial development of the urbanised area*’.

The fractal dimensions of US cities and international cities have values ranging from $`1.2778`$ (Omaha, ) to $`1.93`$ (Beijing, ), where the fractal dimension of large cities tends to cluster around the latter value . Studies of urban growth of London between $`1820`$ and $`1962`$ show that fractal dimensions for this period vary from $`1.322`$ to $`1.791`$ . The fractal dimensions for the growth of Berlin in $`1875`$, $`1920`$ and $`1945`$ are $`1.43`$, $`1.54`$ and $`1.69`$, respectively . The fractal dimension of urban aggregates is a global measure of areal coverage, but detailed measures of spatial distribution are clearly needed to complement adequately the description of the morphology of an urban area . Further, current approaches to data collection and modelling identify cities as fractal only on the urban periphery of the giant urban cluster that grows around the city core (or central business district), as clusters become compact at distances close to the centre of the city . Although remote sensing techniques are promising in extracting urban morphology with greater detail, available studies have, to our knowledge, been limited to individually selected, medium scale cities (see e.g. ).

Hillier and Hanson  suggest an underlying structure to urban open space that is determined by the complexity of buildings which bound the space . Urban space available for pedestrian movement, excluding by definition physical obstacles, is relatively linear. To extract this linearity, a street can be approximated by a set of straight lines. The global set of lines, the so-called *axial map*, is defined as the least number of longest straight lines. An axial map can be derived by drawing the longest possible straight line on a city map, then the next longest line, so-called axial line, until the open space is crossed and ’all axial lines that can be linked to other axial lines without repetition are linked’ . Figure 3 shows several axial maps.

Axial maps are relevant to urban theorists, as they encode the structure of open space in urban settlements and provide a simplified signature of the growth process. One could hope to extract the spatio-temporal dynamics of axial map growth by analysing a sequence of aerial photographs of the urban periphery over a period of growth. Conversely, one could hope to model urban growth as the trajectory of $`N`$ walkers on a plane, the step of the walk being axial line length.

Here we show that we can rescale axial line length and rank to obtain two distinct rank-order curves that provide a classification for several cities independently of city size. We also show that there is a class of cities that does not obey this classification. The collapse of curves suggests that spatial fluctuations in the length of urban linear structures, differing in size and location, are governed by similar statistical rules and supports the hypothesis that the linear dimension of large scale structures in cities reflects generic properties of city growth .

## 2 Structure of urban space

### 2.1 Intermittency in urban space

Let $`𝐥_i=\left\{𝐥_{i,j}\right\}`$, $`j=1,\mathrm{},N_i`$, be the $`N_i`$ axial lines associated with city $`i`$. Each axial line, $`𝐥_{i,j}`$ ($`j=1,\mathrm{},N_i`$) is defined by the coordinates of its extremities

$$𝐥_{i,j}=\{(x_{(i,j),1},y_{(i,j),1}),(x_{(i,j),2},y_{(i,j),2})\}$$

The axial map of city $`i`$, $`𝒞_i`$, is thus a set of $`N_i`$ points on a fourth dimensional space, $`𝒞_i=\left\{(\rho _{i,j},\theta _{i,j},l_{i,j},\phi _{i,j})\right\}`$, where $`(\rho _i,\theta _i)`$ are the polar coordinates of the axial line geometric centre, $`𝐬_{i,j}=(\frac{x_{(i,j),1}+x_{(i,j),2}}{2},\frac{y_{(i,j),1}+y_{(i,j),2}}{2})`$, and $`(\pm \frac{l_{i,j}}{2},\phi _{i,j})`$ are the polar coordinates of the axial line’s extremities on its geometric centre reference system, $`𝐝_{i,j}=\pm (\frac{\left|x_{(i,j),1}x_{(i,j),2}\right|}{2},\frac{\left|y_{(i,j),1}y_{(i,j),2}\right|}{2})`$.

Coordinates $`\rho `$ and $`\theta `$ encode the geographic location of axial lines. The unconditional distribution of $`\phi `$ is multimodal for rather general families of urban settlements. This occurs, for example, when land is partitioned in clusters of randomly oriented orthogonal grids. Nevertheless, the unconditional distribution of $`l`$ is unimodal and skewed to the right (see Figure 1), and, thus, the only coordinate which is a good candidate for inspection of intermittency in urban space. In what follows we will analyse the statistics of line length for a database of cities. We start by fitting the data for Tokyo to a stretched exponential distribution \[21, pp. 153-154\] in Figure 1 *(a)*, but verify that the fit is unsuitable to describe the large events.

### 2.2 Inverse square and cubic laws for the distribution of line length

We analyse the unconditional probability distribution of (axial) line length of $`36`$ cities in $`14`$ different countries (see Table 1). In our analysis we use the rank-order technique . To interpret the apparently unsystematic data in Figure 2(a) effectively, it is instructive to scale the data. Since the rank ranges between $`1`$ and $`\mathrm{max}\left(rank_i\right)`$, we define a scaled relative rank for city $`i`$, $`\widehat{r_{i,j}}rank_j/\mathrm{max}\left(rank_j\right)`$ . Similarly, for the ordinate, it is useful to define a scaled line length by $`\widehat{l_{i,j}}=l_j/l_j`$ .

As shown in Figure 2(d), there is relatively good collapse of the data sets onto two master curves for $`28`$ of the $`36`$ cities under study (the cities plot in red and blue). The other $`8`$ cities do not collapse clearly onto a single curve (see Figure 2(b)). Figure 2(c) is a plot of the exponents from a least-squares fit to the data of Figure 2(b) for $`\mathrm{log}\widehat{l_{i,j}}>0.2`$, where the data are visually the most linear. The relatively small error bars ($`95\%`$ confidence bounds) in Figure 2(c) corroborate our choice of the cutoff. The fits on the rank order plot lead to straight lines with slope $`1/\alpha `$, which suggest that the line length probability density may have a power-law tail, $`P\left(l\right)l^{1\alpha }`$ with exponents close to $`\alpha _B2`$ (cities in blue) or $`\alpha _R3`$ (cities in red). The inverse square and cubic laws have diverging higher moments (larger than $`2`$ and $`3,`$ respectively) and are not stable distributions.

Figure 3 is a plot of several axial maps, where we only plot lines with $`\mathrm{log}\widehat{l_{i,j}}>0.2`$ (the range of data used in the least-squares fit of Figure 2). We suggest that urban growth can be regarded as a process where axial lines are added mainly to the urban periphery (this could, in principle, be monitored through remote sensing techniques for cities undergoing rapid urbanization) and modelled by $`N`$ walkers, which jump along the corresponding $`N`$ axial lines, extending the city. More needs to be known on the distribution of the walkers’ waiting time to correctly model the dynamics of urban growth. Nevertheless, the walkers would generate a non-stable process, as the exponents $`\alpha _B`$ and $`\alpha _R`$ are, apparently, outside of the Lévy stable region.

## 3 Discussion

We have found that the length of urban open space structures displays universal features, largely independent of city size, and is self-similar across morphologically relevant ranges of scales ($`2`$ orders of magnitude) with exponents $`\alpha _B2`$ (cities in blue) and $`\alpha _R3`$ (cities in red). Our results are unexpected as two universality classes appear for a wide range of cities. The power-law tails of the pdfs support the hypothesis that urban space has a fractal structure , but the parallelism to a Sierpinski gasket  may be too simple for an accurate description.

Our findings show that it is important to model in detail the open space geometry of urban aggregates. They also support the hypothesis that it is more useful to model urban morphology as random rather than as the outcome of rational decisions, as previously suggested .

Cities with exponents $`\alpha _B2`$ (cities in blue) display open space alignments which cross the whole structure, whilst cities with $`\alpha _R3`$ (cities in red) tend not to. We propose that the large scale linear structures of the two classes of cities can be explained by two distinct non-stable Lévy processes, where the walker’s jump has a tail that goes as $`\alpha _B2`$ (cities in blue) or $`\alpha _R3`$ (cities in red). The walker’s jumps are much larger for a process with $`\alpha _B2`$ (cities in blue) than for a process with $`\alpha _R3`$ (cities in red), leading to the dominance of global geometric structures for the former as opposed to local geometric structures for the latter. The probability of a large axial line occurring for cities with $`\alpha _G<2`$ (subset of cities in green) is larger than for the other cities (see, e.g., Las Vegas in Figure 2). We suggest that cities in green with $`2<\alpha _G<3`$ have mixed influences of global and local structures and that cities in green with $`\alpha _G<2`$ have very strong patterns of global geometric structures. We note that the Lévy processes which, we propose, imprint the dominant global geometric structures of a city may appear at a late stage of urban growth, but change the city structure drastically (e.g. by generating long axial transportation routes).

Physicists have found exponents $`\alpha 3`$ when studying the distribution of normalized returns in financial markets, both for individual companies  and for market indexes . Our results for the class of cities plot in red on Figure 2 is reminiscent of these studies. We believe that parallels between urban growth and finance may not be too far fetched, as both processes seem to be largely dominated by geometric phenomena . Indeed, there may be similarities between the dynamics of price fluctuations and urban growth, and we propose that axial lines may be seen as the urban equivalent of economic returns.

As more data becomes available through remote sensing, quantitative analyses should provide an improved view of the spatio-temporal dynamics of urban growth, particularly in squatter settlements, where time-scales for growth are much shorter than in conventional cities and one could hope to model growth against observed data. Further studies are still required, but it seems that the impact of local controls on growth (e.g. the green belt policy for London) is, at most, spatially localized. Indeed, at a ’macro’ level, cities display a surprising degree of universality.

## 4 Acknowledgments

We acknowledge financial support from Grant EPSRC GR/N21376/01. We thank Prof. Bill Hillier for most useful discussions and encouragement, Prof. Mike Batty for suggestions on the rank-order plots and Shinichi Iida for invaluable insight on the axial map of Tokyo.


text:
# Paradoxical aspects of the kinetic equations

## Abstract

Two paradoxical aspects of the prevailing kinetic equations are presented. One is related to the usual understanding of distribution function and the other to the usual understanding of the phase space. With help of simple counterexamples and direct analyses, involved paradoxes manifest themselves.

Dept. of Physics, Beijing University of Aeronautics

and Astronautics, Beijing 100083, PRC

Email: cychen@buaa.edu.cn

PACS number: 51.10.+y.

When the Boltzmann equation, thought of as the first equation of the prevailing kinetic theory, came out, strong doubts arose, of which many were related to the fact that whereas Newton’s equations themselves were time-reversible the Boltzmann equation, as a consequence of Newton’s equations, was time- irreversible. As time went by, particularly after the “rigorous” BBGKY theory was formulated in the middle of the last century\[2$``$8\], the philosophical concern of the time-reversal paradox gradually faded out, and it was believed that the ultimate understanding of the related issues had been completely established, at least in the regime of classical mechanics.

However, relatively recent developments of mathematics and physics seem to have brought new elements into the picture. In particular, the studies of fractals reveal that there can, at least in the mathematical sense, exist functions that have structures of self-similarity even at the infinitesimal level. Such functions are intrinsically discontinuous and cannot be described by usual differential apparatuses. Along this line, an increasing number of scientists are surmising that if any similar structures are found in realistic gases, some of conventional concepts in the standard theory need to be revised significantly.

In connection with this, we wish to present and discuss two relatively unknown aspects of the standard framework of kinetic theory. Firstly, it will be shown that a realistic gas contains, almost always, a significant amount of particles whose distribution function does not keep invariant in the six-dimensional phase space and cannot be regarded as a continuous one. Secondly, it will be illustrated that there are inherent difficulties in formulating the particles that enter and leave an infinitesimal phase volume element during an infinitesimal time, which suggests in the sense that the phase space is more sophisticated than the customary thought assumes. By removing abstractness of the matters, all related paradoxes become surprisingly simple and straightforward. We are now convinced that if those paradoxes had been unveiled at the very berginning, kinetic theory would have been renewed several times.

The basic core of the classical kinetic theory says that if particles of a gas do not interact with each other, the distribution function describing them satisfies the collisionless Boltzmann equation

$$\frac{f}{t}+𝐯\frac{f}{𝐫}+\frac{𝐅}{m}\frac{f}{𝐯}=0,$$
(1)
which is, according to well-known textbooks, equivalent to the path-invariance of collisionless distribution function

$$\frac{df}{dt}|_{𝐫(t),𝐯(t)}=0,$$
(2)
where the subindexes $`𝐫(t)`$ and $`𝐯(t)`$ imply that that the differentiation is taken along a particle’s path. Conceptually speaking, equation (2) can be interpreted as saying that such a gas is “incompressible” in the phase space. The standard theory further states that, if collisions between particles are more than negligible, a certain type of collisional operator needs to be introduced to replace the zero term in (1). After having the collisional operator due to Boltzmann himself, equation (1) becomes

$$\frac{f}{t}+𝐯\frac{f}{𝐫}+\frac{𝐅}{m}\frac{f}{𝐯}=_{𝐯_1,\mathrm{\Omega }}[f(𝐯^{})f(𝐯_1^{})f(𝐯)f(𝐯_1)]u\sigma (\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑𝐯_1,$$
(3)
where the first term on the right side represents particles entering the unit phase volume element and the second term particles leaving the unit phase volume element.

Although the formulation briefed in the last paragraph seems stringent and has been accepted unanimously in the community, it is not truly sound. Ironically enough, even a simple glance at the form of the Boltzmann equation (3) offers intriguing things to ponder. On the left side there is a symmetry between $`𝐫`$ and $`𝐯`$ in terms of differentiation operation; whereas on the right side the position vector $`𝐫`$ serves as an inactive “parameter” and all the integral operations are performed in the velocity space. (The solid angle $`\mathrm{\Omega }`$ is defined in terms of the velocities.) This disparity, while seeming a bit too schematic to be fully convincing, should serve a motivation for us to investigate the whole issue more carefully and more thoroughly.

We first look at whether or not the collisionless Boltzmann equation (1) and the path-invariance theorem (2) make sense as they intend to.

In the mathematical sense, the picture provided by the path-invariance is quite clear and simple. Consider a sufficiently small volume element

$$\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z\mathrm{\Delta }v_x\mathrm{\Delta }v_y\mathrm{\Delta }v_z$$
(4)
moving together with a certain particle in the phase space; the theorem asserts that the particle density within the element (in number herein) does not increase or decrease. In well-known textbooks, the theorem is “proved” by applying the rigorous Jacobian approach\[11$``$13\]. Here, for the purpose of this paper, we wish to illustrate the theorem and its derivation in an intuitive way. Fig. 1 describes what takes place in the two-dimensional $`x`$-$`v_x`$ subspace. As time passes from $`t=0`$ to $`t=T`$, the particles distributed inside the rectangle $`\mathrm{\Delta }x\mathrm{\Delta }v_x`$ in Fig. 1a will be distributed inside the parallelogram in Fig. 1b. By excluding all external forces (just for simplicity), it is obvious that the area of the parallelogram is equal to that of the rectangle and the average particle density within the moving volume element keeps invariant. On the understanding that particles distribute continuously within $`\mathrm{\Delta }x\mathrm{\Delta }v_x`$ and the size of $`\mathrm{\Delta }x\mathrm{\Delta }v_x`$ can shrink to zero, the invariance of the average particle density can certainly be interpreted as the path-invariance of distribution function. The discussion above, though simplified somewhat, reflects the essence of the full theorem accurately.

It seems, at this stage, that the validity of the path-invariance (2), together with the picture of incompressible fluid in the phase space, is so solid and so clear that we should discuss it no more. But, somehow, the truth is not that simple: there exist many cases in that gases behave very differently. To get an immediate idea about such behavior, let’s look at the special arrangement shown in Fig. 2, where particles with definite velocity $`𝐯`$ strike a convex solid surface and then get reflected from it elastically. (Later on, more general models about collisions between particles and boundaries will be adopted and examined.) Following a moving particle and counting particles in a definite position range $`d𝐫dxdydz`$ and in a definite velocity range $`d𝐯dv_xdv_ydv_z`$, we find out a clear-cut fact that the reflected particles simply obey

$$\frac{df}{dt}|_{𝐫(t),𝐯(t)}<0\mathrm{and}\frac{f}{t}+𝐯\frac{f}{𝐫}+\frac{𝐅}{m}\frac{f}{𝐯}<0.$$
(5)
That is to say, these particles diverge in the phase space. In a similar way, we can readily see that if the incident particles fall upon a concave solid surface, the reflected particles will, in some region, obey

$$\frac{df}{dt}|_{𝐫(t),𝐯(t)}>0.$$
(6)
This tells us that the particles converge.

A sharp question arises. How can the path-invariance theorem, considered as the very core of all the kinetic equations, possibly suffer from such simple and direct counterexamples? By reviewing the derivation of the Boltzmann equation, we realize that equations (1) and (2) hold only for perfectly continuous distribution functions, while the diverging and converging particles presented above are related to none of them.

To see the point more vividly, let’s go back to the two-dimensional phase space, namely the $`x`$-$`v_x`$ space, and investigate the evolution of the particles marked with the dotted curve in Fig. 2. In the schematic sense, we may say that at $`t=0`$ all these particles are distributed along the diagonal of $`\mathrm{\Delta }x\mathrm{\Delta }v_x`$ shown in Fig. 3a, rather than in $`\mathrm{\Delta }x\mathrm{\Delta }v_x`$ uniformly. After a while, at $`t=T`$, this diagonal is stretched and becomes longer as shown in Fig. 3b. If we set up a small, but definite, phase volume element and let it move together with one of the particles, we will certainly find that the particle density within it decreases. (For discontinuous distributions, if we allow the shape and size of measuring volume element to vary, we can get any value, drastically from zero to infinity, for the particle density.)

Now, it is in order to examine how a realistic boundary surface $`S`$ “reflects” realistic particles that come from all directions and with different speeds. Experimental facts inform us that the reflection cannot be truly elastic and stochastic and dissipative effects must play a certain role. To express particles “produced” by such surface, it is proper (particularly if the gas is a rarefied one) to define the instant emission rate $`\sigma `$ on an infinitesimal surface $`dS`$ in such a way that

$$\sigma dtdSdvd\mathrm{\Omega }$$
(7)
represents the number of particles ejected by the surface element in the speed range $`dv`$ and in the solid angle range $`d\mathrm{\Omega }`$ during the time interval $`dt`$. By allowing the emission rate $`\sigma `$ to be described in probability and to have certain dependence on velocities of incident particles, it can be said that the stochastic and dissipative nature of the reflection has been included. Then, divide $`S`$ into $`N`$ surface elements. Referring to Fig. 4, we find that for the $`i`$th surface element $`(\mathrm{\Delta }S)_i`$ the reflected particles are like ones emitted from a “point particle source” at $`𝐫_{0i}`$, the position vector of $`(\mathrm{\Delta }S)_i`$, and the relevant distribution function at a point $`𝐫`$ in the reflection region takes the form

$$f(𝐫,v,\mathrm{\Omega })=\frac{\sigma (\mathrm{\Delta }S)_i}{|𝐫𝐫_{0i}|^2v^3}V_i(v)\delta (\mathrm{\Omega }\mathrm{\Omega }_{𝐫𝐫_{0i}}),$$
(8)
where $`V_i(v)`$ is a certain function of $`v`$ and $`\delta (\mathrm{\Omega }\mathrm{\Omega }_{𝐫𝐫_{0i}})`$ is the $`\delta `$-function defined on the solid angle in the velocity space. It is very obvious that, regardless of the forms of $`\sigma `$ and $`V(v)`$, this distribution function diverges in the phase space. Another interesting point about the distribution function is that it is perfectly continuous in terms of $`𝐫`$ and $`v`$ and it is like a function defined on a single point in terms of $`\mathrm{\Omega }`$. Based on this observation, we are tempted to say that the function is on a variable domain of $`4+ϵ`$ dimensions. The total distribution function at $`𝐫`$ associated with all reflected particles from $`S`$ is

$$f(𝐫,v,\mathrm{\Omega })=\underset{i=1}{\overset{N}{}}\frac{\sigma (\mathrm{\Delta }S)_i}{|𝐫𝐫_{0i}|^2v^3}V_i(v)\delta (\mathrm{\Omega }\mathrm{\Omega }_{𝐫𝐫_{0i}}).$$
(9)
Though the distribution function expressed by (9) can, as $`N\mathrm{}`$, be regarded as a sum of an infinitely large number of $`(4+ϵ)`$dimensional functions, it is not a continuous function defined in the phase space. Actually, with help of this expression, we can analytically or numerically prove that the particle number within a given moving volume element $`d𝐫d𝐯`$ is capable of decreasing, increasing or keeping invariant. As a limiting case, we may assume $`S`$ to be relatively small, or the distance $`|𝐫𝐫_{0i}|`$ to be relatively large, and find that expression (9) is reduced to expression (8) and the particles related to it always diverge in the phase space.

The above investigation, though formally simple, provides us with a very new picture on gas dynamics. Rather than as a continuous medium or an incompressible fluid in the phase space, a gas should be considered as a special collection of discrete particles, whose distribution function can change from continuous one to discontinuous one, as well as from discontinuous one to continuous one. Another notable fact, which may very much interest scientists carrying out practical studies, is that the changeovers aforementioned occur dramatically near interfaces between fluids and solid boundaries.

We now turn our attention to the validity of the collisional operator on the right side of (3). A widely accepted concept is that in order to formulate collisional effects we are supposed to focus ourselves on a fixed six-dimensional phase volume element and study how particles leave and enter the element due to collisions. In what follows, it will be shown that such concept is, much to our surprise, deceptive.

Firstly, particles leaving a phase volume element $`d𝐫d𝐯`$ during a time interval $`dt`$ are of interest. According to the standard theory, if a particle collides with another within the volume element $`d𝐫d𝐯`$ during $`dt`$, it should be considered as one that leaves the volume element during the time interval because of the velocity change caused by the collision (see Fig. 5). This intuitive, seemingly very reasonable, picture can be challenged in the following way. In deriving kinetic equations, a necessary step is to let $`dt`$, $`d𝐫`$ and $`d𝐯`$ approach zero independently. If it is assumed, in the limiting processes, that the length scale $`|d𝐫|`$ is much smaller than $`|𝐯dt|`$, then virtually all the particles, initially within $`d𝐫d𝐯`$, will leave $`d𝐫d𝐯`$ at the end of $`dt`$, irrespective of suffering collisions or not. If we still want to say that the standard consideration, concerning how many particles stay inside $`d𝐫d𝐯`$ without involving collisions and how many particles get out of $`d𝐫d𝐯`$ with collisions involved, holds its significance, we have no choice but to assume that $`|d𝐫|>>|𝐯dt|`$. An unfortunate fact is that no sound reason can be found out for that we can prefer this assumption to its converse.

Secondly, particles entering a phase volume element $`d𝐫d𝐯`$ during $`dt`$ are under examination. In the standard treatment, two beams with velocities $`𝐯^{}`$ and $`𝐯_1^{}`$ are assumed to collide within a volume element $`d𝐫`$ during $`dt`$ and to give contribution to the particles expressed by $`f(𝐫,𝐯,t)d𝐫d𝐯`$, as shown in Fig. 6. Ironically enough, this treatment involves not one but many paradoxes. As one thing, particles produced by collisions in a small spatial region, like ones emitted from a point particle source, will diverge in the phase space and cannot be treated as an ordinary contribution to $`f(t,𝐫,𝐯)d𝐫d𝐯`$. As another thing, we again need to let $`dt`$, $`d𝐫`$ and $`d𝐯`$ approach zero. If $`|d𝐫|<<|𝐯dt|`$, all produced particles will “instantly” leave $`d𝐫d𝐯`$ and only an insignificant fraction of them can be treated as a contribution to $`f(𝐫,𝐯,t)d𝐫d𝐯`$. See Ref. 16 to get more paradoxes.

In summary, two paradoxical aspects of the standard kinetic theory have been presented. The first aspect is related to the distribution function. A tacit assumption of the existing kinetic equations is that distribution functions, though describing discrete particles, must be mentally and practically continuous. The discussion of this paper, however, shows that distribution functions of realistic gases have, in general, complex local structures and they cannot be described by continuous distribution functions. The second aspect is related to the phase space. A usual concept in the standard theory is that the position space and the velocity space can be separated mentally and practically. The discussion of this paper, however, shows that whenever we investigate the time dynamics of particles in the velocity space we should keep an eye on what takes place in the position space, and vice versa.

A variety of fundamental questions can be raised, of which many are beyond the scope of this brief paper. In some of our recent works, we make more analyses and put forward alternative approaches. With help of a development in quantum mechanics, some of the discussion in this paper can also be extended to the regime of quantum statistical physics.

This paper is supported by School of Science, BUAA, PRC and by Education Ministry, PRC.

Fig. 1: A moving phase volume element (a) at an initial time $`t=0`$ and (b) at a later time $`t=T`$.

Fig. 2: Schematic of reflected particles from a solid surface. It is interesting to note that whether or not the involved collisions are elastic will not alter this picture significantly.

Fig. 3: Schematic of how the particles marked in Fig. 2 spread in the $`xv_x`$ space. (a) These particles are distributed along one diagonal of the rectangle $`\mathrm{\Delta }x\mathrm{\Delta }v_x`$ at $`t=0`$, and (b) the diagonal is stretched at $`t=T`$.

Fig. 4: Particles reflected from a small surface element of realistic boundary.

Fig. 5: A particle involving a collision in a phase volume element.

Fig. 6: Two beams of particles collide with each other and produce particles with velocity $`𝐯`$.


text:
# Reply to C. Tsallis’s “Comment on Critique of q-entropy for thermal statistics by M. Nauenberg”

##

In the abstract of his paper , and in a section entitled “About thermal contact between systems with different values of $`q`$ and the $`0^{th}`$ principle of thermodynamics”, Tsallis acknowledges that the essential point of my critique  is that the q-entropy formalism does not give a prescription to obtain the equilibrium temperature for such systems, as is required by this principle. I paraphrased this failure of the formalism by stating that a conventional thermometer satisfying Boltzmann-Gibbs (BG) statistics could not measure the temperature of a q-entropic system with $`q1`$, and concluded that “the laws of thermodynamics would therefore fail to have general validiy” unless all systems have the same value of $`q`$. Since there are systems in thermal equilibrium for which $`q=1`$, corresponding to BG statistics, my conclusion implies that all systems in thermal equilibrium must have $`q=1`$, and thus correspond to BG statistics. But rather than point out supposed flaws in my criticism, Tsallis responded in  by presenting the results of a molecular dynamics simulation of coupled rotators with long range interactions, which are known to exhibit a quasistationary state below the transition temperature . Tsallis claimed that this state has “non Boltzmannian statistics” without, however, providing any evidence, and then purported to show that a thermometer constructed out of rotators with nearest-neighbor interactions that satisfy BG statistics measures the temperature of this state. Tsallis’s temperature is twice the mean kinetic energy of the rotators, but this definition is only justified if the model satisfies BG statistics. However, Tsallis admits that this is not the case for the rotator-model quasistationary state. Moreover, the general thermodynamic definition of temperature (which is also introduced in the q- entropy formalism) is the derivative of the entropy with respect to the energy, but neither the entropy nor this derivative was evaluated by Tsallis’s in this simulation. The details of this simulation were discussed in a separate paper , but not a shred of evidence was presented that this quasistationary state satisfies the statistics associated with q-entropy. In conclusion, Tsallis’s Comments  cannot in any way be regarded as a rebuttal to my main criticism of the q-entropy formalism in reference .

As an example of one of the many unphysical results which have been obtained during the past decade by applying the q-entropy formalism to problems in thermal statistics, I discussed at some length in my critique the case of black-body radiation. In particular, I pointed out that this formalism leads to a violation of the well-known $`T^4`$ Stefan-Boltzmann law. In his Comments  Tsallis responded that it can be “trivially shown that the $`T^4`$ proportionality law remains the same for all energy statistical distribution (hence not only the BG ones) as long as…for photons the distribution depends on the light frequency $`\nu `$ and the appropriate temperature $`T`$, only through \[the ratio\] $`\nu /T`$”. But as I explained in my critique, for a nonextensive theory both the energy and the entropy density depend not only on $`\nu /T`$, but also (by definition) on the volume $`V`$ of the cavity holding the black-body radiation. Otherwise the theory would be extensive. This volume dependence implies that these distributions must be functions not only of a dimensionless ratio $`h\nu /k_BT`$, but also of another dimensionles ratio $`V/a(T)^3`$, where $`a(T)`$ is a length parameter, characteristic of the system, which depends on the temperature $`T`$. By dimensional arguments or from ordinary statistical mechanics we know that $`a(T)hc/k_BT`$ is the mean thermal wavelength of the photons.

To show the fallacy in numerous publications of Tsallis and others (for references, see ) who argue that a nonextensive formalism lead to the $`T^4`$ Stefan-Boltzmann law, I give here the details of a thermodynamic derivation of this law in such a formalism. Suppose that $`U/V=f(T)V^{\delta /3}`$ and $`S/V=g(T)V^{\delta /3}`$, where $`U`$ is the energy, $`S`$ is the black-body entropy in a cavity of volume $`V`$, $`f(T)`$ and $`g(T)`$ are unknown functions of $`T`$, and the exponent $`\delta `$ is a constant which characterizes a nonextensive formalism. Inserting these expressions into the well known thermodynamic relation

$$dU=TdSpdV,$$
(1)
where $`p`$ is Maxwell’s isotropic electromagnetic radiation pressure $`p=(1/3)U/V`$, one obtains the equations

$$\frac{df(T)}{dT}=T\frac{dg(T)}{dT}$$
(2)
and

$$(4+\delta )f(T)=(3+\delta )Tg(T).$$
(3)
These two equations can be readily solved to give

$$f(T)T^{4+\delta }$$
(4)
and

$$g(T)T^{3+\delta }.$$
(5)
Hence, contrary to Tsallis’s assertions, this general proof demonstrates that to obtain the Stefan-Boltzmann law, both the energy and the entropy of black-body radiation must be extensive, i.e. $`\delta =0`$. Indeed, Boltzmann made this assumption in his original 1884 derivation of this law . I would also like to point out that already in 1916, Einstein showed that Planck’s formula for the black-body energy density, which leads to the Stefan Boltzmann law, can be derived from his quantum transition probabilities for the emission and absorption of radiation . These transition probabilities lead to the conventional Bose-Einstein statistics , and not to the modified versions which have been obtained by various applications of the q-entropy formalism to black-body radiation (for references, see ). Hence these applications are also inconsistent with the laws of quantum electrodynamics.

As I have already indicated in the title of my paper, my critique of the q-entropy formalism was confined only to its applications to systems in thermal equilibrium. Tsallis, however, devoted a major part of his Comments  to applications of this formalism to non-thermal problems such as the logistic map, the standard map, growth models, etc. These applications, however, are not relevant to my critique, and therefore I will not discuss them here.