text:
# Speed and Adaptability of Overlap Fermion Algorithms

## I Introduction

A lack of consistent definition of chiral fermions on the lattice has hampered definitive and convincing investigations of chiral aspects of quantum chromodynamics (QCD) until now. Thus important physics issues, such as the spontaneous breaking of the chiral symmetry at low temperatures and its restoration at finite temperature, have remained hostages to technical questions such as the fine-tuning of the bare quark mass (Wilson fermions) or the precise number of massless flavours (staggered fermions). Recent developments in defining exact chiral symmetry on the lattice have therefore created exciting prospects of studying an enormous amount of physics in a cleaner manner from first principles. However, the corresponding Dirac operators are much more complicated. Without good control of the algorithms needed to deal with them, one is unlikely to derive the full benefit of their better chiral properties. Our goal in this paper is to evaluate the efficiency of the most widely used, or most promising, algorithms. By efficiency we mean both the speed of the algorithm, which is measured by the computer time required to achieve a certain accuracy in the solution, and the adaptability, which is measured by how the speed scales as the problem becomes harder. This study is made for various values of the required accuracy along with the corresponding analysis on the accuracy obtained for the expected properties of the resulting Dirac operator such as the Ginsparg-Wilson relation, the central circle relation, $`\gamma _5`$ hermiticity or normality. In particular, we have observed that these properties can be satisfied accurately only if the sign computation of the Wilson Dirac operator has high enough precision.

### I.1 The overlap Dirac operator

One version of chiral fermions on the lattice is the overlap formalism. The overlap Dirac operator ($`D`$) is defined neu2  in terms of the Wilson-Dirac operator ($`D_w`$) by the relation

$$D=1+D_w(D_w^{}D_w)^{1/2}.$$
(1)
In this paper we shall use the shorthand notation

$$M=D_w^{}D_w.$$
(2)
The Wilson-Dirac operator $`D_w`$ (for lattice spacing $`a=1`$) is given by

$$D_w=\frac{1}{2}\underset{\mu }{}\left[\gamma _\mu \left(_\mu +_\mu ^{}\right)_\mu _\mu ^{}\right]+m,$$
(3)
where $`_\mu `$ and $`_\mu ^{}`$ are (gauge covariant) forward and backward difference operators respectively. It has been shown that as long as the mass $`m`$ is in the range $`2<m<0`$, the above overlap Dirac operator is well-defined, and corresponds to a single massless fermion. Furthermore, it satisfies the Ginsparg-Wilson relation gw

$$\gamma _5D+D\gamma _5=D\gamma _5D,$$
(4)
which leads to a good definition of chirality on the lattice and has been shown to correspond to an exact chiral symmetry on the lattice.

The overlap Dirac operator, $`D`$, enjoys many nice properties in addition to the Ginsparg-Wilson relation in Eq. (4). In particular, it satisfies $`\gamma _5`$-Hermiticity—

$$D^{}=\gamma _5D\gamma _5.$$
(5)
Together with the Ginsparg-Wilson relation, this implies that $`D`$ is normal, i.e.,

$$[D,D^{}]=0.$$
(6)
Normality clearly means that $`D`$ and $`D^{}`$ have the same eigenvectors. Eqs. (4,5) also imply

$$D+D^{}=D^{}D,$$
(7)
and hence the eigenvalues of $`D`$ and $`D^{}`$ lie on the unit circle centered at unity on the real line, implying that $`D1`$ is unitary. We define measures of numerical errors on each of these quantities, and relations between them in Section II.

### I.2 Numerical algorithms

All computations of hadronic correlators involve the determination of the fermion propagator $`D^1`$, and need a nested series of two matrix iterations for their evaluation, since each step in the numerical inversion of $`D`$ involves the evaluation of $`M^{1/2}`$. This squaring of effort makes a study of QCD with overlap quarks very expensive.

This problem defines for us the properties that an efficient algorithm to deal with $`M^{1/2}`$ must have. First, it should achieve a given desired accuracy as quickly as possible. The need for accuracy is clear: the accuracy to which the Ginsparg-Wilson relation in Eq. (4) is satisfied depends on the accuracy achieved in the computation of $`M^{1/2}`$. The second, and equally important, requirement is that the method should adapt itself easily to matrices with widely different condition numbers—

$$\kappa =\frac{\lambda _{\mathrm{max}}}{\lambda _{\mathrm{min}}},$$
(8)
where $`\lambda _{\mathrm{min}}`$ and $`\lambda _{\mathrm{max}}`$ are, respectively, the minimum and the maximum eigenvalues of $`M`$. Adaptability is needed because in QCD applications the eigenvalue spectrum of $`M`$ can fluctuate over many orders of magnitude from one configuration to the next. Since the condition number on a configuration is a priori unknown, a method with low adaptability will end up either being inefficient or inaccurate or even both. Algorithms, which are adaptable in principle, may require tuning of parameters by hand, or they may incorporate a procedure for self-tuning. Clearly, self-tuning algorithms are the ones that can best deal with fluctuating condition numbers in any real situation.

In this paper we examine the speed and adaptability of several different algorithms to compute the inverse square root, $`M^{1/2}`$, of Hermitean matrices (in our applications the eigenvalues of $`M`$ are non-negative ) acting on a vector. Several algorithms for this have been proposed in the literature.

We do not consider the first algorithm to be proposed, since this requires a matrix inversion to be performed at each step of an iteration to determine $`M^{1/2}`$ chiu . Later algorithms are more efficient. These fall into two classes— expansions of $`1/\sqrt{z}`$ in appropriate classes of functions (rational functions neu  or orthogonal polynomials hepilu ), and iterative methods borici . We have analyzed four derived algorithms, namely the Optimized Rational Approximation (ORA) nara , the Zolotarev Approximation (ZA, which is also a rational expansion) wupp ; ch , the Chebychev Approximation (CA, a polynomial expansion) hepile  and the Conjugate Gradient Approximation (CGA, an iterative method) ragula . We find that an expansion in Chebychev polynomials is the fastest when $`\kappa `$ is moderately large, but it suffers from various instabilities including a lack of adaptability. Rational expansions cure many of the instabilities of polynomial expansions; indeed the ZA is the fastest and is adaptable but not self-tuning. An iterative method is the only fully self-tuned algorithm, and it turns out to be reasonable also from the point of view of speed.

### I.3 Algorithmic costs and adaptability

We make two different estimates of the cost of each algorithm. The complexity, $`𝒞`$, counts the number of arithmetic operations required to achieve the solution to the problem and is a measure of speed independent of the specific machine on which the algorithm is implemented. The spatial complexity, $`𝒮`$, is the memory requirement for the problem. While timing runs on particular machines on chosen test configurations are instructive, the scaling of speed for each algorithm with physical and algorithmic parameters is provided by our estimates of $`𝒞`$.

Our estimate of the adaptability, $`𝒜`$, of each algorithm is the following. If the scalar algorithm for $`1/\sqrt{z}`$ is tuned to have maximum relative error $`\epsilon `$ in the range $`[z_{\mathrm{min}},z_{\mathrm{max}}]`$, then we find the smallest range $`[z_{\mathrm{min}}^{},z_{\mathrm{max}}^{}]`$ where the error is at most $`10\epsilon `$. Note that the second interval cannot be smaller than the first. In terms of these quantities we define

$$𝒜=\frac{\mathrm{log}(z_{\mathrm{max}}^{}/z_{\mathrm{min}}^{})}{\mathrm{log}(z_{\mathrm{max}}/z_{\mathrm{min}})}1.$$
(9)
The least adaptable algorithms have small values of $`𝒜`$ ( $`𝒜>0`$ ). $`𝒜`$ is a measure of the relative accuracy achieved in a fixed CPU time for the same algorithm running on two different configurations with condition numbers $`\kappa =z_{\mathrm{max}}/z_{\mathrm{min}}`$ and $`\kappa ^{}=z_{\mathrm{max}}^{}/z_{\mathrm{min}}^{}`$. In conjunction with estimates of $`𝒞`$, it also contains information about the scaling of CPU time required to achieve the same accuracy on the two configurations.

### I.4 Numerical tests

Our numerical tests were performed with three typical $`SU(3)`$ gauge configurations on a $`4\times 12^3`$ lattice at $`\beta =5.80`$ (i.e., $`T=1.25T_c`$). The configuration A had eigenvalues of $`M`$ in the range $`[0.032,32]`$ so that $`\kappa =10^3`$. The configuration B had eigenvalues in the range $`[7.2\times 10^5,32]`$, giving $`\kappa =4.4\times 10^5`$. Finally, configuration C had eigenvalues in the range $`[8.9\times 10^9,32]`$ and hence $`\kappa =3.6\times 10^9`$. Configuration A is one of the easiest configuration we found in our simulations, and there were several configurations with $`\kappa `$ of order $`10^5`$$`10^9`$ in the data set we worked with in ragula . If there is a single scale in the level spacing of the eigenvalues of $`M`$, then we expect $`\kappa =𝒪(V)7\times 10^3`$ on our test configurations. We conclude that A is indeed a little easier than the generic configuration and B and C are successively harder. The CPU times we quote in our tables are obtained on a Fujitsu VPP5000, which is a vector computer. Our computations ran on this machine at a speed of around 4.1 Gigaflops.

In sections IIIVI we describe and analyze the four algorithms and also present data on precision and time measurements on these three test configurations. In this work we have not investigated the performance of the algorithms with deflation (explicit subtraction) of some eigenvectors. However, we do remark on the precision required for deflation and the propagation of errors due to such a subtraction. Section VII contains a comparison of the numerical results and our conclusions.

## II Errors

In general, every numerical method to compute $`M^{1/2}`$ constructs an operator $`L`$ which applied to a vector $`\mathrm{\Phi }`$ gives the vector

$$X=L[\mathrm{\Phi }]\mathrm{with}L[\mathrm{\Phi }]=M^{1/2}\mathrm{\Phi }+E[\mathrm{\Phi }],$$
(10)
where $`E`$ is the error in the approximation, $`L`$, to the matrix $`M^{1/2}`$. Typically, these operators $`L`$ and $`E`$ are not matrices because the algorithms can introduce a dependence on $`\mathrm{\Phi }`$ which is not linear. The error $`E[\mathrm{\Phi }]`$ on the computation of $`M^{1/2}`$ leads to the violation of the properties in Eqs. (4-7). In our numerical tests we have investigated five measures of the accuracy of the algorithms through norms of the following operators—

$`Z_{1/2}=ML^21,`$ $`Z_{GW}=D\gamma _5+\gamma _5DD\gamma _5D,`$ $`Z_N=DD^{}D^{}D,`$
$`Z_H=D^{}\gamma _5D\gamma _5,`$ $`Z_{CC}=D+D^{}D^{}D,`$ (11)
where $`D=1+D_wL`$ and $`D^{}=1+LD_w^{}`$. Each of these operators is zero when $`E[\mathrm{\Phi }]`$ = 0. With Gaussian random vectors $`\mathrm{\Phi }`$, we have measured the deviations from this exact value through

$$ϵ_i=|Z_i\mathrm{\Phi }|/|\mathrm{\Phi }|\mathrm{and}ϵ_i^{}=\mathrm{\Phi }^{}Z_i\mathrm{\Phi }/|\mathrm{\Phi }|^2.$$
(12)
Here, and later, we use the notation $`|v|=\sqrt{v^{}v}`$ for any complex vector $`v`$. Note that $`ϵ_i`$ are real and non-negative whereas $`ϵ_i^{}`$ are complex in general.

### II.1 Polynomial approximation

The polynomial approximation consists of writing

$$L=\underset{i=1}{\overset{N_O}{}}b_iM^i,$$
(13)
where $`b_i`$ are constants. It is clear that both $`L`$ and $`E`$ are matrices in this case and

$$[L,M]=0.$$
(14)
Since in numerical implementations of $`D_w`$, $`\gamma _5`$-Hermiticity is accurate to machine precision, i.e., $`|(D_w^{}\gamma _5D_w\gamma _5)\mathrm{\Phi }|=0`$, one has the following relation:

$$\gamma _5D_wM^n\gamma _5=M^nD_w^{}.$$
(15)
Its direct consequence is

$$Z_H=LD_w^{}\gamma _5D_wL\gamma _5=0.$$
(16)
Using Eqs.(14-15), one can easily obtain the following relations between the various $`Z`$’s :

$`Z_{CC}`$ $`=`$ $`Z_{1/2},`$
$`Z_{GW}`$ $`=`$ $`\gamma _5Z_{CC},`$
$`Z_N`$ $`=`$ $`Z_{CC}\gamma _5Z_{CC}\gamma _5.`$ (17)
As a consequence,

$`|ϵ_{1/2}^{}|`$ $`=`$ $`|ϵ_{CC}^{}|,`$
$`ϵ_{1/2}`$ $`=`$ $`ϵ_{CC}=ϵ_{GW},`$
$`ϵ_H^{}`$ $`=`$ $`ϵ_H=0.`$ (18)
Expanding $`Z_{1/2}=ML^21=Z_{1/2}^{}`$ in powers of $`E`$ and retaining only the leading order terms, we find that $`Z_{1/2}^{}Z_{1/2}=4E^2M`$ and $`Z_{1/2}=2EM^{1/2}`$. Defining the averages of $`ϵ_i`$ and $`ϵ_i^{}`$ over an ensemble of $`\mathrm{\Phi }`$ as $`\overline{ϵ_i^2}=\mathrm{Tr}Z_i^{}Z_i`$ and $`\overline{ϵ_i^{}}=\mathrm{Tr}Z_i`$, one obtains,

$`\overline{ϵ_{1/2}^2}=\overline{ϵ_{CC}^2}=\overline{ϵ_{GW}^2}=4{\displaystyle 𝑑\lambda \rho (\lambda )\lambda ϵ^2(\lambda )},`$
$`\overline{ϵ_{1/2}^{}}=\overline{ϵ_{CC}^{}}=2{\displaystyle 𝑑\lambda \rho (\lambda )\sqrt{\lambda }ϵ(\lambda )},`$
$`\overline{ϵ_{GW}^{}}=2{\displaystyle 𝑑\lambda \mathrm{\Delta }\rho (\lambda )\sqrt{\lambda }ϵ(\lambda )},`$ (19)
where $`\rho (\lambda )`$ is the density of eigenvalues of $`M`$, $`ϵ(\lambda )`$ is the error in the approximation and $`\mathrm{\Delta }\rho (\lambda )=\rho _+(\lambda )\rho _{}(\lambda )`$, the difference between the spectral densities (of $`M`$) in the chiral positive and negative sectors. Note that $`\sqrt{\lambda }ϵ(\lambda )`$ is the relative error in the determination of the inverse square root, and all the integrals depend only on this relative error. Since there is no reason for $`\rho (0)`$ to vanish, it is clear that the error in the expansion must remain under control even as $`\lambda 0`$. Clearly, this is impossible to arrange in polynomial expansions for $`1/\sqrt{\lambda }`$. However, a finite sample of gauge configurations does not need full control over $`ϵ(0)`$, but only for $`ϵ(\lambda _<)`$, where $`\lambda _<`$ is the smallest eigenvalue encountered in the sample. To achieve this while optimizing CPU costs on configurations where all the eigenvalues are much larger requires the algorithm to be adaptive.

It was assumed above that no deflation has been performed, or that deflation has been performed with no arithmetic errors. We comment on the effects of deflation in Section VII.

### II.2 Rational approximation

In case of a rational function approximation to $`M^{1/2}`$, one writes the operator $`L`$ as

$$L=\underset{i=1}{\overset{N_O}{}}\frac{b_i}{M+d_i}$$
(20)
$`E`$ here depends on the order $`N_O`$ and the accuracy of the inversion of $`(M+d_i)`$. If the inversion can be achieved with infinite precision, then $`L`$ is a matrix again which commutes with $`M`$ and the analysis of the previous subsection applies in full. If, on the other hand, the error due to the inversion dominates, then for many algorithms, such as the Conjugate Gradient, $`L`$ depends explicitly on the vector $`\mathrm{\Phi }`$ in a complicated way and it is not a matrix. One has to compute the different errors explicitly and study their behavior as in iterative methods. Thus the behavior of errors from rational approximation case interpolates between that of the polynomial approximation and an iterative method according to the relation between the order and the precision of the inversion.

## III Fixed order: Chebychev approximation

The first use of the polynomial approximation utilized Legendre polynomials hepilu . Later the same group proposed a more robust version using the Chebychev approximation (CA) hepile . As is well known, when expanding any function, $`f(z)`$ in a fixed range $`z_{\mathrm{min}}zz_{\mathrm{max}}`$, to a given order $`N_O`$ through orthogonal polynomials, the use of Chebychev polynomials minimizes the maximum error on the function to be approximated.

For the function $`1/\sqrt{z}`$, the coefficients of the Chebychev expansion for $`z_{\mathrm{min}}zz_{\mathrm{max}}`$ to order $`N_O`$ are

$$C_k=\frac{2\sqrt{2}}{N_O}\underset{j=1}{\overset{N_O}{}}\frac{\mathrm{cos}((k1)(j\frac{1}{2})\pi /N_O)}{[z_{\mathrm{min}}+z_{\mathrm{max}}+(z_{\mathrm{min}}z_{\mathrm{max}})\mathrm{cos}((j\frac{1}{2})\pi /N_O)]^{1/2}}.$$
(21)
Applying this approximation to a matrix $`M`$ corresponds to finding $`L\mathrm{\Phi }`$ (see Eq. 10) by the expansion—

$$L\mathrm{\Phi }=\frac{1}{2}c_1\mathrm{\Phi }^{(0)}+\underset{k=2}{\overset{N_O}{}}c_k\mathrm{\Phi }^{(k1)},$$
(22)
where the successive vectors $`\mathrm{\Phi }^{(n)}`$ can be found by the iteration $`\mathrm{\Phi }^{(0)}=\mathrm{\Phi }`$ and

$$\mathrm{\Phi }^{(n)}=\frac{2}{z_{\mathrm{max}}z_{\mathrm{min}}}\left[2M\mathrm{\Phi }^{(n1)}(z_{\mathrm{max}}+z_{\mathrm{min}})\mathrm{\Phi }^{(n1)}\right]\mathrm{\Phi }^{(n2)},$$
(23)
for $`n2`$. For $`n=1`$ the term $`\mathrm{\Phi }^{(1)}`$ is dropped from the recursion.

There are various sources of error, which we now analyze in succession. For the scalar version of the algorithm (or for each eigenvector of $`M`$ separately), with fixed $`z_{\mathrm{min}}`$, $`z_{\mathrm{max}}`$ and parameter $`N_O`$, there is an error in the approximation of $`1/\sqrt{z}`$ which we call $`ϵ(z)`$. For computations at arbitrary precision with a fixed range, $`[z_{\mathrm{min}},z_{\mathrm{max}}]`$, $`ϵ(z)`$ depends entirely on $`N_O`$. To keep the absolute relative error bounded, $`|ϵ(z)|\sqrt{z}\epsilon `$, as $`z_{\mathrm{min}}`$ changes, we must tune

$$N_O\mathrm{log}\left(\frac{1}{\epsilon }\right)\sqrt{\frac{1}{z_{\mathrm{min}}}}\mathrm{log}\left(\frac{1}{\epsilon }\right)\sqrt{\kappa },$$
(24)
as shown in Figure 1. The last expression follows if we choose $`z_{\mathrm{min}}=\lambda _{\mathrm{min}}`$ and $`z_{\mathrm{max}}=\lambda _{\mathrm{max}}`$.

However, as shown in the right panel of the figure, tuning $`N_O`$ in this manner causes the relative error outside the range to blow up. The CA has no tolerance to violations of the requirement on the range. This is easy to understand. In any polynomial approximation, with decreasing $`z_{\mathrm{min}}`$ larger values of $`N_O`$ are required for keeping the error fixed within the interval $`[z_{\mathrm{min}},z_{\mathrm{max}}]`$. However, outside this interval, the error then increases as

$$ϵ(z)\left(\frac{2zz_{\mathrm{min}}z_{\mathrm{max}}}{z_{\mathrm{max}}z_{\mathrm{min}}}\right)^{N_O},\mathrm{for}z<z_{\mathrm{min}}\mathrm{or}z>z_{\mathrm{max}},$$
(25)
and hence the error increases faster as $`z_{\mathrm{min}}`$ decreases. Solving this for $`z`$, given some fixed value of $`ϵ(z)/\epsilon `$, we can easily see that for large $`\kappa `$ and fixed precision $`\epsilon `$,

$$𝒜\mathrm{exp}(\alpha \sqrt{\kappa })$$
(26)
where $`\alpha `$ is some number.

The effect of finite arithmetic precision can also be analyzed easily since the iteration in Eq. (23) is linear. Any arithmetic error, $`\delta ^{(m)}`$, in $`\mathrm{\Phi }^{(m)}`$ of the order of the machine precision remains in control whenever all eigenvalues, $`\lambda `$ of $`M`$ satisfy $`z_{\mathrm{min}}\lambda z_{\mathrm{max}}`$. If any eigenvalue of $`M`$ lies outside this range, then the iteration magnifies the error geometrically—

$$\delta ^{(m+n)}\left(\frac{2\lambda _{\mathrm{min}}z_{\mathrm{min}}z_{\mathrm{max}}}{z_{\mathrm{max}}z_{\mathrm{min}}}\right)^n\delta ^{(m)},$$
(27)
This is the vector version of the low adaptability of this algorithm. If estimates of $`\lambda _{\mathrm{min}}`$ and $`\lambda _{\mathrm{max}}`$ for $`M`$ are available, then, in view of this instability, it is best to choose $`z_{\mathrm{min}}<\lambda _{\mathrm{min}}`$ and $`z_{\mathrm{max}}>\lambda _{\mathrm{max}}`$.

The complexity of this algorithm is clearly dominated by the time required to operate upon a vector by the matrix $`M`$ in the iterations in Eq. (23). For the Wilson-Dirac matrix this time is of order $`V`$. Neglecting the time taken for scalar operations in the remainder of the algorithm, and also the order $`V`$ time for vector additions, in comparison with this, we have—

$$𝒞_{CA}wN_OVw^{}V\mathrm{log}\left(\frac{1}{\epsilon }\right)\sqrt{\kappa },$$
(28)
where $`wV`$ is the complexity of operating upon a vector by $`M`$ and $`w,w^{}`$ are constants. Apart from the gauge configuration, in QCD applications the storage required is for the three vectors needed for the iteration in Eq. (23). The space complexity is therefore

$$𝒮_{CA}=8N_c(N_c+3)V,$$
(29)
(for $`N_c`$ colors) neglecting storage for scalars.

With a fixed $`N_O`$, the precision of the algorithm deteriorates sharply when one or more of the eigenvalues of $`M`$ lie outside the interval $`[z_{\mathrm{min}},z_{\mathrm{max}}]`$, as shown in Figure 1 and by the low value of $`𝒜`$ in Eq. (26). This is often sought to be corrected for by deflating, i.e., explicitly dealing with the eigenspace of the lowest eigenmodes, and applying the algorithm to the orthogonal space. This would keep the accuracy constant as the condition number changes. If $`N_D`$ vectors need to be deflated, then the contribution to $`𝒞`$ clearly increases as $`(N_DV)^2`$, since each vector has to be orthogonalized with respect to every other. Also, $`𝒮`$ increases as $`N_DV`$ due to the necessity of storing the vectors. In order to achieve a target precision, $`\epsilon `$, on all gauge configurations, we are forced to deflate all vectors with $`\lambda <z_{\mathrm{min}}`$. Working with a fixed $`N_D`$ forces us to do unnecessarily large amount of work on most configurations, while still failing on a small set of configurations. As a result, $`N_D`$ has to be chosen appropriately for each configuration. With deflation then we have

$`𝒞_{CA}`$ $`=`$ $`w^{}V\mathrm{log}\left({\displaystyle \frac{1}{\epsilon }}\right)\sqrt{\kappa _{eff}}+w^{\prime \prime }V^2N_{D}^{}{}_{}{}^{2},`$
$`𝒮_{CA}`$ $`=`$ $`8N_c(N_c+3)V+8N_cN_DV,`$ (30)
where the angular brackets denote averages over gauge configurations, $`w^{\prime \prime }`$ is a constant independent of $`V`$ and $`N_D`$, and $`\kappa _{eff}`$ is the condition number of the matrix after deflating $`N_D`$ vectors. With careful programming we can arrange to make $`w^{\prime \prime }<w^{}`$, although they cannot differ by orders of magnitude ($`w^{\prime \prime }`$ can depend on $`\epsilon `$ and $`\lambda _{\mathrm{min}}`$).

We have not investigated the expectation values of $`N_D`$. However, when we change the volume at fixed physics, we expect that $`N_DV\overline{\rho }`$, where $`\overline{\rho }`$ is the average density of eigenvalues of $`M`$ near $`\lambda _{\mathrm{min}}`$. Since deflation is designed to hold $`\kappa _{eff}`$ fixed, this means that for large enough $`V`$ the complexity $`𝒞_{CA}V^4`$. On the other hand, $`\kappa `$ should generically grow linearly in $`V`$. Hence, on sufficiently large volumes, without deflation we would have $`𝒞_{CA}V^{3/2}`$. For best actual performance, one would have to tune the playoff between these two limits.

A further complication arises in CA, and indeed, in any method which utilizes a polynomial expansion. For any fixed finite precision, the deflation of $`\mathrm{\Phi }^{(0)}`$ is inaccurate since each component of the deflated vector is in error at least in the least significant bit. Since the CA iteration of Eq. (23) is not stable on the deflated eigenspace, this error blows up geometrically as in Eq. (27). Consequently, more and more bits are corrupted, as the iteration proceeds, and the process may eventually render the whole computation unusable. Note that this problem becomes more acute with decreasing $`\kappa `$, even if several eigenvectors are deflated. To prevent the error from swamping the result in this fashion, one has to reorthogonalize $`\mathrm{\Phi }^{(n)}`$ repeatedly (this process itself is not free of complications, see golub ). This involves finding $`N_D`$ dot products and subtracting $`N_D`$ vectors. While this is crucial in maintaining the accuracy of the result, it does not change the complexity, and we still have the results in Eq. (30).

A different approach, and the one we have adopted for our tests, is to start the algorithm by making an estimate of $`\lambda _{\mathrm{min}}`$ and $`\lambda _{\mathrm{max}}`$ and then to select $`N_O`$ accordingly. This obviates any need for deflation, and controlling the rounding errors in such a method. The algorithm is well-behaved, both with respect to precision and propagation of rounding errors, whenever $`z_{\mathrm{min}}\lambda _{\mathrm{min}}\lambda _{\mathrm{max}}z_{\mathrm{max}}`$. However, in this case the complexity rises as $`\sqrt{\kappa }`$.

As shown in Table 1, the algorithm performs well on configuration A. $`N_O`$ needed to achieve a given value of $`ϵ_{1/2}`$ is seen to rise logarithmically, as argued above. Consequently, $`ϵ_{1/2}`$ can be made very small and the required chiral properties obtained at any desired precision. Note also that the equalities (18) are exactly satisfied, as expected, and remain so for configurations B and C too, as seen from the Tables 2 and 3 respectively. However, the algorithm was found to be extremely slow for these configurations and saturated at $`N_O1500`$ with $`ϵ_{1/2}=4\times 10^3`$ on the configuration C, leading to the same precision for both the GW relation and the unit circle property of $`D`$.

## IV Fixed order: Optimized Rational Approximation

The first algorithm to compute $`M^{1/2}`$ through a rational expansion was a polar formula introduced by Neuberger neu . An improved version nara , called the optimized rational approximation, uses coefficients obtained numerically to give the approximation

$$L\mathrm{\Phi }=\sqrt{\alpha }\left[c_0+\underset{k=1}{\overset{N_O}{}}\frac{c_k}{\alpha M+d_k}\right]\mathrm{\Phi },$$
(31)
where $`\alpha `$ is an arbitrary scale whose choice we describe later, and the values of $`c_k`$ and $`d_k`$ for a given $`N_O`$ are obtained by optimizing the fit on some fixed interval $`[z_{\mathrm{min}},z_{\mathrm{max}}]`$ using the Remez algorithm (see nara  for details). The inversion of $`(\alpha M+d_k)`$ is made by a multimass conjugate gradient stopped according to a value $`ϵ_{CG}`$ for the residual. This version is used in particular in nara2 , liu , and degrand .

Taking the values of $`c_k`$ and $`d_k`$ for $`N_O=14`$, $`z_{\mathrm{min}}=0.01`$ and $`z_{\mathrm{max}}=1`$ from nara , we show in Figure 2 the relative accuracy, $`|ϵ(z)|\sqrt{z}`$ for the expansion both inside and outside the fitted range. It is clear from the figure that the algorithm rapidly degrades outside the chosen range. A numerical computation shows that

$$𝒜1,$$
(32)
so that it has much higher adaptability than the CA. Nevertheless, it makes large errors on the eigenvalues of $`M`$ which are greater than $`z_{\mathrm{max}}`$. It is easy to see that many eigenvalues are greater than unity, and a scaling factor $`\alpha `$ is therefore needed to bring these into range. The tuning of $`\alpha `$ is shown in Table 4. For the conjugate gradient inversion of each term, when $`ϵ_{CG}`$ is large, it determines the accuracy of the solution. Hence $`ϵ_{CG}`$ must be kept small enough so that the accuracy is $`ϵ`$.

In order to specify the scaling of the CPU time in each algorithm with various parameters of the problem, we count the complexity of the method. For this algorithm, $`𝒞`$ is proportional to the number of steps of the Conjugate Gradient inversion, $`N_{CG}`$. It can be proved that $`N_{CG}`$ grows no faster than $`\sqrt{\kappa }`$. However, from the observed convergence rate of the CG iterations (shown in Section VII), we see that $`N_{CG}\mathrm{log}(1/ϵ_{CG})\mathrm{log}\kappa `$. This expression can be used when $`ϵ_{CG}`$ is tuned to be smaller than the error shown in Figure 2. For each step of the multimass CG inversion, the complexity is dominated by the time required to operate upon a vector by the matrix $`M`$. For the Wilson-Dirac matrix this is of order $`V`$. Neglecting the time taken for scalar operations and also the order $`V`$ time for vector additions, we have—

$$𝒞_{ORA}wN_{CG}Vw^{}V\mathrm{log}\left(\frac{1}{ϵ_{CG}}\right)\mathrm{log}\kappa ,$$
(33)
where $`wV`$ is the complexity of operating upon a vector by $`M`$, and $`w^{}`$ is a constant. The dependence of $`N_{CG}`$ on $`V`$ is very weak for realistic $`N_{CG}`$ and is neglected. The memory requirement is essentially for the storage of the gauge configuration and for $`2+2N_O`$ vectors required to build up the approximation. The space complexity is therefore

$$𝒮_{ORA}=8N_c(N_c+2+2N_O)V,$$
(34)
(for $`N_c`$ colors) neglecting storage for scalars.

With a fixed value of $`N_O`$, the precision of the algorithm deteriorates when the eigenvalues of $`M`$ lie outside the range $`[z_{\mathrm{min}},z_{\mathrm{max}}]`$, as shown in Figure 2. For $`D_w`$, the highest eigenvalue remains in the vicinity of 32 for most configurations, whereas the lowest eigenvalue may fluctuate by several orders of magnitude. The large eigenvalues are brought into range by tuning $`\alpha <1`$ as shown already. However, this drives the lowest eigenvalues further out of the range, thus degrading performance. As a result, it is necessary to deflate, i.e., explicitly deal with the eigenspace of the lowest eigenmodes, and apply ORA on the orthogonal space in order to keep a constant accuracy as the condition number changes. As before, this changes the complexity to

$`𝒞_{ORA}`$ $`=`$ $`w^{}V\mathrm{log}\left({\displaystyle \frac{1}{ϵ_{CG}}}\right)\mathrm{log}\kappa _{eff}+w^{\prime \prime }N_{D}^{}{}_{}{}^{2}V^2,`$
$`𝒮_{ORA}`$ $`=`$ $`8N_c(N_c+2+2N_O)V+8N_cN_DV.`$ (35)
where the angular brackets denote averaging over the sample of configurations used. From available data on the growth of $`N_O`$ required to keep the relative error fixed with growing $`\kappa `$ nara  it seems that it is better to increase $`N_O`$ rather than to keep $`\kappa _{eff}`$ fixed by increasing $`N_{D}^{}{}_{}{}^{2}`$ with increasing volume at fixed physics.

In Tables 5, 6 and 7 we show the results of our numerical tests of the ORA using the set of $`c_k`$ and $`d_k`$ for $`N_O=14`$ from nara  for the three configurations already described. One sees that with higher precision of inversion, $`ϵ_{CG}`$, the relations (18) indeed get satisfied well. The gradual deterioration of the performance of ORA with fixed $`N_O`$ in going to larger $`\kappa `$ is clear from the tables. This indicates that the performance of ORA may improve if a degree of adaptability can be built into the algorithm by, for example, allowing for changes in $`N_O`$ because with $`N_O=14`$, we cannot obtain better precision than 0.00001, 0.00002 and 0.003 for the GW relation for configurations A, B and C respectively.

## V Fixed order: Zolotarev approximation

The Zolotarev algorithm was explored in wupp  and used in ch . It is a rational expansion defined by

$$L\mathrm{\Phi }=\underset{l=1}{\overset{N_O}{}}\left(\frac{b_l}{M+d_l}\right)\mathrm{\Phi },$$
(36)
in the range $`[z_{\mathrm{min}},z_{\mathrm{max}}]`$ (the smallest and largest eigenvalues of $`M`$ must satisfy the conditions $`z_{\mathrm{min}}\lambda _{\mathrm{min}}\lambda _{\mathrm{max}}z_{\mathrm{max}}`$), and the expansion coefficients are

$$d_l=c_{2l1}\mathrm{and}b_l=d_0\frac{_{i=1}^{N_O1}(c_{2i}c_{2l1})}{_{i=1,il}^{N_O1}(c_{2i1}c_{2l1})}.$$
(37)
The $`c_l`$’s are

$$c_l=\frac{\mathrm{sn}^2(lK/2N_O;\sqrt{1z_{\mathrm{min}}/z_{\mathrm{max}}})}{\mathrm{cn}^2(lK/2N_O;\sqrt{1z_{\mathrm{min}}/z_{\mathrm{max}}})},$$
(38)
where the values of the Jacobi elliptic functions, $`\mathrm{sn}(u,k)=\mathrm{sin}\eta `$ and $`\mathrm{cn}(u,k)=\mathrm{cos}\eta `$, are defined by the integral

$$u(\mathrm{sin}(\eta ))=_0^{\mathrm{sin}(\eta )}\frac{dt}{\sqrt{(1t^2)(1k^2t^2)}}.$$
(39)
The constant in Eq. (38), $`K=u(1)`$, is the complete elliptic integral. When sn is near 0 or 1, high precision in the expressions of the coefficients of the corresponding $`c`$’s is essential. The constant $`d_0`$ in Eq. (37) can be expressed in term of elliptic theta function chiu2 , or equivalently, fixed by the condition wupp

$$\underset{z}{\mathrm{min}}\underset{i=1}{\overset{N_O}{}}\left(\frac{\sqrt{z}b_l}{z+d_l}\right)+\underset{z}{\mathrm{max}}\underset{i=1}{\overset{N_O}{}}\left(\frac{\sqrt{z}b_l}{z+d_l}\right)=2.$$
(40)
As in the ORA, the multimass CG which is used to invert the terms in Eq. (36) should have a stopping criterion, $`ϵ_{CG}`$. One advantage of ZA over ORA is that the quantities $`d_l`$ in Eq. (36) are larger than those in Eq. (31). As a result, a multimass CG inverter can evaluate this approximation somewhat faster. Another advantage, emphasized in wupp  is that $`N_O`$ required for a certain accuracy is smaller for ZA than for ORA. It was found that a relative accuracy of better than 1 part in $`10^5`$ is obtained for the interval $`[0.01,1]`$ with $`N_O6`$ in ZA, as compared to 14 in ORA. As shown in Figure 3, the relative error for ZA in the range $`[10^4,32]`$ does not require significantly higher $`N_O`$ for similar control over error. However, the low adaptability, $`𝒜0.01`$, means that the coefficients should be computed over a range appropriate to the condition number of the matrix. The main effect of increasing the range for a fixed $`N_O`$ is to change the coefficients $`b_l`$ and $`d_l`$ in such a way that the logarithmic range of $`d_l`$ increases. We found that a factor 100 decrease in $`z_{\mathrm{min}}`$ (for fixed $`z_{\mathrm{max}}=32`$) led to a factor 20–30 decrease in the ratio of the minimum and maximum values of $`d_l`$.

The complexity and spatial complexity of ZA are very similar to that of ORA. The complexity is dominated by the matrix-vector multiplication in the CG inversions, and the memory requirement is dominated by the vectors in the multimass CG. Hence

$`𝒞_{ZA}`$ $``$ $`w^{}V\mathrm{log}\left({\displaystyle \frac{1}{ϵ_{CG}}}\right)\mathrm{log}\kappa ,`$
$`𝒮_{ZA}`$ $`=`$ $`8N_c(N_c+2+2N_O)V,`$ (41)
where $`w^{}`$ is some constant, $`N_c`$ is the number of colors and $`V`$ is the lattice volume. Since the effect of deflation is also similar to that in ORA, we do not repeat that discussion here.

The performance of the Zolotarev algorithm in numerical tests is summarized in Tables 8, 9 and 10. One needs to tune two algorithmic parameters, $`N_O`$ and $`ϵ_{CG}`$ for better efficiency. For a given value of $`ϵ_{CG}`$, we increase $`N_O`$ until a saturation in the value of error is evident. For $`N_O68`$, the performance is similar to that of the ORA. The improvement with increasing $`N_O`$ as $`\kappa `$ increases further indicates that the ZA and ORA should both improve if $`N_O`$ is allowed to change algorithmically with configuration. Such a method can be constructed from the results of wupp , when $`N_O`$ is increased until the maximum relative error (see Figure 3) attains a fraction ($`<`$1/2) of the desired accuracy. This can be implemented at the initialization step from the knowledge of the minimum and maximum of the relative error (defined in LHS of Eq.(40) or from the elliptic theta functions).

## VI Adaptive algorithm: conjugate gradient approximation

The first adaptive method used to compute $`M^{1/2}`$ was based on Lanczös algorithm borici . In the original suggestion, the number of Lanczös steps to be taken in order to reach a given precision was investigated in terms of the variation of the eigenvalues of $`M`$ with the number of Lanczös steps. A stopping criterion à la Conjugate Gradient was proposed but its relation to the precision was not direct. A related adaptive method based on the Conjugate Gradient algorithm was used in ragula . Here the stopping criterion is put on the residual vector in the inversion of $`M`$. This enables a direct control over the precision.

The CGA starts with an iteration which is almost the same as the usual CG algorithm for the inversion of $`M`$

1. Start from $`r_1=\mathrm{\Phi }`$, $`p_1=r_1`$ and $`\beta _1=0`$,
2. Iterate as in regular CG, $`\alpha _i=|r_i|^2/(p_i^{}Mp_i)`$, $`r_{i+1}=r_i\alpha _iMp_i`$, $`\beta _{i+1}=|r_{i+1}|^2/|r_i|^2`$, and $`p_{i+1}=\beta _{i+1}p_i+r_{i+1}`$.
3. Stop when $`|r_{i+1}|<ϵ_{CG}`$.

Note that the the only difference from the usual CG is that the vector which is $`M^1\mathrm{\Phi }`$ does not need to be obtained during the iteration.

In the orthonormal basis of $`q_i=r_i/|r_i|`$, the matrix $`M`$ is the composition of the matrix $`Q`$ whose $`i`$-th column is $`q_i`$, and a symmetric tridiagonal matrix, $`T`$,

$$M=Q^{}TQ,\mathrm{where}T_{ii}=\frac{1}{\alpha _i}+\frac{\beta _i}{\alpha _{i1}},\mathrm{and}T_{i,i+1}=\frac{\sqrt{\beta _{i+1}}}{\alpha _i},$$
(42)
where $`\alpha _i`$ and $`\beta _i`$ are defined in the iteration above. Then compute the eigenvalues and eigenvectors of this truncated tridiagonal matrix T,

$$T=U\mathrm{\Lambda }U^{}$$
(43)
where $`\mathrm{\Lambda }`$ is the diagonal matrix of the eigenvalues and $`U`$ the matrix of the eigenvectors in the basis $`Q`$. The CGA solution is

$$L[\mathrm{\Phi }]=Q^tU\mathrm{\Lambda }^{1/2}U^{}Q\mathrm{\Phi }/|\mathrm{\Phi }|.$$
(44)
The adaptability of the algorithm arises from the fact that we retain only the vectors $`q_i`$ which contribute significantly to the inverse of $`M`$ and we stop the iterations for $`i=N_{CG}`$ when $`|r_{N_{CG}+1}|<ϵ_{CG}`$.

The contribution to $`M^{1/2}`$ of the smallest eigenvalue $`1/\lambda _{\mathrm{min}}`$ of $`M`$ will be only $`1/\sqrt{\lambda _{\mathrm{min}}}`$. Since the stopping criterion $`|r_{i+1}|<ϵ_{CG}`$ is meant to compute $`M^1`$ it is more stringent than required. One can be more generous for $`M^{1/2}`$, and use instead the stopping criterion

$$|r_{i+1}|<ϵ_{CG}/\sqrt{\lambda _0^{(i)}}$$
(45)
where $`\lambda _0^{(i)}`$ is an upper bound of $`\lambda _{\mathrm{min}}`$. Fortunately a reasonable estimate can be obtained at each iteration $`i`$ without large overheads. For any tridiagonal matrix $`T`$ of order $`N_O`$, the number of eigenvalues greater than a fixed number $`\mu `$ is the number of positive values of $`d^{(j)}`$, where this set of numbers is defined by $`d^{(1)}=T_{11}\mu `$ and

$$d^{(j)}=T_{jj}\lambda _0^{(i)}(T_{j1,j})^2/d^{(j1)},$$
(46)
for $`2jN_O`$ book2 . An upper bound for $`\lambda _0^{(i)}`$ can always be fixed by searching for a number for which at least one of $`d^{(j)}`$ is non-positive. This can be done by bisection, starting from the initial estimate at the first step, $`\lambda _0^{(1)}=T_{11}`$. While this procedure increases the complexity by order $`\mathrm{log}N_{CG}`$, the new stopping criterion in Eq. (45) has two advantages over the usual CG stopping criterion— first, $`N_{CG}`$ is reduced and, second, the method becomes better adaptable since the observed $`ϵ_{1/2}`$ for a given $`ϵ_{CG}`$ becomes independent of $`\lambda _{min}`$.

Practically, to do the computation without storing the orthonormal basis $`Q`$, one makes $`N_{CG}`$ iterations to get the truncated matrix $`T`$, computes the matrix $`U`$ and the diagonal $`\mathrm{\Lambda }`$ using standard methods golub , and then repeats the $`N_{CG}`$ iterations to compute the solution $`L[\mathrm{\Phi }]`$. The most stringent restriction on the algorithm seems to be that one cannot use any pre-conditioning and must always start the iterations from $`p_1=r_1=\mathrm{\Phi }`$. This algorithm has only one parameter, $`ϵ_{CG}`$. The algorithm automatically adjusts the number of iterations to achieve the specified precision irrespective of the condition number. Thus, no configuration dependent tuning of algorithmic parameters is necessary when employing the CGA for QCD applications.

The complexity of the CGA is

$$𝒞_{CGA}2wN_{CG}V+\omega N_{CG}^{}{}_{}{}^{2}2w^{}V\mathrm{log}\left(\frac{1}{ϵ_{CG}}\right)\mathrm{log}\kappa $$
(47)
where $`\omega `$ is a number independent of $`V`$. The $`N_{CG}^{}{}_{}{}^{2}`$ term comes from the handling of the tridiagonal matrix, and can be neglected since $`N_{CG}V`$. The space complexity is the same as that of a standard CG—

$$𝒮_{CGA}8N_c(N_c+3)V.$$
(48)
Since the method is adaptive, no deflation is necessary. However, deflation reduces the condition number of the matrix, and hence could improve the complexity by reducing $`N_{CG}`$. Nevertheless, for reasons that we have already discussed in connection with CA and ORA, deflation is unlikely to improve the performance at fixed physics when taking the limit of large $`V`$.

The results of our numerical tests for this algorithm are collected in Tables 1113. Note that for all three test configurations there is a threshold in $`ϵ_{CG}`$ above which $`ϵ_{CC}10ϵ_{1/2}`$ and below which $`ϵ_{CC}`$ is roughly constant. The threshold value of $`ϵ_{CG}`$ is somewhat larger than $`\lambda _{\mathrm{min}}`$ for the configuration. Similar thresholds are also seen for the ORA and ZA. This behaviour possibly reflects the existence of a large unconverged subspace in the CG iterations.

## VII Comparing the algorithms

In Figure 4 we have collected different measures of performance of the four algorithms we investigated in this paper, namely, the Optimized Rational Approximation (ORA) nara , the Zolotarev Approximation (ZA, which is also a rational expansion) wupp ; ch , the Chebychev Approximation (CA, a polynomial expansion) hepile  and the Conjugate Gradient Approximation (CGA, an iterative method) ragula . Further details can be found in Tables 113.

It is clear that for modest values of the condition number of $`M`$, $`\kappa 10^3`$, the CA is the preferred algorithm. This is clear from the figure, as well as our results for the algorithmic complexities in Eqs. (28), (33), (41) and (47). However, with increasing condition numbers the performance of CA rapidly degrades. This is visible in the figure as well as in our analysis of the adaptability in Eq. (26). We have argued earlier that these drawbacks of the CA are generic to all polynomial expansions.

The ORA, in its present form with fixed $`N_O`$, also suffers from a lack of adaptability. In principle, this can be alleviated if the order of the approximation can be chosen adaptively. We have implemented the ZA, which is another rational approximation, for several different choices of order. As can be seen from Tables 8\- 10, and from Figure 4, this improves the performance tremendously. For comparable CPU times, corresponding to low order ZA, the performance is at least one order better than that of ORA on all configurations. The key to improving the performance of rational approximations is the automatic variation of the order $`N_O`$ with the condition numbers. In our tests we have simulated adaptability by working with several different orders and retained the one corresponding only to a small fraction of the inversion error. The so-tuned order is only slightly higher than that obtained in an automatic procedure defined at the end of section V.

The CGA depends on only one parameter $`ϵ_{CG}`$. For a given value of $`ϵ_{CG}`$, the corresponding errors $`ϵ_{1/2}`$ and $`ϵ_{GW}`$ are almost independent of the condition number of the matrix, thanks to the relaxed stopping criterion. The price for such a good adaptability is a computing time which is 50% higher than ZA for a given accuracy (70% excess if $`No`$ is small, 20% for large $`No`$). The price, however, ensures that for all the configurations one guarantees the same order of accuracy from a given value of $`ϵ_{CG}`$ and with a predicted value of $`ϵ_{GW}`$.

The variation in the number of conjugate gradient iterations, $`N_{CG}`$, as the stopping criterion, $`ϵ_{CG}`$, is changed for the three configurations is shown in Figure 5. The data for the ORA are not shown in the figure because they are very similar to those of the ZA. Note that the curves for the CGA lie below that for the ZA (despite the shift in ZA as compared to CGA), which is the influence of the relaxed stopping criterion discussed in Section VI. Ref. wupp  has devised a similar modification for ZA which can reduce $`N_{CG}`$ in that case. Note that our above results for ZA did not use any such modifications; using it will further enhance the performance of ZA reported above.

We have noted in Section II that the relations (18) between the errors are valid for those approximations to $`M^{1/2}`$ which commute with $`M`$. In particular, we noted that for the iterative algorithms these relations become valid, provided that $`ϵ_{CG}`$ is sufficiently small. In Figure 6 we demonstrate this for $`ϵ_H`$, which is expected to be zero when $`ϵ_{CG}`$ is small enough. For the CGA and ZA (data for the ORA are not shown because they almost coincide with that for ZA), $`ϵ_H`$ decreases with $`ϵ_{CG}`$. The slopes in this plot correspond to linear decrease when $`ϵ_{CG}`$ is sufficiently small. Clear non-linearities are present for larger $`ϵ_{CG}`$ when the condition number is large. We believe that these non-linearities are due to large non-converged subspaces, implying a need for high accuracy.

For fixed order algorithms the adaptability, $`𝒜`$, quantifies the configuration dependence of speed. The numerical study can be used more directly to illustrate the adaptability by studying the slowdown in going from configuration A to B (i.e., from $`\kappa =10^3`$ to $`4.4\times 10^5`$) or from A to C ($`\kappa `$ changes from $`10^3`$ to $`3.6\times 10^9`$). As shown in Figure 7, both ZA and CGA are adaptable algorithms over a wide range of $`ϵ_{1/2}`$. Since ZA is faster, as seen in Figure 4, it is thus the method of choice. Note, however, that CGA is very comparable to it, and may be preferred for its self-tuning ability.

We emphasise that a fair test of relative performance of algorithms is to work without deflation. First, deflation improves the performance of each of the algorithms we have investigated. Details are given in the sections on each algorithm. Nevertheless the algorithms based on rational approximation seems to be less sensitive to deflation than other ones because of the positive shifts introduced in the matrix. Second, since the computation of the eigensystem of $`M`$, necessary to deflation, is done at finite accuracy, it introduces extra errors. If the error in the computation of the eigenvalue $`\lambda _i`$ is $`\delta _i`$, then the contribution to $`ϵ_{1/2}`$ is $`\delta _i/\lambda _i`$. Thus, if we want to achieve a given $`ϵ_{1/2}`$, then we must keep $`|\delta _i|ϵ_{1/2}|\lambda _i|`$. When the condition number $`\kappa `$ increases, this criterion becomes impossible to satisfy, leading to catastrophic loss of accuracy.

We feel it worth pointing out that deflation is only one of many possible methods to decrease the effective condition number of the problem. Other preconditioning methods have not been seriously explored for overlap Fermions. The cost of accurate numerical methods seems to suggest that numerically stable preconditioning methods will pay a big dividend in this problem.

This work was supported by the Indo-French Centre for Promotion of Advanced Research under project number 2104-2.


format: text

text:
# FORWARD-BACKWARD MULTIPLICITY CORRELATIONS AND LEAKAGE PARAMETER BEHAVIOUR IN ASYMMETRIC HIGH ENERGY COLLISIONS

## 1 Introduction: the weighted superposition mechanism of different classes of events

It has been found quite recently  that the weighted superposition of different classes of events with negative binomial (NB) properties (events with and without mini-jets in $`pp`$ collisions, 2-, 3- or more-jets samples of events in $`e^+e^{}`$ annihilation) reproduces, in the GeV energy range, the available experimental data on forward-backward (FB) multiplicity correlations of the two types of collisions; it is an intriguing result which sheds new light on long range properties in multiparticle dynamics and outlines the experimental consistency of the weighted superposition mechanism of different classes of events. In addition, the same superposition mechanism has been shown to provide in the GeV region an interesting phenomenological tool in order to describe the observed shoulder effect in $`n`$-charged particle multiplicity distributions (MD’s) and $`n`$-oscillations in the ratios of $`n`$-particle factorial moments to the $`n`$-particle factorial cumulant moments . Effects on which QCD has not (up to now) satisfactory predictions .

In the case of $`pp`$ collisions, essential part of the theoretical background in the GeV region has been the assumption that, in each substructure or component described by a negative binomial (Pascal) multiplicity distribution \[NB(P)MD\], independently produced clans (they follow a Poisson distribution) are binomially distributed in the two hemispheres and that logarithmically produced charged particles within each clan distribute themselves again binomially in the two hemispheres with an energy independent ‘leakage parameter.’ This parameter controls the number of particles, generated by clans in one hemisphere, falling in the opposite one and was determined in  from the data at 63 GeV and 900 GeV. Accordingly, the correct reproduction in the GeV region of the experimentally observed increase with c.m. energy of the forward-backward multiplicity correlation strength as well as of the relation between the average charged particle multiplicity of particles lying in the backward hemisphere versus the number of charged particles lying in the forward hemisphere (and vice versa) support strongly our approach.

In the case of $`e^+e^{}`$ annihilations, in addition, it should be pointed out that under the same assumptions of binomial distributions of clans and of particles generated by clans in the forward and backward hemispheres, the relatively small value of the forward-backward correlation strength in the total charged particle MD and the absence of FB multiplicity correlations in the separate 2- and 3-jet samples of events, measured by OPAL collaboration , have been also correctly reproduced. Forward-backward correlations in the total charged particle multiplicity distribution are indeed only due here to the superposition effect of the two samples of events. In our approach, in fact, forward-backward multiplicity correlations in the two separate sample of events, each of NB type, turn out to be zero, as no leakage is found from one hemisphere to the other in the two separate samples, and in the total sample resulting from the weighted superposition of two NB(P)MD we predict a correlation which coincides with the experimental one within experimental error .

The striking difference in FB charged particle multiplicity correlation strength for the total sample of events between proton proton collisions in the GeV energy domain and $`e^+e^{}`$ annihilation at LEP energy outlines the deep link between FB multiplicity correlations and long range correlations, which are expected to be quite strong in the first case and relatively weak in the second one. In addition, the observed lack of FB multiplicity correlations in the two separate 2-jet and 3-jet samples of events in $`e^+e^{}`$ suggests that the weak FB correlations seen in this reaction in the total sample of events are entirely due to the superposition of the two separate samples. This fact can be considered indeed an experimental evidence of the weighted superposition mechanism characteristic of our approach (justified up to now only on the basis of a quite successful description of collective effects in high energy energy phenomenology which QCD is unable at present to describe) and of the presence or absence of charged particle leakage from one hemisphere to the opposite one.

In terms of clan structure analysis of the NB behaviour for each class of events, it is clear that charged particle leakage as well as FB multiplicity and long range correlations are expected to be stronger when particle population per clan is larger, a phenomenon which usually goes together with the reduction of the average number of clans.

It should be pointed out that the occurrence of a larger average number of smaller size clans and of a smaller average number of larger size clans has a suggestive interpretation at parton level, in the framework of a two-step mechanism, and could be related to smaller and larger colour flow densities respectively. In a NB description of final charged particle MD for a single component (soft or semi-hard in $`pp`$ collisions, 2- and 3-jet samples of events in $`e^+e^{}`$ annihilation), after using generalised local parton hadron duality (GLPHD), the first step of the parton production process is dominated by the $`A_{qq+g}(N_{\text{color}},ϵ)`$ vertex and the average number of partonic clans corresponds to the average number of bremsstrahlung gluon jets (BGJ), which are effective independent intermediate gluon sources (IIGS), and the second step is controlled by the vertex $`A_{gg+g}(N_{\text{color}},ϵ)`$, the gluon self-interaction vertex, whose increase corresponds to an enhancement of parton cascading from the IIGS . $`N_{\text{color}}`$ is the number of colours and $`ϵ`$ the fixed cut-off regularization prescription of the theory.

Coming to the difference between $`e^+e^{}`$ and $`pp`$ reactions, it is likely to expect that the observed increase of the average number of clans in $`e^+e^{}`$ with respect to the $`pp`$ case is a consequence of a stronger activity of the first vertex with respect to the second one, and that the opposite will occur for the increase of the average number of partons per partonic clans. It is quite clear that the production of a large average number of partons per clan is a consequence of a longer cascading process, originated by IIGS generated at relatively high virtuality in regions where, being the coupling constant smaller, stronger colour flow between partons should be at work (a situation favoured in $`pp`$ collisions). When the average number of partons per clan is relatively small, IIGS are expected to start to be effective at lower virtuality, their cascading becomes shorter and colour exchanges reduced with respect to the previous case (a situation favoured in $`e^+e^{}`$ annihilation). It seems therefore that the occurrence of stronger FB multiplicity correlations, long range correlations and related particle leakage enhancement from one hemisphere to the opposite one at hadron level are a specular image of stronger cascading from high virtuality IIGS and of larger colour exchange in this region at parton level.

We believe in fact that the understanding of FB multiplicity correlations at hadron level is a possible starting point in order to study new effects of colour quantum number exchanges in multiparticle dynamics.

Accordingly, we decided to continue our search initiated in Ref. , where FB correlations have been understood in the GeV region in symmetric reactions (like $`pp`$ collisions and $`e^+e^{}`$ annihilation) and for symmetric definition of the hemispheres by assuming at hadron level:

* NB behaviour for the (forward plus backward) MD of each component or substructure (class of events), i.e., clan structure analysis is assumed to be applicable to each of them. The generalisation to the class of compound Poisson distributions (CPD) is of course possible: the generating function $`G_{\text{CPD}}(z)`$ in this case can be written as follows

$$G_{\text{CPD}}(z)=\mathrm{exp}\left\{\overline{N}_g[g_c(z)1]\right\}$$
(1)
where $`\overline{N}_g`$ is the average number of generalised clans and $`g_c(z)`$ the generating function of the MD of charged particles originated by a generalised clan (e.g., it is a logarithmic distribution when the MD for a single component is a NB(P)MD).
* clans are therefore independently emitted and their distribution is Poissonian. In case of a generic CPD we talk, as previously stated, of generalised clans.
* clans are binomially distributed in the forward and backward hemispheres (as are generalised clans).
* logarithmically produced charged particles in each clan are also binomially distributed in the forward and backward hemispheres but with different probabilities $`p`$ and $`q`$ ($`p+q=1`$, $`p`$ different from $`q`$); $`p`$ controls the leakage from one hemisphere to the other: $`p=1`$ means that no particle leaks, while $`p`$ larger or equal to 0.5 and smaller than 1 indicates leakage. In case of a generic CPD, the $`n`$-charged particle MD is generated by $`g_c(z)`$, but $`p`$ and $`q`$ retain their meaning as described.

In order to perform our calculations in $`pp`$ collisions one extra assumption has been added, i.e., that the particle leakage parameter is constant throughout the GeV region for the two separate classes of events.

Possible scenarios in the TeV region based on extrapolations from data in the GeV region and on the weighted superposition mechanism of soft (without mini-jets) and semi-hard (with mini-jets) events have been indeed studied and predictions given on $`n`$-charged particle multiplicity distributions and on the ratios of $`n`$-particle factorial moments to the $`n`$-particle factorial cumulants moments general properties in the TeV region .

An attempt to predict the energy dependence of the forward-backward multiplicity correlation strength as well as $`\overline{n}_B(n_F)`$ vs $`n_F`$ general trends in the new energy range available at CERN with Alice detector has been considered . The word is now to experiment which is supposed to test all these three sets of predictions at 14 TeV.

The problem we want to face in this paper, on the theory side, is the generalisation of the approach discussed in Ref.  for determining forward-backward multiplicity correlations properties for a single component, in more complex asymmetric reactions (like heavy ion AB and $`p`$A collisions) and to provide a general framework for the study of forward-backward multiplicity correlations which includes symmetric reactions (like $`pp`$ and AA collisions and $`e^+e^{}`$ annihilation) as a particular case.

Corner stones of our argument remain of course assumptions a, b, c, d.

They could hardly be abandoned in view of their success in giving a good phenomenological description of available experimental data on FB multiplicity correlations in symmetric reactions.

The generalisation of the approach to asymmetric reactions will concern therefore mainly how to implement, in the framework defined by assumptions a.b.c.d, the asymmetry of the reaction and the asymmetric definition of the forward and backward hemispheres.

In view of the lack of sound experimental data on asymmetric reactions and of the related analyses in terms of two or more component substructures both in full phase space and in rapidity intervals, the present paper should be considered as a stimulus to experimentalists of the new generation machines for a deeper analysis of total $`n`$-charged particles MD and a more satisfactory understanding of forward-backward multiplicity correlations as the c.m. energy increases and at fixed c.m. energy within rapidity intervals.

## 2 The general asymmetric case

Generalisation of the study performed in Ref.  on FB multiplicity correlations consists at the present stage of investigation in assuming that

* the particle leakage percentage from the backward (B) to the forward (F) hemisphere, $`q_B`$, is different from the particle leakage percentage from F to B hemisphere, $`q_F`$, with $`p_F+q_F=1`$ and $`p_B+q_B=1`$, $`p_F`$ and $`p_B`$ being the corresponding percentages of particles not leaking from one hemisphere to the opposite one and remaining in the F and B hemisphere respectively; notice that in Ref. , $`p_F=p_B=p`$.
* Poissonianly generated clans are asymmetrically (but binomially) distributed in the two hemispheres with the asymmetry parameter $`r`$ different from 1/2, $`s`$ being its complement to 1 in the opposite hemisphere; in Ref. , notice that $`r=s=1/2`$ in addition to $`p_F=p_B=p`$.

### 2.1 The generating function

In general, the joint MD $`P(n_F,n_B)`$ for $`n_F`$ particles in the forward hemisphere and $`n_B`$ particles in the backward one is related to the global MD $`P(n)`$ though the probability distribution $`f(n_F|n)`$, which gives the probability to have $`n_F`$ F-particles when the total number of particles is $`n`$, as follows:

$$P(n_F,n_B)=P(n_F+n_B)f(n_F|n_F+n_B).$$
(2)
The generating function $`G(z_F,z_B)`$ for the joint distribution then satisfies

$$G(z_F,z_B)\underset{n_F,n_B}{}z_F^{n_F}z_B^{n_B}P(n_F,n_B)=\underset{n}{}z_B^nP(n)g_f(z_F/z_B;n),$$
(3)
where $`g_f(z;n)`$ is the generating function for $`f(n_F|n)`$; in case it is the binomial distribution with parameter $`p`$,

$$f(n_F|n)=\left(\genfrac{}{}{0pt}{}{n}{n_F}\right)p^{n_F}(1p)^{nn_F},$$
(4)
then its GF is (defining $`q1p`$):

$$g_f(z;n)=(q+pz)^n,$$
(5)
and one obtains a considerable simplification of Eq. (3):

$$G(z_F,z_B)=g(z_Fp+z_Bq),$$
(6)
where $`g(z)`$ is the generating function of the global distribution $`P(n)`$:

$$g(z)\underset{n}{}z^nP(n).$$
(7)
Accordingly, we proceed now to calculate the GF for the general case. The formulae below are heavily based on previous work . We consider a clan with logarithmic MD produced in the F hemisphere and assume that each particle has the same probability $`p_F`$ not to leak into the B hemisphere. Then the distribution at fixed number of particles is binomial, and we have immediately an application of the just-explained scheme: the GF for the joint distribution within a F-clan is thus

$$g_{c,F}(z_F,z_B)=g_{\text{log}}(z_Fp_F+z_Bq_F);$$
(8)
for a clan produced in the B hemisphere we have the corresponding GF:

$$g_{c,B}(z_F,z_B)=g_{\text{log}}(z_Fq_B+z_Bp_B).$$
(9)
Here the GF for the logarithmic distribution of parameter $`\beta `$ is

$$g_{\text{log}}(z)\mathrm{log}(1z\beta )/\mathrm{log}(1\beta ),$$
(10)
where, in terms of standard NB parameters $`\overline{n}`$ and $`k`$,

$$\beta =\frac{\overline{n}}{\overline{n}+k}$$
(11)
and is related to the average number of particles per clan via the formula:

$$\overline{n}_c=\frac{\beta }{(\beta 1)\mathrm{log}(1\beta )}.$$
(12)
Recalling now Eq. (34) of , it is easy to calculate the GF of the joint distribution of F and B particles at given numbers $`N_F`$ of F clans and $`N_B`$ of B clans:

$$g(z_F,z_B|N_F,N_B)=\left[g_{c,F}(z_F,z_B)\right]^{N_F}\left[g_{c,B}(z_F,z_B)\right]^{N_B}:$$
(13)
due to the fact that all clans are by definition independent from each other, we can convolute the respective MD’s, which corresponds to multiply together the GF’s.

In order to sum over the clan MD, let us remember that we have assumed that clans are Poisson distributed and independent of each other, thus the joint distribution can be written as:

$$𝒫(N_F,N_B)=\frac{\overline{N}^{N_F+N_B}}{(N_F+N_B)!}e^{\overline{N}}\left(\genfrac{}{}{0pt}{}{N_F+N_B}{N_F}\right)r^{N_F}(1r)^{N_B},$$
(14)
where $`\overline{N}=\overline{N}_F+\overline{N}_B`$ is the average number of clans, and $`r=\overline{N}_F/\overline{N}`$ is the fraction of clans emitted in the F hemisphere. The corresponding GF is

$$𝒢(z_F,z_B)=\mathrm{exp}\left\{\overline{N}[rz_F+(1r)z_B1]\right\}.$$
(15)
We can now perform the last step in the calculation, again exploiting the general properties of the binomial distribution:

$$g(z_F,z_B)=\underset{N_F}{}\underset{N_B}{}g(z_F,z_B|N_F,N_B)𝒫(N_F,N_B)=𝒢(g_{c,F}(z_F,z_B),g_{c,B}(z_F,z_B)).$$
(16)
$$g(z_F,z_B)=\mathrm{exp}\left\{r\overline{N}[g_{\text{log}}(z_Fp_F+z_Bq_F)1]\right\}\mathrm{exp}\left\{(1r)\overline{N}[g_{\text{log}}(z_Fq_B+z_Bp_B)1]\right\}.$$
(17)
Notice that the above formula is valid for the NB(P)MD, but can easily be extended, as anticipated, to any compound Poisson distribution, provided all correlations are exhausted within a (generalised) clan: it is sufficient to replace $`g_{\text{log}}`$ with the appropriate GF within a clan, $`g_c(z)`$.

### 2.2 The correlation strength and the marginal distributions

We can now use the above result to deduce the correlation strength $`b`$:

$$b=\frac{n_Fn_B\overline{n}_F\overline{n}_B}{\left[(n_F^2\overline{n}_F^2)(n_B^2\overline{n}_B^2)\right]^{1/2}},$$
(18)
since everything can be read from the GF:

$$n_Fn_B=\frac{^2g}{z_Fz_B}|_{z_F=z_B=1}$$
(19)
$$\overline{n}_i=\frac{g}{z_i}|_{z_F=z_B=1}n_i^2=\frac{^2g}{z_i^2}|_{z_F=z_B=1}+\overline{n}_i$$
(20)
with $`i=`$ F, B.

The average number of F-particles at fixed number of B-particles is also easily obtained though differentiation, since the MD of $`n_F`$ at fixed $`n_B`$ is given by

$$p(n_F|n_B)=\frac{p(n_F,n_B)}{p_{\text{marg}}^{(B)}(n_B)}=\frac{p(n_F,n_B)}{_{n_F}p(n_F,n_B)},$$
(21)
where $`p_{\text{marg}}^{(B)}(n_B)`$ is the marginal distribution in the B hemisphere; the GF is

$$g(z_F|n_B)\underset{n_F}{}z_F^{n_F}p(n_F|n_B)=\frac{{\displaystyle \frac{^{n_B}}{z_B^{n_B}}}g(z_F,z_B)}{{\displaystyle \frac{^{n_B}}{z_B^{n_B}}}g(1,z_B)}|_{z_B=0}.$$
(22)
being $`g(1,z_B)`$ the GF of $`p_{\text{marg}}^{(B)}`$. The first moment is obtained from the GF by differentiating once:

$$\overline{n}_F(n_B)\frac{g(z_F|n_B)}{z_F}|_{z_F=1}=\frac{{\displaystyle \frac{^{n_B}}{z_B^{n_B}}}{\displaystyle \frac{}{z_F}}g(z_F,z_B)}{{\displaystyle \frac{^{n_B}}{z_B^{n_B}}}g(z_F,z_B)}|_{z_B=0,z_F=1}.$$
(23)
Analogous formulae for the B hemisphere quantities can be easily obtained in the same way.

The corresponding marginal distributions are $`g(z_B=z,1)`$ and $`g(1,z_F=z)`$:

$$g(z,1)=\mathrm{exp}\left\{r\overline{N}[g_{\text{log}}(zp_F+q_F)1]\right\}\mathrm{exp}\left\{s\overline{N}[g_{\text{log}}(zq_B+p_B)1]\right\};$$
(24)
$`g(1,z)`$ can be obtained from Eq. (24) by interchanging parameter $`p_i`$ with $`q_i`$ (i=F,B).

The marginal distribution of Eq. (24) is the product of the GF’s of two NB(P)MD’s with characteristic NB parameters {$`\overline{n}rp_F`$ , $`kr`$} and {$`\overline{n}sq_B`$, $`ks`$} respectively as can be seen immediately by noticing that

$$g(z,1)=\left\{1+\frac{\overline{n}rp_F}{kr}(1z)\right\}^{kr}\left\{1+\frac{\overline{n}sq_B}{ks}(1z)\right\}^{ks};$$
(25)
$`g(1,z)`$ can be defined in an analogous way.

$`g(z,1)`$ and $`g(1,z)`$ are the GF’s of the MD’s obtained by convoluting the MD’s for particles generated by F (B) clans which stay in the F (B) hemisphere and do not leak in the opposite hemisphere with the MD’s for particles generated by B (F) clans which are leaking in the F (B) hemisphere.

Although $`g(z,1)`$ (and $`g(1,z)`$) is the product of the GF’s of NB(P)MD’s it is not the GF of a NB(P)MD. This consideration notwithstanding, it is interesting to remark that $`g(z,1)`$ (and $`g(1,z)`$) is an infinitely divisible distribution (IDD), a fact which allows to define the generalised clan concept, as already remarked.

In fact Eq. (24) can be rewritten as follows

$$g(z,1)=\mathrm{exp}\left\{\overline{N}[rg_{\text{log}}(zp_F+q_F)+sg_{\text{log}}(zq_B+p_B)]\overline{N}\right\}$$
(26)
i.e., as

$$g(z,1)=\mathrm{exp}\left[\overline{N}(A(z)1)\right]$$
(27)
with $`A(z)=[rg_{\text{log}}(zp_F+q_F)+sg_{\text{log}}(zq_B+p_B)]`$ and $`A(0)`$ different from zero. This last remark implies that the probability of generating zero particles is different from $`e^{\overline{N}}`$. Therefore the probability of generating empty clans is also different from zero, a result in contrast with the standard clan definition (each clan contains at least one particle), definition which we would like of course to enforce.

In order to do that, let us add and subtract in the exponent of Eq. (27) the term $`A(0)`$, i.e., we rewrite Eq. (27) as follows

$$g(z,1)=\mathrm{exp}\left\{\overline{N}[A(z)A(0)+A(0)1]\right\}.$$
(28)
From Eq. (28) one obtains the compound Poisson distribution belonging to the class of IDD

$$g(z,1)=\mathrm{exp}\left\{\overline{N}_g[G_g(z)1]\right\},$$
(29)
with $`\overline{N}_g=[1A(0)]\overline{N}`$ and $`G_g(z)=[A(z)A(0)][1A(0)]^1`$.

Accordingly, one can define the average number of generalised clans, $`\overline{N}_g`$, in terms of the standard NB parameters $`\overline{n}`$ and $`k`$, and of $`p_F`$ and $`p_B`$:

$$\begin{array}{cc}\hfill \overline{N}_g& =\overline{N}[1A(0)]=\overline{N}\{1[rg_{\text{log}}(q_F)+sg_{\text{log}}(p_B)]\}\hfill \\ & =\mathrm{ln}g(0,1)=rk\mathrm{ln}(\overline{n}+k\overline{n}q_F)+sk\mathrm{ln}(\overline{n}+k\overline{n}p_B)k\mathrm{ln}k.\hfill \end{array}$$
(30)
The backward marginal multiplicity distribution and related properties can be be obtained from the forward one by interchanging parameter $`p_i`$ with $`q_i`$ (i=F,B).

## 3 The cases of partially removed symmetry

### 3.1 Symmetry for clans only ($`r=s=1/2`$, $`p_Fp_B`$, $`q_Fq_B`$)

The symmetry of the reaction can be broken partially by assuming that binomially distributed clans go fifty per cent in the F hemisphere and fifty per cent in the B hemisphere ($`r=s=1/2`$) but particle leakage from clans in F to B hemisphere, $`q_F`$, and particle leakage from clans in B to F hemisphere, $`q_B`$, (as their complements $`p_B`$ and $`p_F`$ with $`p_F+q_F=1`$ and $`p_B+q_B=1`$) are different. In other words, assumption B from the previous section is valid but A is not.

The forward-backward joint charged particle multiplicity GF for one single component turns out to be

$$g(z_B,z_F)=\mathrm{exp}\left\{\frac{1}{2}\overline{N}[g_{\text{log}}(z_Fp_F+z_Bq_F)1]\right\}\mathrm{exp}\left\{\frac{1}{2}\overline{N}[g_{\text{log}}(z_Bp_B+z_Fq_B)1]\right\},$$
(31)
and for the corresponding marginal B distribution ($`z_B=z`$, $`z_F=1`$)

$$g(z,1)=\mathrm{exp}\left\{\frac{1}{2}\overline{N}[g_{\text{log}}(p_F+zq_F)1]\right\}\mathrm{exp}\left\{\frac{1}{2}\overline{N}[g_{\text{log}}(zp_B+q_B)1]\right\}.$$
(32)
In terms of standard NB parameters Eq. (32) becomes

$$\begin{array}{cc}\hfill g(z,1)& =\left[1+\frac{\overline{n}p_F}{k}(1z)\right]^{k/2}\left[1+\frac{\overline{n}q_B}{k}(1z)\right]^{k/2}\hfill \\ & =\left[1+\frac{\overline{n}}{k}(1z)+\frac{\overline{n}^2}{k^2}p_Fq_B(1z)\right]^{k/2}.\hfill \end{array}$$
(33)
The symmetric definition of the two hemispheres is removed here in one component by assuming that particle leakage parameters in the forward and backward hemispheres are different.

The GF of the forward marginal multiplicity distribution thanks to the quadratic term as already shown in the general case is no more a NB(P)MD, although it is the product of the GF’s of two NB(P)MD’s with characteristic parameters {$`\frac{1}{2}\overline{n}p_F`$, $`\frac{1}{2}k`$} and {$`\frac{1}{2}\overline{n}q_B`$, $`\frac{1}{2}k`$} with asymmetric average charged particle multiplicities. The GF of the marginal distribution becomes the GF of a NB(P)MD in some special cases, i.e., for $`q_B=0`$, corresponding to no particle leakage from the backward to the forward hemisphere, and for $`\overline{n}k`$, a situation which occurs in the Poissonian limit for average charged particle multiplicity much less than the $`k`$ parameter and almost coinciding with the average number of clans. The GF $`g(z,1)`$, although not NB, is still an infinitely divisible distribution, as expected.

### 3.2 Symmetry for particles within clans only ($`rs1/2`$, $`p_F=p_B`$, $`q_F=q_B`$)

Another way to remove partially the symmetry is to use assumption A without B, i.e., to assume that binomially distributed clans are not symmetrically subdivided between F and B hemispheres, but particle leakage from clans in one hemisphere to the other is the same ($`rs1/2`$ with $`r+s=1`$, and $`p_F=p_B=p`$, $`q_F=q_B=q`$ with $`p+q=1`$, and $`p`$ larger or equal than 1/2 and smaller than 1).

The FB joint particle multiplicity distribution GF for one component, $`g(z_F,z_B)`$, becomes in this case

$$g(z_F,z_B)=\mathrm{exp}\left\{\overline{N}r[g_{\text{log}}(z_Fp+z_Bq)1]\right\}\mathrm{exp}\left\{\overline{N}s[g_{\text{log}}(z_Bp+z_Fq)1]\right\},$$
(34)
and the corresponding forward marginal charged particle MD GF

$$g(z,1)=\mathrm{exp}\left\{\overline{N}r[g_{\text{log}}(zp+q)1]\right\}\mathrm{exp}\left\{\overline{N}s[g_{\text{log}}(zq+p)1]\right\},$$
(35)
which in terms of standard NB parameters of the total charged particle MD of the component under investigation is

$$g(z,1)=\left[1+\frac{\overline{n}p}{k}(1z)\right]^{kr}\left[1+\frac{\overline{n}q}{k}(1z)\right]^{ks}.$$
(36)
The forward marginal MD GF, although again not of NB type, is the product of the GF’s of two NB(P)MD’s with parameters {$`\overline{n}rp`$, $`kr`$} and {$`\overline{n}qs`$, $`ks`$} respectively. The two sets of parameters characterise the group of non-leaking particles lying in the F hemisphere and those which leak from the backward hemisphere to the forward one. The backward marginal distribution can be easily obtained from the above equations by interchanging parameter $`r`$ with $`s`$.

## 4 The symmetric case ($`r=s=1/2`$, $`p_F=p_B`$, $`q_F=q_B`$)

In this case, assumptions A and B are both rejected. The formula studied in  for the joint charged particle MD GF in the symmetric case follows from Eq. (17) by taking $`p_F=p_B=p`$, $`q_F=q_B=q`$ with $`p`$ different from $`q`$, $`p+q=1`$, and $`r=s=1/2`$.

We get

$$\begin{array}{cc}\hfill g(z_B,z_F)& =\mathrm{exp}\left\{\frac{1}{2}\overline{N}[g_{\text{log}}(z_Fp+z_Bq)1][g_{\text{log}}(z_Bp+z_Fq)1]\right\}\hfill \\ & =k\frac{[k+\overline{n}(1pz_Bqz_F)]^{k/2}}{[k+\overline{n}(1pz_Fqz_B)]^{k/2}}\hfill \end{array}$$
(37)
and for the corresponding forward marginal distribution $`g(z,1)`$

$$\begin{array}{cc}\hfill g(z,1)& =[1+\frac{\overline{n}q}{k}(1z)]^{k/2}[1+\frac{\overline{n}p}{k}(1z)]^{k/2}\hfill \\ & =\left\{1+\frac{\overline{n}}{k}(1z)+pq\left[\frac{\overline{n}}{k}(1z)\right]^2\right\}^{k/2}.\hfill \end{array}$$
(38)
In conclusion, Eq. (17) can be considered the wanted generalisation of Eq. (37) to the case of asymmetric reactions and an asymmetric definition of the forward and backward hemispheres.

## 5 Behaviours of $`\overline{n}_F(n_B)`$ vs. $`n_B`$ and of $`\overline{n}_B(n_F)`$ vs. $`n_F`$

In this section we examine the relation between the average number of particles in one hemisphere, $`\overline{n}_F(n_B)`$ or $`\overline{n}_B(n_F)`$, at fixed value of the number of particles, respectively $`n_B`$ or $`n_F`$, in the other hemisphere, according to Eq. (23).

Let us start by recalling that the von Bahr-Ekspong theorem  implies, since we assume the total MD to be of NB type, that there is no linearity in $`\overline{n}_B(n_F)`$ vs $`n_F`$ (and in $`\overline{n}_F(n_B)`$ vs $`n_B`$) unless the MD for F-particles (and for B-particles!) at fixed total number of particles is binomial, in which case the GF is simply:

$$\mathrm{exp}[g_{\text{log}}(z_Fp+z_Bq)1].$$
(39)
Notice that it does not make sense to distinguish $`p_F`$ from $`p_B`$ in this case; however, if one did not distinguish F-clans from B-clans, i.e., if each clan emitted the same fraction $`p`$ in one hemisphere, then in Eqs. (8) and (9) above one would put $`p_F=q_B=p`$ and Eq. (17) would reduce to Eq. (39).

In the following the symmetry in the clan distribution is contrasted with the asymmetric clan distribution. In addition, leakage parameters for the two components are taken either equal or different. The resulting pictures of forward-backward multiplicity correlations are shown in Fig.s 1 and 2: of course our choice of the involved parameters is arbitrary. The word is again to experiment on symmetric collisions ($`r=s=1/2`$) in order to understand eventual deviations from identity of particle leakage parameters from clans in the forward and backward hemispheres and to experiment on asymmetric collisions ($`rs1/2`$) in order to measure the different leakage parameters in the two hemispheres. All together, expected different values of $`r`$ and $`s`$ as well as $`p_F`$ and $`p_B`$ parameters could lead to a new intriguing classification, in terms of FB multiplicity correlations, of high energy collisions and their substructures.

### 5.1 The case $`r=s=1/2`$ (symmetry in the clan distribution)

We start by examining the case in which $`r=1/2`$, illustrated in Figure 1. When $`p_F=p_B=1/2`$, there is perfect linearity in $`\overline{n}_B(n_F)`$ versus $`n_F`$. This is in agreement with the mentioned theorem, because here we are saying that, within each clan, particles are binomially distributed in F and B with the same probability 1/2, thus the fact that particles are produced in clans does not make any difference. When $`1/2<p_F=p_B<1`$, on the contrary, the fact that the production happens in two steps becomes again important and the relation between $`\overline{n}_B(n_F)`$ and $`n_F`$ is non-linear. Linearity is recovered again in the limiting case of no correlations ($`p_F=p_B=1`$), see below.

When $`p_Fp_B`$, the symmetry is lost because the particles leaking from the F hemisphere into B are not compensated (if $`p_F<p_B`$), or are over-compensated (if $`p_F>p_B`$), by those leaking from B into F. Thus one obviously finds that $`\overline{n}_B(n_F)`$ vs $`n_F`$ is not the same as $`\overline{n}_F(n_B)`$ vs $`n_B`$. Indeed the two marginal distributions differ from each other, but can be obtained one from the other by exchanging $`p_F`$ with $`p_B`$. In general for $`p_F>p_B`$ there is less leakage from F than from B: particles prefer to stay in the F hemisphere thus one has the line of $`\overline{n}_F(n_B)`$ vs $`n_B`$ above $`\overline{n}_B(n_F)`$ vs $`n_F`$. Vice versa, the opposite is true for $`p_F<p_B`$.

One can notice that when $`p_F=1`$, $`\overline{n}_F(n_B)`$ vs $`n_B`$ is a straight line, and when $`p_B=1`$, $`\overline{n}_B(n_F)`$ vs $`n_F`$ is linear. The fact is that when $`p_F=1`$, clans that fall in the F hemisphere do not contribute to the B hemisphere. The only relation between $`n_B`$ and $`\overline{n}_F`$ is given by B clans leaking to the F hemisphere. And this relation is by definition binomial for each clan: the number of clans is Poisson distributed, thus the overall is still binomial with the same parameter, as was shown above. Indeed, the GF with $`p_F=1`$ is:

$$g(z_F,z_B)=\mathrm{exp}\left\{r\overline{N}[g_{\text{log}}(z_F)1]\right\}\mathrm{exp}\left\{(1r)\overline{N}[g_{\text{log}}(z_Fq_B+z_Bp_B)1]\right\};$$
(40)
applying Eq. (23), the first term gives just the constant $`r\overline{n}`$, while the second one is of the same type of Eq. (39) and thus also gives a linear behaviour. Notice that this happens independently of the value of $`r`$. A corresponding result can be obtained for $`p_B=1`$ by appropriately exchanging the roles of F and B.

### 5.2 The case $`rs1/2`$ (asymmetry in the clan distribution)

When $`r`$ increases over 1/2 (e.g., in Figure 2 the case $`r=3/4`$ is illustrated), the excess of clans in the F hemisphere implies an increase of $`\overline{n}_F(n_B)`$ vs $`n_B`$ and a decrease of $`\overline{n}_B(n_F)`$ vs $`n_F`$, enlarging the differences between the curves when $`p_F>p_B`$ and reducing them when $`p_F<p_B`$ ($`\overline{n}_F(n_B)`$ vs $`n_B`$ can even become larger than $`\overline{n}_B(n_F)`$ vs $`n_F`$, depending on the values of $`r`$, $`p_F`$ and $`p_B`$). Exchanging $`p_F`$ with $`p_B`$ does not restore the symmetry. If $`p_F=p_B=1/2`$, then the value of $`r`$ has no influence on the outcome: if particles within each clan have a 50% chance of going into the other hemisphere, it makes no difference if clans prefer one hemisphere over the other, and the two curves coincide again. This is not true if $`p_F=p_B>1/2`$. Finally, in agreement with the remark at the end of the last subsection, we notice that even in the case $`r>1/2`$ when $`p_F=1`$, $`\overline{n}_F(n_B)`$ vs $`n_B`$ is a straight line, and when $`p_B=1`$, $`\overline{n}_B(n_F)`$ vs $`n_F`$ is linear.

### 5.3 Asymmetry in the average number of particles per clan

All the above has been studied in the framework of a NB(P)MD describing in one component the forward plus backward MD. However, there are low energy data on pA collisions in which NB behaviour has been found to hold separately in the forward and in the backward hemisphere, with different NB parameters . Because clans are independently emitted, the only correction to make to our formulae is to allow for different average numbers of particles per clan in the two hemispheres: in other words, the parameter $`\beta `$ of the logarithmic distributions in Eq.s (8) and (9) will be different, say $`\beta _F`$ and $`\beta _B`$. Now, we immediately get

$$g(z_F,z_B)=\mathrm{exp}\left\{r\overline{N}[g_{\text{log}}(z_Fp_F+z_Bq_F;\beta _F)1]\right\}\mathrm{exp}\left\{s\overline{N}[g_{\text{log}}(z_Fq_B+z_Bp_B;\beta _B)1]\right\}$$
(41)
in place of Eq. (17). The product of these two NB(P)GF is not a NB(P)MD, but it is still, by construction, an infinitely divisible distribution.

## 6 Conclusions

It has been shown that assuming different particle leakage percentages ($`p_Bp_F`$) for binomially generated particles from clans in one hemisphere to the opposite one and asymmetric ($`rs`$) distribution in the two hemispheres of binomially generated clans, a general formula for the generating function of the joint $`(n_F,n_B)`$-charged particle multiplicity distribution for each class of events (or substructure) can be obtained when the total MD GF is of NB type. The formula reduces to that discussed in Ref.  for $`p_B=p_F`$ and $`r=s=1/2`$. Of particular interest are also the cases in which the symmetry is only partially removed (assumption A without B and B without A). All above-mentioned results, although explicitly derived for substructures of NB type, can be easily generalised to any discrete infinitely divisible MD. This search is relevant for the study of forward-backward multiplicity correlations in non-identical heavy ion and in proton-nucleus collisions. Accordingly, the newly introduced particle leakage and asymmetry parameters can be considered as effective indeces classifying different classes of collisions.


text:
# Unphysical Predictions of Some Doubly Special Relativity Theories

## Abstract

A kind of doubly special relativity theory proposed by J. Magueijo and L. Smolin \[Phys. Rev. Lett. 88, 190403 (2002)\] is analysed. It is shown that this theory leads to serious physical difficulties in interpretation of kinematical quantities. Moreover, it is argued that statistical mechanics and thermodynamics cannot be resonably formulate within the model proposed in the mentioned paper.

Recently J. Magueijo and L. Smolin Magueijo and Smolin (2002) have proposed a modification of special relativity with intrinsically built in Planck length. Their construction is based on a non-linear form of the Lorentz group action in the momentum space. This is related to a special choice of the Lie algebra of the Lorentz group in the enveloping algebra of the nilpotent Lie algebra generated by $`/p_\mu `$, $`p_\mu `$ and $`I`$. Accordingly to Magueijo and Smolin (2002) the Lorentz boosts in the $`x`$-direction read

$`p_0^{}`$ $`={\displaystyle \frac{p_0\mathrm{cosh}\xi +p_1\mathrm{sinh}\xi }{1+lp_0(\mathrm{cosh}\xi 1)+lp_1\mathrm{sinh}\xi }},`$ (1a)
$`p_1^{}`$ $`={\displaystyle \frac{p_1\mathrm{cosh}\xi +p_0\mathrm{sinh}\xi }{1+lp_0(\mathrm{cosh}\xi 1)+lp_1\mathrm{sinh}\xi }},`$ (1b)
$`p_2^{}`$ $`={\displaystyle \frac{p_2}{1+lp_0(\mathrm{cosh}\xi 1)+lp_1\mathrm{sinh}\xi }},`$ (1c)
$`p_3^{}`$ $`={\displaystyle \frac{p_3}{1+lp_0(\mathrm{cosh}\xi 1)+lp_1\mathrm{sinh}\xi }},`$ (1d)
where $`l`$ is interpreted as the Planck length $`l_P`$. The above transformations together with rotations close to the Lorentz group and leave unchanged the following invariant built from the particle four-momentum $`p_\mu `$:

$$m^2c^2=\frac{\eta ^{\mu \nu }p_\mu p_\nu }{(1lp_0)^2},$$
(2)
when $`c`$, as usually, denotes the speed of light.

The scheme proposed in Magueijo and Smolin (2002) belongs to the class of recently investigated models known as “doubly special relativity” Magueijo and Smolin (2002); Amelino-Camelia and Piran (2001); Amelino-Camelia (2001a, b, 2002); Lukierski and Nowicki ; Kowalski-Glikman and Nowak (a, b). In this paper we analyze this scheme in details and we show that it is plagued by physically unacceptable features.

To present our discussion as simply as possible, in the following we restrict ourselves to the two-dimensional case. We will denote the energy and momentum as $`cp_0`$ and $`p`$, respectively. Furthermore, to separate the dimensional and scale parameters of the model we use the quantity $`\lambda /m_Pc`$ instead of $`l`$, where $`\lambda `$ is dimensionless and $`m_P=\sqrt{\mathrm{}c/G_N}`$ is the Planck mass. It is enough to consider positive values of $`\lambda `$. Therefore the Lorentz transformations (1) take the form

$`p_0^{}`$ $`={\displaystyle \frac{p_0\mathrm{cosh}\xi +p\mathrm{sinh}\xi }{1+\frac{\lambda }{m_Pc}p_0(\mathrm{cosh}\xi 1)+\frac{\lambda }{m_Pc}p\mathrm{sinh}\xi }},`$ (3a)
$`p^{}`$ $`={\displaystyle \frac{p_0\mathrm{sinh}\xi +p\mathrm{cosh}\xi }{1+\frac{\lambda }{m_Pc}p_0(\mathrm{cosh}\xi 1)+\frac{\lambda }{m_Pc}p\mathrm{sinh}\xi }},`$ (3b)
while the invariant (2) reads

$$\frac{p_0^2p^2}{\left(1\frac{\lambda }{m_Pc}p_0\right)^2}=m^2c^2.$$
(4)
Let us consider firstly the possible solutions of (4). It follows that we have singularity at $`p_0=m_Pc/\lambda `$, so possible range of $`p_0`$ is restricted to the two distinct areas: $`p_0<m_Pc/\lambda `$ or $`p_0>m_Pc/\lambda `$.

Now, for $`p_0m_Pc/\lambda `$ Eq. (4) determines algebraic conics; namely

$`\text{hyperbola},`$ $`\text{when }m`$ $`<m_P/\lambda ,`$ (5a)
$`\text{parabola},`$ $`\text{when }m`$ $`=m_P/\lambda ,`$ (5b)
$`\text{ellipse},`$ $`\text{when }m`$ $`>m_P/\lambda ,`$ (5c)
$`\text{the pair of half-lines},`$ $`\text{when }m`$ $`=0.`$ (5d)
Under physical condition that for the free motion the energy $`cp_0`$ is bounded from below we obtain four possible situations showed in the Fig. 1.

The line $`p_0=m_Pc/\lambda `$ divides each curve into distinct parts: upper and lower ones. The momentum manifold must be an orbit of the Lorentz group realized as in (3). We can easily prove that this is possible only for the lower parts of conics showed in the Fig. 1(a,b,d); on the upper parts of these conics the Lorentz group cannot be globally realized. Moreover, the upper part of the curve showed in the Fig. 1(c), despite of the fact that it forms an orbit of the Lorentz group, is also unphysical because the energy of the particle reaches the maximum when its momentum vanishes. Furthermore, if we restricts ourselves to these physically acceptable parts, denoted in the Fig. 1 by the solid lines, we can linearize the transformations (3) by the following map Lukierski and Nowicki

$$k_\mu =\frac{p_\mu }{1\frac{\lambda }{m_Pc}p_0},$$
(6)
where the momentum components $`k_\mu `$ form the usual Minkowski vector. Note, that the Jacobian determinant of this transformation vanishes only for $`p_0=m_Pc/\lambda `$. under the action of the Lorentz group $`k_\mu `$’s transform linearly

$`k_0^{}`$ $`=k_0\mathrm{cosh}\xi +k\mathrm{sinh}\xi ,`$ (7a)
$`k^{}`$ $`=k\mathrm{cosh}\xi +k_0\mathrm{sinh}\xi ,`$ (7b)
and leaves the invariant $`k_0^2k^2=m^2c^2`$ unchanged. Note, that the minima of $`p_0`$ and $`k_0`$ (i.e. the masses) are related by $`\mu =m/(1+\lambda m/m_P)`$.

Concluding the above discussion: The nonlinear transformation law (3) is a consequence of the choice of a nonlinear coordinate system in the momentum space and it is not *essentially nonlinear realization of the Lorentz group* 111Strictly speaking, the nonlinearity is an illusion in this case: As it is well known Coleman et al. (1969); Callan et al. (1969) essentially nonlinear group realizations are connected with the group action in the homogeneous coset spaces $`G/H`$, where $`H`$ is a subgroup of the group $`G`$., because it can be linearized by the appropriate choice of coordinates (6).

The question if the momentum coordinates $`p_0`$ and $`p`$ are physically admissible is open. In particular we do not have the notion of the inertial frame (observer) defined operationally. Therefore, up to now, we cannot properly answer this question (in contrast to the statement in e.g. Magueijo and Smolin (2002); Judes and Visser ). However, if we agree that canonical formalism can be used in such a simple kinematical problem (i.e. the free motion), we can try to partially answer this question. According to Magueijo and Smolin (2002); Judes and Visser  let us identify the energy with the generator of translation in time (Hamiltonian), i.e. $`H=cp_0`$, while $`p`$ with the canonical momentum. Therefore the particle velocity $`v`$ can be calculated from one of the Hamilton equations

$$v=\frac{H}{p}=\frac{pc}{p_0\left(1\frac{\lambda ^2m^2}{m_P^2}\right)+\frac{\lambda ^2m^2c}{m_P}}.$$
(8)
Now, with help of (3), (4) and (8) we are able to find the transformation law for the particle velocity under the action of the Lorentz group

$$v^{}=\frac{c\left[\frac{v}{c}\left(\mathrm{cosh}\xi +\frac{v}{c}\mathrm{sinh}\xi \right)+\frac{m_P}{\lambda m}\sqrt{1\frac{v^2}{c^2}\left(1\frac{\lambda ^2m^2}{m_P^2}\right)}\left(\frac{v}{c}\mathrm{cosh}\xi +\mathrm{sinh}\xi \right)\right]}{1\frac{v^2}{c^2}+\frac{v}{c}\left(\frac{v}{c}\mathrm{cosh}\xi +\mathrm{sinh}\xi \right)+\frac{m_P}{\lambda m}\sqrt{1\frac{v^2}{c^2}\left(1\frac{\lambda ^2m^2}{m_P^2}\right)}\left(\mathrm{cosh}\xi +\frac{v}{c}\mathrm{sinh}\xi \right)}.$$
(9)
The transformation law (9) goes to the standard Lorentz transformation law for velocities when $`\lambda 0`$ or when $`m=0`$ ($`\mu =0`$) and $`v=c`$. However, for $`\lambda 0`$ and $`m0`$ it depends on the particle mass $`\mu `$ (recall that $`m=\mu /(1\lambda \mu /m_P)`$). Such a situation is obviously in conflict with our physical space-time intuition: Indeed, let us imagine two bodies $`A`$ and $`B`$ with mass $`m_A`$ and $`m_B`$, respectively, moving in an inertial frame with the same velocity $`v`$ (for simplicity one can assume $`v=0`$). From the point of view of Lorentz boosted observer (i.e. from another inertial frame), if $`m_Am_B`$ and $`\lambda 0`$, these two bodies *have different velocities*! This is very undesirable feature from the physical point of view.

To understand better this issue let us look for the relationship of $`v`$ and the Lorentz velocity (i.e. the velocity which appears in the usual, *linear*, Lorentz transformations) $`v_L=ck/k_0`$. By means of (6) and (8) we get that

$$v=\frac{v_L}{1+\frac{\lambda m}{m_P}\sqrt{1\frac{v_L^2}{c^2}}}.$$
(10)
The velocity $`v_L`$ has the proper physical interpretation as the velocity canonically related to $`k`$, i.e. $`v_L=k_0/k=k/k_0`$. For the other hand it is defined as $`v_L=dx/dt`$, where $`t`$ and $`x`$ are usual Minkowskian time and coordinate. However, $`v_L`$ cannot be canonically related to the momentum $`p`$. Notice that $`v_L=k/k_0=p/p_0`$ and it has the standard transformation law.

Concluding, we have serious problems with the velocity transformation law and with the definition of the inertial observers within the model presented in Magueijo and Smolin (2002).

The next difficulty occurs when we try to formulate statistical mechanics and thermodynamics within the framework of theory proposed by Magueijo and Smolin Magueijo and Smolin (2002). Let us return to the three-dimensional case. It is easy to see that the invariant measure in the momentum space is given by the following formula

$$d\mathrm{\Gamma }=\frac{\left(1\frac{\lambda ^2m^2}{m_P^2}\right)^3d^3p}{2\left[1\frac{\lambda }{m_Pc}\sqrt{m^2c^2+|\stackrel{}{p}|^2\left(1\frac{\lambda ^2m^2}{m_P^2}\right)}\right]^3\sqrt{m^2c^2+|\stackrel{}{p}|^2\left(1\frac{\lambda ^2m^2}{m_P^2}\right)}},$$
(11)
(this holds also for the case $`m=\mu =0`$) and, since the momentum space is bounded (see Fig. 1), the one-particle partition function is

$$Z_1=V_{\frac{m_Pc}{\lambda }}^{\frac{m_Pc}{\lambda }}𝑑\mathrm{\Gamma }e^{\beta E}.$$
(12)
But the measure (11) is singular at $`|\stackrel{}{p}|m_Pc/\lambda `$ and the energy approaches the limit $`m_Pc^2/\lambda `$ when $`|\stackrel{}{p}|m_Pc/\lambda `$, so the integral (12) is divergent. Consequently the partition function does not exist in the considered case, as well as the internal energy and the entropy. The same holds in arbitrary $`N`$-particle case.

If we take non-relativistic measure in the momentum space $`d^3p`$ instead of $`d\mathrm{\Gamma }`$, the integral (12) becomes convergent, but on the other hand the energy of the free gas of $`N`$ identical particles (either massive or massless) is equal to Judes and Visser

$$E=m_Pc^2\frac{_{i=1}^N\frac{E_i}{m_Pc^2\lambda E_i}}{1+_{i=1}^N\frac{\lambda E_i}{m_Pc^2\lambda E_i}}$$
(13)
and

$$E\frac{m_Pc^2}{\lambda }$$
(14)
for large $`N`$ 111Note, that for the choice $`\lambda =1`$, as it is assumed in Magueijo and Smolin (2002), the mass of any large enough system is nearly equal to the Planck mass, i.e. to $`2.176\times 10^8\mathrm{kg}`$. On the other hand, there is no other natural physical choice for the parameter $`\lambda `$.. The partition function $`Z`$ of the $`N`$-particle free gas is then of the form

$$Ze^{\beta \frac{m_Pc^2}{\lambda }}\left[8\left(\frac{m_Pc}{\lambda }\right)^3V\right]^N$$
(15)
and, consequently, we come at the curious conjecture that the internal energy in the thermodynamical limit (i.e. for $`N\mathrm{}`$)

$$U=\frac{\mathrm{ln}Z}{\beta }=\frac{m_Pc^2}{\lambda }$$
(16)
does not depend on temperature!

Concluding this point of the discussion we have to state that statistical mechanics and/or thermodynamics do not exist within the model presented in Magueijo and Smolin (2002).

We would like to point out that the non-additivity of energy seems to be a common feature of doubly special relativity theories (see Judes and Visser ). Therefore, the question of the possiblity of formulation of statistical mechanics or thermodynamics in such theories is still open 222In the paper Kowalski-Glikman  only the one-particle partition function for the massless case is discussed within the model proposed in Amelino-Camelia and Piran (2001); Amelino-Camelia (2001a, b, 2002)..

We can conclude that the model proposed by Magueijo and Smolin Magueijo and Smolin (2002) encounters some serious difficulties discussed above. First of all, they are interpretational problems of space-time quantities: coordinates, velocities, etc. Moreover, it seems to be impossible to formulate resonably statistical mechanics and thermodynamics within this framework.

###### Acknowledgements.

The Authors are grateful to J. Lukierski and J. Jędrzejewski for interesting discussions. This work is supported by Łódź University grant No. 267.