text:
# Spectral properties of the overlap Dirac operator in QCD

## 1 Introduction

*Chiral perturbation theory* is a powerful tool to analyze QCD at low energy. In this context, the pions as the lightest hadrons carry the relevant degrees of freedom. They are considered as quasi-Nambu-Goldstone bosons obtained from chiral symmetry breaking — their mass is provided by the masses of the $`u`$ and $`d`$ quark, which are small compared to $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. One then constructs an effective theory for the expansion in the pion mass and momenta, which is strongly constrained by the requirements of chiral symmetry.

In particular the so-called *$`ϵ`$-expansion* is designed to describe the system in a volume which is so small that it cannot even include a pion Compton wave length . In this situation the rôle of the topological sectors of the gauge field is far more important than it is the case in a large volume .

In view of the comparison to lattice data, it is an important virtue that a finite and even small volume can be interpreted directly — in contrast to most other lattice simulations one does not need to extrapolate to large volumes. Chiral perturbation theory in the $`ϵ`$-expansion describes the finite volume and quark mass dependence of various quantities such as the scalar condensate or 2-point correlation functions. These dependencies are parametrized by *the same* low energy constants as in the effective chiral lagrangian in the infinite volume. Therefore the unphysical small volume can provide results for physical parameters. The finite size effects allow for a determination of the low energy constants in the effective chiral lagrangian. Examples of this procedure are given in ref.  for the O(4) symmetric non-linear $`\sigma `$-model, and in ref.  for (quenched) QCD.

In this paper we are going to investigate *Ginsparg-Wilson fermions* , which allow for simulations at small fermion masses. In particular, the evaluation of the corresponding lattice Dirac operator can be done directly in the chiral limit, where the Ginsparg-Wilson fermion has an exact (lattice modified) chiral symmetry at finite lattice spacing . Moreover, since the Ginsparg-Wilson fermions have exact zero modes, we can use them to define the topological charge of the lattice gauge field  by means of the Atiyah-Singer Index Theorem .

A particularly simple solution of the Ginsparg-Wilson relation is obtained if the Wilson operator is inserted into the overlap formula . As long as the gauge coupling is not too strong, such a lattice Dirac operator is local  in the sense of an exponential localization, hence it is conceptually correct in view of the continuum limit. However, its simulation is tedious: for the time being, only quenched QCD simulations are feasible. Since also the quenched chiral perturbation theory in the $`ϵ`$-regime has been analyzed , the simulation data can be compared to these analytical predictions.

Even without considering hadrons, there are interesting predictions for the spectrum of the Dirac operator which are based on Random Matrix Theory (RMT), following the ideas of ref. . Chiral RMT simplifies QCD to a gaussian distribution of the elements in the fermion matrix, with the global symmetries of the QCD Dirac operator. For a review on this field, see for instance ref. . We only mention that although on the lagrangian level RMT is equivalent to chiral perturbation theory at lowest order in the $`ϵ`$-expansion, the predictions for the eigenvalue distributions involve additional assumptions. Thus numerical simulations can provide a test of the predictions from RMT beyond chiral perturbation theory. Even if the fermion mass vanishes, there is still a *lower limit* for the volume where RMT applies; in some steps, the latter assumes the volume (mathematically) to be infinite. Theoretically one often requires the temperature $`T`$ to be well below the Thouless energy $`TE_c=F_\pi ^2/(\sqrt{V}\mathrm{\Sigma })`$. However, it is not predicted explicitly at which lower bound for the volume RMT really collapses.

There are explicit predictions for the statistical distribution of the low-lying eigenvalues  which can be confronted with lattice data. We refer here to predictions in full QCD, although one might be worried that quenching causes logarithmic corrections . In this comparison, the chiral condensate $`\mathrm{\Sigma }`$ enters as a free parameter, hence, in turn, the fit to the functions predicted by RMT provide a value for $`\mathrm{\Sigma }`$. Another theoretical conjecture refers to the unfolded level spacing distribution which is not sensitive to the energy scale, but which allows for the inclusion of bulk eigenvalues.

Such studies were performed before, using staggered fermions in QCD , as well as Ginsparg-Wilson fermions in the Schwinger model , in QED  and in QCD on $`4^4`$ lattices . Here we are going to present data for the overlap Dirac operator spectrum in quenched QCD on larger lattices than used before for Ginsparg-Wilson fermions.

## 2 The spectrum of the overlap Dirac operator compared to the RMT predictions

### 2.1 Theoretical aspects

In QCD with two massless quark flavors, the chiral symmetry breaks spontaneously as

$$SU(2)_LSU(2)_RSU(2)_{L+R},$$
(1)
which generates three Nambu-Goldstone bosons. If we deal with light $`u`$ and $`d`$ quarks, these bosons pick up small masses because the quark mass breaks the chiral symmetry explicitly. The resulting quasi-Nambu-Goldstone bosons are identified with the pions. They can be represented by $`SU(2)`$ matrix fields $`U(x)`$. At low energy they are effectively described by chiral perturbation theory. To the leading order, the corresponding lagrangian reads

$$[U]=\frac{F_\pi ^2}{4}Tr\left[_\mu U(x)_\mu U^{}(x)\right]\frac{\mathrm{\Sigma }m}{2}Tr\left[e^{i\theta /2}U(x)+e^{i\theta /2}U^{}(x)\right].$$
(2)
Here $`m`$ is the quark mass (for simplicity we assume it to coincide for the two flavors) and $`\theta `$ is the vacuum angle. The coupling constants in leading order are the pion decay constant $`F_\pi `$ and the chiral condensate $`\mathrm{\Sigma }`$.

We put the system in a periodic box of size $`L^4`$. Then the $`ϵ`$-regime is characterized by the condition

$$\frac{1}{m_\pi }L\frac{1}{4\pi F_\pi },m_\pi =\text{pion mass}.$$
(3)
In this regime, RMT can be applied to the Dirac operator $`D`$. In the continuum it leads to explicit formulae for the statistical distribution of the low eigenvalues, i.e. for the *“microscopic regime”*. There the momenta are counted as $`O(1/L)`$. One introduces the spectral density

$$\rho (\lambda )=\underset{n}{}\delta (\lambda \lambda _n),$$
(4)
where the sum runs over all modes, and $`\lambda _n`$ are the eigenvalues of the operator $`iD`$. It is favorable to use the dimensionless variable

$$z=\lambda V\mathrm{\Sigma },$$
(5)
which is of $`O(1)`$. In these terms, the RMT prediction addresses the microscopic spectral density

$$\rho _s(z)=_V\mathrm{}^{lim}\frac{1}{V\mathrm{\Sigma }}\rho \left(\frac{z}{V\mathrm{\Sigma }}\right),$$
(6)
which can be decomposed into the contributions of different topological sectors,

$$\rho _s(z)=\underset{\nu =\mathrm{}}{\overset{\mathrm{}}{}}\rho _s^{(\nu )}(z),\nu =\text{topological charge}.$$
(7)
Refs.  present expressions for the leading contributions $`\rho _n^{(\nu )}(z)`$ in

$$\rho _s^{(\nu )}(z)=\underset{n}{}{}_{}{}^{}\rho _{n}^{(\nu )}(z),$$
(8)
where the sum runs over the non-zero modes. The predictions for the lowest eigenvalue at $`|\nu |=0,`$ 1 and 2 are depicted in figure 1. We see in particular that the density peak moves to larger values of $`z`$ as $`|\nu |`$ increases.

### 2.2 Lattice results in the microscopic regime

In order to verify the viability of the $`ϵ`$-expansion and of chiral RMT, it is of interest to compare these predictions to the corresponding eigenvalue distributions obtained from lattice simulations. The Wilson Dirac operator seems hardly suitable in this respect; because of the mass renormalization it would be highly problematic to identify the zero modes and the relative values of the remaining modes. In addition, simulations at small quark masses are problematic due to exceptional configurations.

Staggered fermions do not suffer from additive mass renormalization because they have an exact remnant chiral symmetry $`U(1)U(1)`$. Indeed such comparisons exist. The staggered fermion spectrum agrees well with the RMT prediction in the sector $`\nu =0`$. However, it turns out that all other sectors, $`\nu 0`$, yield the same distributions, in particular the same histogram for $`\rho _1^{(\nu )}(z)`$ in QCD . It seems that staggered fermions are generally insensitive to the topology, and therefore not adequate for this purpose, at least for moderate and strong gauge coupling.

Another type of lattice fermions without additive mass renormalization is the Ginsparg-Wilson fermions. Their lattice Dirac operator $`D`$ fulfills the Ginsparg-Wilson relation (GWR), which reads (in its simplest form)

$$D\gamma _5+\gamma _5D=\frac{a}{\mu }D\gamma _5D,$$
(9)
i.e. the continuum condition for the right-hand-side to vanish is relaxed to a term of $`O(a)`$, where $`a`$ is the lattice spacing. $`\mu `$ is a mass parameter, see below. The GWR implies that $`D^1\gamma _5+\gamma _5D^1=a\gamma _5/\mu `$ is local, therefore the poles in the propagator are not shifted away from zero. Thus the Ginsparg-Wilson fermions have exact zero modes, which occur with positive or negative chirality. We now define the topological charge $`\nu `$ simply by the index

$$\nu =n_+n_{},$$
(10)
where $`n_+`$ ($`n_{}`$) is the number of positive (negative) chiral zero modes. Here one adapts the continuum Index Theorem and uses it to define the topological charge of the lattice gauge configurations . The applicability of chiral RMT to Ginsparg-Wilson fermions has been discussed theoretically in ref. .

Many lattice Dirac operators, such as the Wilson operator $`D_W`$, obey $`\gamma _5`$ Hermiticity, $`D_W^{}=\gamma _5D_W\gamma _5`$. This allows for a simple solution of the GWR by inserting $`D_W`$ (say with Wilson parameter 1) into the overlap formula ,

$$D_{\mathrm{ov}}=\frac{\mu }{a}\left[1+A/\sqrt{A^{}A}\right],A=aD_W\mu .$$
(11)
$`D_{\mathrm{ov}}`$ is a solution to condition (9) and it is $`\gamma _5`$-hermitean as well. (There are alternatives to replace $`D_W`$ by lattice Dirac operators that are more involved but lead to improved properties, as suggested in refs. .) $`\mu `$ represents a negative mass of the Wilson fermion, which can be chosen in some interval as long as the gauge fields are smooth. In the free case $`\mu =1`$ is optimal with respect to locality, but at $`\beta =6`$ and $`\beta =5.85`$ we move to $`\mu =1.4`$ resp. $`\mu =1.6`$ to compensate the mass renormalization of the Wilson fermion.

Note that the operator $`A/\sqrt{A^{}A}`$ is unitary, hence the spectrum of $`D_{\mathrm{ov}}`$ is located on a circle in the complex plane through zero, with center and radius $`\mu /a`$. In order to relate the eigenvalues found on this Ginsparg-Wilson circle to the continuum eigenvalues $`\lambda _n`$, we map the circle stereographically onto the imaginary axis. Requiring $`f(z)=z+O(z^2)`$ and $`f(2\mu )=\mathrm{}`$ singles out the Möbius transform

$$f(z)=\frac{z}{1z/(2\mu )}.$$
(12)
This mapping has been suggested before in ref. , in connection with the Leutwyler-Smilga sum rules .

The eigenvalues of Ginsparg-Wilson operators were compared to the RMT prediction for the Schwinger model in ref. , where the overlap operator (11) as well as a truncated fixed point Dirac operator were considered. In both cases the results agreed with the RMT formulae.

Such studies were also performed with the overlap operator in 4d QED  and again the predictions were successful within the statistical errors.

Finally QCD was considered, but only on small lattices of size $`4^4`$. Ref.  used the overlap operator at strong coupling of $`\beta =5.1`$. Ref.  applied a truncated fixed point action and obtained a decent agreement in a volume of $`V=(1.2\mathrm{fm})^4`$. However, when the physical lattice spacing is decreased so that the volume shrinks to $`(0.88\mathrm{fm})^4`$, the leading non-zero eigenvalue distributions of the different topological sectors are on top of each other, in contrast to the RMT prediction.

Our results are obtained by using again the overlap operator (11) in QCD with the standard plaquette gauge action. We have chosen a rather weak gauge coupling, $`\beta 5.85`$, so that the overlap formula can safely be applied to construct a Ginsparg-Wilson fermion.111From our point of view, it is not clear if the overlap fermion is well-controlled at stronger coupling. For instance, typical spectra at $`\beta =5`$ and even $`\beta =5.4`$ on $`4^4`$ lattices have their eigenvalues spread over a broad area , so that it is hardly possible to find a good value for $`\mu `$ which would split the (nearly) real eigenvalues into a small branch (to be mapped to the vicinity of 0), and a large branch (to be mapped onto the opposite arc of the Ginsparg-Wilson circle). Due to the computational effort required by the overlap operator we had to use the quenched approximation, as it was also the case in all previous studies mentioned above. We approximated the inverse square root by Chebyshev polynomials to an accuracy of $`10^{12}`$. Then we used the PARPACK routines to evaluate up to 100 eigenvalues of $`D_{\mathrm{ov}}`$. We focused on the eigenvalues with the least real parts resp. absolute values. Of course the non-zero eigenvalues occur in complex conjugate pairs, hence we only consider one sign for the imaginary part.

First we measured the leading non-zero eigenvalue on a $`12^4`$ lattice at $`\beta =6`$. Here the volume amounts to $`(1.12\mathrm{fm})^4`$. The result can be represented by a histogram, which could then be compared to figure 1. However, such a picture depends on the (arbitrary) choice of the bin size in the histogram. This can be avoided by plotting the “cumulative density”, see e.g. ref. , which sums up all the entries up to the considered value of $`z`$ — it is normalized so that the full number of entries corresponds to the cumulative density 1. This way to plot the data can be compared to the cumulative density according to the RMT prediction,

$$\rho _{1,c}^{(\nu )}(z)=\frac{_0^z\rho _1^{(\nu )}(z^{})𝑑z^{}}{_0^{\mathrm{}}\rho _1^{(\nu )}(z^{})𝑑z^{}}.$$
(13)
This comparison is done in figure 2 (top). We recognize a satisfactory agreement, especially for $`|\nu |=1`$, where we have the largest statistics, see table 1. In particular, we do see the effect that the peak of the density — resp. the interval of steepest ascent of the cumulative density — moves to larger values of $`z`$ for increasing topological charge. In this plot the chiral condensate $`\mathrm{\Sigma }`$ enters as the one free parameter, which sets the energy scale. Our plot is shown for the value

$$\mathrm{\Sigma }=(256\mathrm{MeV})^3,$$
(14)
which provides optimal agreement with the RMT curves.

Next we studied a $`10^4`$ lattice at $`\beta =5.85`$, which corresponds to a somewhat larger physical volume of $`V=(1.23\mathrm{fm})^4`$. Again the results for the leading non-zero eigenvalue are shown for the topological sectors $`|\nu |=0,1`$ and $`2`$ in figure 2 (bottom). Also here we find a reasonable agreement with the RMT predictions, and the optimal value of the chiral condensate is modified to $`\mathrm{\Sigma }=(253\mathrm{MeV})^3`$. RMT predicts a (logarithmically) increasing $`\mathrm{\Sigma }(V)`$ due to quenching . A conclusive verification of this behavior (beyond possible lattice artifacts) would require further simulations. In figure 3 we show our results for the cumulative density of the next-to-leading eigenvalue at $`|\nu |=0`$, $`1`$ and $`2`$, for both lattices discussed before. Here we do not find a convincing agreement with the RMT prediction in the first ref. of ; note that the relevant values of $`z`$ are larger compared to figure 2. Also the second non-zero eigenvalue moves to larger values of $`z`$ if $`|\nu |`$ increases.

Now we consider the *eigenvalue density* without focusing particularly on the first or on the second non-zero eigenvalue (though we still omit the zeros). Then we obtain for instance at $`|\nu |=1`$ the densities shown in figure 4. RMT predicts an oscillating behavior (see first ref. in ), which is also plotted for comparison. We find a good agreement roughly up to the second peak (as expected from figures 2 and 3). Then we are leaving the microscopic regime and turn to the bulk; for the latter the density is supposed to rise as $`\lambda ^3`$ resp. $`z^3`$. This is nicely confirmed by our data; an example is illustrated in figure 5.

Finally we consider a smaller volume: we use a $`8^4`$ lattice at $`\beta =5.85`$, which corresponds to $`V=(0.98\mathrm{fm})^4`$. Figure 6 shows that in this case there is a clear disagreement with the RMT prediction; in particular the sensitivity of the peak to $`|\nu |`$ is lost to a large extent (beyond very small values of $`z`$). This is in qualitative agreement with the results of ref.  on the $`4^4`$ lattice. The observation that the volume has to be $`V(1.2\mathrm{fm})^4`$ can be viewed as an empirical determination of the Thouless energy.

At this point, we would like to add that the above results tend to look similar for the leading eigenvalues of the hermitean overlap Dirac operator

$$H_{\mathrm{ov}}=\gamma _5D_{\mathrm{ov}}.$$
(15)
(Of course, in this case the stereographic mapping is not needed.) As an example we show the cumulative density of the first positive (non-zero) eigenvalue of $`H_{\mathrm{ov}}`$ at $`L=10`$, $`\beta =5.85`$ in figure 7.

For completeness we report our statistics in table 1 and summarize the results in table 2. In that table we also include an estimate for the topological susceptibility, where some configurations with higher charges contribute as well. We present the dimensionless quantity $`\chi _tr_0^4`$, where $`\chi _t=\nu ^2/V`$ is the susceptibility and $`r_0`$ is the Sommer parameter obtained from the static quark-antiquark potential . The results are well consistent with the values in the literature (for an overview, see ref. \[20, figure 18\]), although our statistics is clearly too small for that quantity. (Hence we do not try to estimate errors on $`\mathrm{\Sigma }`$ and $`\chi _tr_0^4`$.) Finally table 2 also gives the results of the Kolmogorov-Smirnov test  for the confidence level of the results obtained, if one assumes the RMT probability distribution to be correct. The corresponding number described the probability for the observed deviation from the theory to occur with the given statistics. This test is just designed for the cumulative density. The results suggest that the volume $`(1.12\mathrm{fm})^4`$ is still somewhat small, so that the topologies are still not fully separated, but we see that this process sets in at least for the first non-zero eigenvalue on volumes larger than that.

### 2.3 Unfolded distribution

We now turn to another way of comparing our lattice data to a conjecture from RMT. This new evaluation allows us to take all our non-zero eigenvalues into account (for one sign of the imaginary part, i.e. up to 50 for each configuration).

We build from all the eigenvalues of all configurations the “unfolded distribution” as described for instance in ref. 222A different, more systematic notion of unfolding is described in ref. . (for generalities, see ref. ). To this end, we first numerate all available non-zero eigenvalues with positive imaginary part in each configuration in ascending order, given by the angle in the Ginsparg-Wilson circle. Then we put all these eigenvalues from all configurations together and numerate them again in ascending order. Now we consider pairs of eigenvalues from the same configuration, which follow immediately one after the other in the original numeration. If they differ in the global numeration by $`k`$, then $`k/N_{\mathrm{conf}}`$ is the “unfolded level spacing”, where $`N_{\mathrm{conf}}`$ is the number of configurations involved.

In figure 8 we show our results for $`L=12`$, at $`\beta =6`$ and for $`L=8`$, at $`\beta =5.85`$. The histograms are compared to RMT conjectures for different groups and representations, and we can confirm a good agreement with the Wigner distribution predicted for $`SU(3)`$ in the fundamental representation. Considering also ref.  we conclude that this property holds over a large range of volumes. Since we include all eigenvalues here, our statistics is much larger than in the plots for the microscopic regime. Note, however, that this analysis is not sensitive to the topology any more.

Moreover the unfolded distribution only tests the correct symmetry group that is to be taken in RMT. There is no physics information such as values for the condensate or the pion decay constant, because one does not keep track of a physical scale. Nevertheless, the unfolded distribution provides a non-trivial test for RMT.

## 3 Conclusions

We performed a lattice study of quenched QCD using the overlap operator to test RMT predictions for the eigenvalue distributions of the Dirac operator. In this pilot study we worked at $`\beta `$ values large enough so that the overlap formula can safely be applied as a solution of the GWR. By using lattices of size $`8^4`$, $`10^4`$ and $`12^4`$ we were able to reach physical volumes which are large enough to test RMT.

In the microscopic regime of the very small eigenvalues we observe agreement with the predictions by RMT applied to the $`ϵ`$-regime of chiral perturbation theory, *if* the physical volume is large enough; $`V(1.2\mathrm{fm})^4`$ seems to be sufficient to capture at least the leading non-zero eigenvalue. On smaller physical volumes this agreement — and in particular the sensitivity to the topology — disappears.

If we consider the unfolded level spacing distribution, however, we obtain a good agreement with the RMT conjecture, even down to a small physical volume.

The eigenvalue distributions in figure 2 show interesting characteristic features. In topological charge sector zero the probability for a very low non-zero eigenvalue is non-negligible. Therefore such very low-lying modes will appear in simulations. This might have serious consequences for numerical studies when the quark mass is lowered too much. In particular, this behavior might render simulations in the $`ϵ`$-regime problematic. Indications of these difficulties were presented in ref. , where also an analysis of the influence of the lowest mode in the $`\nu =0`$ sector on the measurement of the chiral condensate was given. It was reported that the condensate is very hard to measure in topological charge sector zero, in agreement with the results presented here.

On the other hand, the situation is clearly better in a non-trivial topology. Eigenvalues close to zero are strongly suppressed and a much better signal can be expected.

Note that this phenomenon is not restricted to the quenched approximation considered here but holds also in the full theory. Thus, it can be expected that simulations close to physical values of the pion mass may become extremely expensive. Clearly, contact to chiral perturbation theory is then mandatory to extrapolate to physical pion masses.

### 3.1 Outlook

We still want to proceed to larger physical volumes. This will hopefully lead to further precision in the agreement for the leading non-zero eigenvalues in the different topological sectors. The larger volume might provide a histogram for the spectral density, which follows the microscopic RMT for several oscillations. With an increased statistics — which may be accessible using the optimized algorithmic techniques described in ref.  — we could evaluate the spectral correlators as well. Further values for $`\mathrm{\Sigma }(V)`$ would also allow us to extract an estimation for $`F_\pi `$ . Finally we are also about to extend this analysis to a hypercube overlap fermion, which is described in ref. . This is a different Ginsparg-Wilson fermion, and its study will provide a more complete picture of the spectral behavior of this class of lattice Dirac operators.

###### Acknowledgments.

We would like to thank Poul Damgaard for his help with the interpretation of the formulae in ref. . We are also indebted to him as well as Gernot Akemann, Michael Müller-Preußker, Jacques Verbaarschot and Tilo Wettig for useful discussions and many helpful comments. K.J. wants to thank Leonardo Giusti, Pilar Hernández, Martin Lüscher, Peter Weisz and Hartmut Wittig for helpful discussions about RMT and the interpretation of the eigenvalue distributions. We further thank Kei-ichi Nagai for communication regarding the numerical work. The computations have been performed at the Konrad Zuse Zentrum in Berlin and at the Forschungszentrum Jülich. This work was supported in part by the DFG Sonderforschungsbereich Transregio 9, “Computergestützte Theoretische Teilchenphysik”.


format: text

text:
# Comment on Conservative force fields in nonextensive kinetic theory

## Abstract

We discuss an improper application of the Maxwell-Boltzmann theory to a very unrealistic model. The authors of the original paper claim that a generalized theory would solve a problem that really does not exists. This Comment was submited to Physica A and was not accepted as it is.

and

In a recent paper Lima et al.  proposed to generalize the expression that describes how the density of an ideal gas varies with the height $`z`$, measured from the planet surface, in an isothermal planetary atmosphere. It is also assumed that the gravitational field is constant in the $`z`$-direction, so that the potential energy is simply $`U(z)=mgz`$. These simplified hypotheses lead to well-known barometric formula

$$\rho (z)=\rho _0\mathrm{exp}\left[\frac{mgz}{k_BT}\right],$$
(1)
where $`\rho _0`$ is the gas density at the planet surface level ($`z=0`$), $`m`$ is the gas mass, $`k_B`$ represents the Boltzmann constant and $`T`$ is the absolute temperature.

In the wake of Tsallis  non-extensive $`q`$-Thermostatistics Lima et al. found Eq. (24) presented in their paper . These authors claim that for suitable values of $`q`$ their expression would solve the problem of having an atmosphere that would extend to infinite. Actually, it is well known from elementary physics that there is really no such problem with respect to the Boltzmann distribution which leads to Eq. (1). Any book of physics for undergraduate students explains that there are two fundamental hypotheses behind Eq. (1), namely, that both $`g`$ and $`T`$ are constant as $`z`$ varies from $`0`$ to $`\mathrm{}`$. For typical heights $`z10`$ to $`100km`$ the error in $`g`$ is of the order of $`1\%`$, so there is not a serious problem to assume a constant $`g`$. On the other hand, it is empirically known that the temperature falls approximately $`5K`$ as we go $`1km`$ higher, if we start from the earth surface. Thus one notes that as $`z`$ goes from $`0`$ up to $`z=25km`$ the isothermal condition is violated, and we could not apply an equilibrium Boltzmann distribution to deduce Eq. (1). On this basis, the Table 1 presented in  is completely void of significance (we think there is problably a misprinting for the $`z_{max}`$ of Oxygen when $`q=0.8`$). Also, that Table gives figures for the Hydrogen and Oxygen gases. If the authors had considered the most abundant gases in our atmosphere, they would have found a somewhat better results considering $`q=0.7`$, namely, $`z_{max}=18km`$, and $`29km`$ for $`CO_2`$ and $`N_2`$, respectively. However, we were not able to find any criterion in their paper about the right choice for the value of the parameter $`q`$.

For pedagogical reasons we thing that any discussion relative to distances in this sort of problem should begin by considering the natural length scale for this problem, namely,

$$\xi =\frac{k_BT}{mg}.$$

Taking, for instance, $`T=300K`$ we can easily obtain $`\xi _{H_2}100km`$, whereas $`\xi 10km`$ for $`O_2`$, $`N_2`$ and $`CO_2`$.

In summary, contrary to the authors claim, there is really no problem with the Maxwell-Boltzmann distribution in the present context. The problem we are faced with in this case lies in the model which is completely inappropriated to discuss the question raised Lima et el. .


text:
# The origin of peculiar jet-torus structure in the Crab nebula

## 1 Introduction

The Crab Nebula is a prototype compact synchrotron nebula continuously powered by ultra-relativistic, magnetized wind from young rapidly rotating pulsar. As this nebula is confined within a nonrelativistic surrounding, the wind must terminate at a shock wave. It is the wind plasma heated to relativistic temperatures at this shock that fills the nebula and produces the observed nonthermal electromagnetic emission from the radio to the gamma-ray band (Rees & Gunn 1974, Kennel & Coroniti 1984). Indeed, early observations of the Crab Nebula revealed an apparent central hole in the nebula brightness distribution (Scargle 1969), which was identified with the termination shock in the theoretical models of the nebula. However, the assumption of spherical symmetry utilized in these models in order to simplify calculations was shattered when recent X-ray and optical observations (Hester et al. 1995, 2002; Weisskopf et al. 2000) revealed the inner structure of the nebula in much greater details. They discovered two jets aligned with the rotational axis of the pulsar and a luminous torus or disc in its equatorial plane. As similar structures have been found in other pulsar wind nebulae (Gaensler, Pivovaroff & Garmire 2001; Helfand, Gotthelf & Halpern 2001; Pavlov et al. 2001; Gaensler et al. 2002; Lu et al. 2002), we must be dealing with quite generic a phenomenon.

The interaction between a relativistic MHD wind from the pulsar and a dense, nonrelativistic surrounding is a much more challenging problem in the case of anisotropic wind and there is no much hope to find its analytical solution. Fortunately, recent progress in numerical methods for relativistic gas dynamics and MHD (Marti & Muller 1999; Komissarov 1999) has made possible to approach this problem numerically. The main results of the first attempt of such a study are briefly described in this letter. Full details will be presented elsewhere.

## 2 Structure of the pulsar wind in the far zone

Although a self-consistent solution to the problem of pulsar wind remains to be found, it is commonly accepted that far away from the pulsar the wind can be considered as almost radial supermagnetosonic outflow with purely azimuthal magnetic field and anisotropic angular distribution of the energy flux (Michel 1982; Beskin, Kuznetsova & Rafikov 1998; Chiueh, Li &Begelman 1998; Bogovalov & Tsinganos 1999). According to the simplified split-monopole models of pulsar magnetospheres (Michel 1973; Bogovalov 1999), the total energy flux density of the wind, $`f_w`$, satisfies the following equation

$$f_w=\frac{f_0}{r^2}(\mathrm{sin}^2\theta +1/\sigma _0),$$
(1)
where $`r`$ and $`\theta `$ are the spherical coordinates whose polar axis is aligned with the rotation axis of the pulsar. The first term in the brackets represents the Poynting flux, whereas the second one accounts for the small initial contribution of particles, $`\sigma _01`$. It was shown recently that the termination shock of such a wind is not at all spherical. In fact, it is located significantly closer to the pulsar along its rotational axis than in the equatorial plane (Lyubarsky 2002; Bogovalov & Khangoulyan 2002).

Crab’s jets, as well as jets of other pulsars, appear to originate from the pulsar (Weisskopf et al. 2000; Helfand, Gotthelf & Halpern 2001; Pavlov et al. 2001; 2003; Gaensler et al. 2002; Hester et al. 2002;). This seems to indicate that they are formed within the pulsar wind and the collimation by magnetic hoop stress suggests itself. However, a closer look reveals a number of problems with this explanation. First of all, such a collimation is found to be extremely ineffective in ultra-relativistic flows \[Beskin et al. 1998, Chiueh & Begelman 1998, Bogovalov & Tsinganos 1999, Lyubarsky & Eichler 2001\]. Moreover, the direct observations of proper motions in Crab’s and Vela’s jets indicate rather moderate velocities of only $`0.3÷0.7c`$ (Hester et al. 2002; Pavlov et al. 2003). In order to overcome these problems, Lyubarsky (2002) proposed that the jets are formed downstream of the termination shock where velocities are no longer ultra-relativistic and the magnetic collimation is much more effective. Because the termination shock is much closer to the pulsar along its rotational axis, this could give an impression of jets being produced by the pulsar itself.

Properties of the pulsar wind nebula strongly depend on the wind magnetization. It is widely accepted that a typical pulsar wind is launched as a Poynting-dominated flow but most of the electromagnetic energy is transferred to particle along the way to the termination shock (see, however, Begelman 1998). Although the problem of energy conversion remains a subject of intensive debate (e.g. Melatos 2002), a number of mechanisms have been proposed in recent years. Since the pulsar magnetic axis is inclined with respect to its rotational axis, a significant fraction of the Poynting flux is carried out by the component of electromagnetic field oscillating with the rotational period of the pulsar. These small-scale waves can decay via various dissipation processes (Lyubarsky & Kirk 2001; Melatos 2002; Lyubarsky 2003a; Kirk & Skjæraasen 2003). Even if the dissipation time-scale is larger than the time of travel from the pulsar to the termination shock, these waves rapidly decay at the termination shock. In this case, the post-shock parameters are still the same as if the alternating fields had already annihilated in the upstream flow (Lyubarsky 2003b). For this reasons, we assume in our model of the pulsar wind that all these waves have already decayed and transferred their energy to particles. The rest of the Poynting flux is transported by the large-scale azimuthal magnetic field. The exact latitudinal distribution of this field remains to be found but it must vanish both along the rotational axis, as any axisymmetrical magnetic field does, and in the equatorial plane, because the average magnetic field in the obliquely rotating magnetosphere is zero at the equator. We take the residual magnetic field in the wind in the form

$$B=\sqrt{\frac{4\pi f_0}{c}}\frac{\xi }{r}\mathrm{sin}\theta \left(1\frac{2\theta }{\pi }\right).$$
(2)
The free parameter $`\xi 1`$ controls the wind magnetization.

## 3 Numerical simulations

The initial solution includes 1) a spherically symmetric shell of cold dense gas expanding radially with velocity of 5000 km/s and 2) a radial ultra-relativistic wind with the energy flux density given by equation (1) with $`\sigma _0=100`$, and the magnetic field described by equation (2). Provided the wind is ultra-relativistic, the solution is not sensitive to the exact value of the wind Lorentz factor, which we set to $`\gamma _w=10`$. To carry out these axisymmetric simulations we used the Godunov-type scheme for relativistic magnetohydrodynamics constructed recently by Komissarov (1999).

The results for $`\xi =0.3`$ are presented in Figure 1 (overall structure of the flow) and Figures 2, 3 (its central part). One can see that instead of a single termination shock, the numerical solution displays a whole complex of shocks. The equatorial, weakly magnetized part of the pulsar wind terminates at almost cylindrical shock crossing the equator, “the Mach belt”. At higher latitudes, the flow passes first through a highly oblique ”arch shock”. Downstream of this shock the flow remains supermagnetosonic with the typical velocity of $`0.8÷0.9c`$. Then it passes through another shock, “the rim shock”, that originates from the edge of the Mach belt. In contrast to the spherically symmetric expansion assumed in the current theories of pulsar nebulae, most of the downstream flow is confined to the equatorial plane. Its typical velocity, $`v0.6c`$, agrees with the measurements of proper motions in Crab’s torus (Hester et al. 2002). This equatorial outflow is eventually pushed back by the pressure force and the magnetic hoop stress and forms a backflow.

This backflow reaches the central region of the nebula and forms large-scale vortices at intermediate latitudes.

Another backflow is launched from the surface of the equatorial disc as the result of strong magnetic braking. This inner backflow eventually forms magnetically collimated polar jets. The typical velocity of these jets is $`v_j0.5c`$, which is well in agreement with observed velocities of Crab’s and Vela’s jets (Hester et al. 2002; Pavlov et al. 2003). Similar jets are produced in the model with $`\xi =0.5`$ but not found in the solution with $`\xi =0.2`$. The pulsar wind magnetization is traditionally measured by the parameter $`\sigma `$, defined as the ratio of the Poynting flux to the kinetic energy flux of the wind (Kennel & Coroniti 1984). It turns out that in our model the mean $`\sigma 0.1\xi ^2`$ Thus, the polar jets are produced for $`\sigma 0.01`$, which is somewhat higher than the previous estimates for the Crab Nebula based on simplified theoretical models (Kennel & Coroniti 1984; Emmering & Chevalier 1987; Begelman & Li 1992).

Figure 4 shows the simulated synchrotron X-ray emission from the inner part of the nebula obtained for the same orientation of the wind relative to the observer as in the Crab nebula. In order to create this image, we assumed that electrons (and positrons) are injected at the termination shock and then suffer synchrotron energy losses at a rate determined by the typical value of magnetic field in the numerical solution. Because the magnetic field in our simulations is purely azimuthal, it vanishes on the symmetry axis and so does the synchrotron emissivity. In a real jet, the strong velocity shear would generate the poloidal component of the magnetic field, which does not have to vanish on the axis. To take this into account, we added, only inside the jet, the poloidal field aligned with the flow velocity at the level of 30% of the gas pressure.

## 4 Discussion

The obtained image is excitingly similar to the X-ray image of the Crab nebula (Weisskopf et al. 2000). In addition to the polar jets, one can see a system of rings which makes an impression of a disc-like or even a toroidal structure. Well in agreement with the observations these rings are brighter on the counter-jet side, which is entirely due to the relativistic beaming effect. The bright spot within the inner ring may correspond to the so-called “sprite” of the Crab nebula (Hester et al. 2002). This emission originates in the high velocity plasma flowing just above the arch shock where its velocity vector points directly toward the observer.

On the whole, our model captures the main properties of the Crab nebula remarkably well. There are, however, some qualitative and quantitative differences which demand further investigation. Crab’s jet is not so straight and well collimated and it cannot be traced that far away from the pulsar. It bends, spreads, and eventually merges into the surrounding plasma (Weisskopf et al. 2000). The most likely reason for such a behavior is the development of the kink instability (Begelman 1998), which is suppressed in our simulations by the condition of axisymmetry. Full 3-dimensional simulations are needed to overcome this restriction.

The brightness contrast between the jet side and the counter-jet side of the rings in the simulated maps is too high. Since this asymmetry is entirely due to the relativistic beaming, it strongly depends on the velocity field very close to the termination shock and hardly depends on anything else. The main factor determining this velocity field is the angular structure of the pulsar wind. Thus, the brightness asymmetry of Crab’s “torus” imposes strong observational constraints on the pulsar wind models. Future simulations will be used to determine the model parameters providing the best fit to the observational data.