text:
# A new class of binning-free, multivariate goodness-of-fit tests: the energy tests11footnote 1Supperted by Bunderministerium für Bildung und Forschung, Deutschland

## 1 Introduction

There exist a large number of binning-free goodness-of-fit tests for univariate distributions. The extension to the multivariate case is difficult, because there a natural ordering system is missing. The popular $`\chi ^2`$ test on the other hand suffers from the arbitrariness of the binning and from lack of power for small samples.

The test proposed in this article avoids both, ordering and binning. It is especially powerful in the multivariate case. It is called energy test, because the definition of the test statistic is closely related to the energy of electric charge distributions. To facilitate the computation of the test statistic, we approximate the probability distribution function (p.d.f.) $`f_0`$ of the null hypothesis by simulated observations. The statistic “energy” can also be used to test whether two samples belong to the same parent distribution.

In Section 2 we define the test and study its properties. Section 3 contains a comparison of the power of the energy test with that of other popular tests in one dimension. In Section 4 we apply the energy test to two-dimensional problems. We discuss possible modifications and extensions of the test in Section 5 and conclude with a summary in Section 6.

## 2 The test statistic

### 2.1 The energy function

We define a quantity $`\varphi `$, the *energy*, which measures the difference between two p.d.f’s $`f_0(𝐱)`$ and $`f(𝐱)`$, $`𝐱`$ $`^n`$, by

$$\varphi =\frac{1}{2}\left[f(𝐱)f_0(𝐱)\right]\left[f(𝐱^{})f_0(𝐱^{})\right]R(𝐱,𝐱^{})\mathrm{𝐝𝐱𝐝𝐱}^{}.$$
(1)
The integrals extend over the full variate space. The weight function $`R`$ is a monotonically decreasing function of the Euklidian distance $`|𝐱𝐱^{}|`$. Relation (1) with $`R=1/|𝐱𝐱^{}|`$ is the electrostatic energy of two charge distributions of opposite sign which is minimum if the charges neutralize each other. Setting the function $`R=\delta (𝐱𝐱^{})`$, where $`\delta (𝐱𝐱^{})`$ is the Dirac Delta function, $`\varphi `$ reduces to the integrated quadratic difference of the two p.d.f.s

$$\varphi _\delta =\frac{1}{2}\left[f(𝐱)f_0(𝐱)\right]^2\mathrm{𝐝𝐱}.$$

Since we want to generalize (1) in such a way that we can apply it to a comparison of a sample with a distribution or of two samples with each other, $`\varphi _\delta `$ is not suitable. Functions like Gaussians which correlate different locations have to be used.

The test which we will introduce is based on the fact that for fixed $`f_0(𝐱)`$, $`\varphi `$ is minimum for the null hypothesis $`f(𝐱)f_0(𝐱)`$ for all $`𝐱`$ $`^n`$. We sketch a prove of this assertion in the Appendix.

Expanding (1)

$$\varphi =\frac{1}{2}[f(𝐱)f(𝐱^{})+f_0(𝐱)f_0(𝐱^{})2f(𝐱)f_0(𝐱^{})R(|𝐱𝐱^{}|)]\mathrm{𝐝𝐱𝐝𝐱}^{},$$

we obtain three terms which have the form of expectation values of $`R`$. They can be estimated from two samples $`𝐗_1,𝐗_2,\mathrm{},𝐗_N`$ and $`𝐘_1,𝐘_2,\mathrm{},𝐘_M`$ drawn from $`f`$ and $`f_0`$, respectively. Since in the first and the second term a product of identical distributions occurs, there it is not necessary to draw two different samples of the same p.d.f.. The expectation value of $`R`$ in these expressions is given by the mean value of $`R`$ computed from all combinations of the sample observations

$$\varphi _{NM}=\frac{1}{N(N1)}\underset{j>i}{}R(|𝐱_i𝐱_j|)+\frac{1}{M(M1)}\underset{j>i}{}R(|𝐲_i𝐲_j|)\frac{1}{NM}\underset{i,j}{}R(|𝐱_i𝐲_j|).$$

Throughout this article unless specified differently, all sums run from $`1`$ to the maximum value of the index. For the energy statistic $`\varphi _N`$ of a sample $`𝐗_1,𝐗_2,\mathrm{},𝐗_N`$ relative to the p.d.f. $`f_0`$ is given by:

$`\varphi _N`$ $`={\displaystyle \frac{1}{N(N1)}}{\displaystyle \underset{j>i}{}}R(|𝐱_i𝐱_j|)+{\displaystyle f_0(𝐱)f_0(𝐱^{})R(|𝐱𝐱^{}|)\mathrm{𝐝𝐱𝐝𝐱}^{}}+`$
$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{i}{}}{\displaystyle f_0(𝐱^{})R(|𝐱_i𝐱^{}|)\mathrm{𝐝𝐱}^{}}.`$
To test whether a statistical sample $`𝐗_1,𝐗_2,\mathrm{},𝐗_N`$ of size $`N`$ is compatible with the null hypothesis $`H_0`$, we could in principle use the test statistic $`\varphi _N`$ of the sample with respect to the associated p.d.f. $`f_0`$ according to $`\varphi _N`$ but since the evaluation of $`\varphi _N`$ usually requires a sum over difficult integrals, we prefer to represent $`f_0`$ by a sample $`𝐘_1,𝐘_2,\mathrm{},𝐘_M`$, usually generated through a Monte Carlo simulation. We drop in $`\varphi _{NM}`$ the term depending only on $`f_0`$ which is independent from the sample observations and obtain:

$$\varphi _{NM}=\frac{1}{N^2}\underset{j>i}{}R(|𝐱_i𝐱_j|)\frac{1}{NM}\underset{i,j}{}R(|𝐱_i𝐲_j|)$$
(2)
Furthermore, we have replaced the denominator of the first term $`1/N(N1)`$ by $`1/N^2`$ which has superior small sample properties (see Appendix). Statistical fluctuations of the simulation are negligible if $`M`$ is large compared to $`N`$, typically $`M10N`$.

### 2.2 The weight function

We have investigated three different types of weight functions, power laws, a logarithmic dependence and Gaussians.

$`R_{pow}(r)`$ $`=\{\begin{array}{c}\frac{1}{r^\kappa }\text{for }r>d_{\mathrm{min}}\hfill \\ \frac{1}{d_{\mathrm{min}}^\kappa }\text{for }rd_{\mathrm{min}}\hfill \end{array}`$ (5)
$`R_{\mathrm{log}}(r)`$ $`=\{\begin{array}{c}\mathrm{ln}r\text{for }r>d_{\mathrm{min}}\hfill \\ \mathrm{ln}d_{\mathrm{min}}\text{for }rd_{\mathrm{min}}\hfill \end{array}`$ (8)
$`R_G(r)`$ $`=\mathrm{exp}\left(r^2/(2s^2)\right)`$ (9)
The first type is motivated by the analogy to electrostatics, the second is long range and the third emphasizes a limited range for the correlation between different observations. The power $`\kappa `$ of the denominator in (5) and the parameter $`s`$ in (9) may be chosen differently for different dimensions of the sample space and different applications. For slowly varying $`f_0`$ a small value of $`\kappa `$ around $`0.1`$ is recommended. For short range variations the test quantity with larger values around $`0.3`$ is more sensitive.

Also the logarithmic function $`R_{\mathrm{log}}`$ is well adapted to slowly varying $`f_0`$. The corresponding test quantiles are invariant under the transformation $`rar`$. For more rapidly varying $`f_0`$, the Gaussian weight function $`R_G`$ is recommended. It permits to adjust the range of the smearing to the shape of $`f_0`$.

The inverse power law and the logarithm exhibit a pole at $`r`$ equal to zero which could produce infinities in the double discrete version of the energy function $`\varphi _{NM}`$. Very small distances, however, should not be weighted too strongly since deviations from $`f_0`$ with sharp peaks are not expected and usually inhibited by the finite experimental resolution. We eliminate the poles by introducing a lower cutoff $`d_{\mathrm{min}}`$ for the distances $`r`$. Distances less than $`d_{\mathrm{min}}`$ are replaced by $`d_{\mathrm{min}}`$. The value of this parameter is not critical at all, it should be of the order of the average distance of the simulation points in the region where $`f_0`$ is maximum.

### 2.3 The distribution of the test statistic

The distribution of $`\varphi `$ depends on the distribution function $`f_0`$ and on the weight function $`R`$. Figure 1 shows the distributions for uniform $`f_0`$ with Gaussian and logarithmic weight functions. The distributions are well described by a generalized extreme value distribution

$$f(x)=\frac{1}{\sigma }\left(1+\xi \frac{x\mu }{\sigma }\right)^{1/\xi 1}\mathrm{exp}\left\{\left(1+\xi \frac{x\mu }{\sigma }\right)^{1/\xi }\right\}$$

depending on three parameters, a scale parameter $`\sigma `$, a location parameter $`\mu `$, and a shape parameter $`\xi `$. Rather than computing these parameters from the moments of the specific $`\varphi `$ distributions, we propose to generate the distribution of the test statistic and the quantiles by a Monte Carlo simulation. As a consequence of the dramatic increase of computing power during the last decade, it has become possible to perform the calculations on a simple PC within minutes. There is no need to publish tables of percentage points, to distinguish between simple and composite hypotheses and also censoring can be incorporated into the simulation.

### 2.4 Relation to the Pearson $`\chi ^2`$ test

Let us assume that we have in a certain $`\chi ^2`$-bin $`n_i`$ experimental observations and $`m_i`$ Monte Carlo observations and $`\beta =M/N`$ $`1`$. Then the $`\chi ^2`$ contribution of that bin is

$`\chi _i^2`$ $`={\displaystyle \frac{(n_im_i/\beta )^2}{m_i/\beta }}`$
$`={\displaystyle \frac{n_i^22n_im_i/\beta +m_i^2/\beta ^2}{m_i/\beta }}`$
For a goodness-of-fit test, an additive constant is irrelevant. Thus we can drop $`m_i^2`$ in the last relation. If in addition the theoretical distribution is uniform and the bins have constant size, we can also ignore the denominator.

$$\chi _i^2=\frac{n_i^2}{2}\frac{n_im_i}{\beta }$$

Up to a constant factor this last relation corresponds nearly to the energy defined in (2) for a weight function which is constant inside the bin and zero outside:

$$\varphi _{NM}=\frac{n_i(n_i1)}{2N^2}\frac{n_im_i}{NM}\frac{1}{N^2}\chi _i^2$$

The $`\chi ^2`$ test applied to a uniform distribution $`f_0`$ is equivalent to an energy test with a box shaped weight function and fixed box locations.

Replacing the box function by the Gaussian weight function, the sharp cut at the bin boundaries and the arbitrary location of the boundary for fixed binning are avoided but the idea underlying the $`\chi ^2`$ test is retained.

### 2.5 Optimization of the test parameters in a simple case

To study the dependence of the power of the test on the choice of the weight function and its parameters, we have chosen for $`H_0`$ a uniform, univariate p.d.f.$`f_0`$ restricted to the unit interval $`[0,1]`$. We determined the rejection power with respect to contaminations of $`f_0`$ with a linear and two different Gaussian distributions which represent a wide and a more local distortion.

$`f_0`$ $`=1`$
$`f_1(X)`$ $`=2X`$ (10)
$`f_2(X)`$ $`=c_2\mathrm{exp}\left(64(X0.5)^2\right)`$ (11)
$`f_3(X)`$ $`=c_3\mathrm{exp}\left(256(X0.5)^2\right)`$ (12)
We required 5% significance level and computed the rejection power which is equal to one minus the probability for an error of the second kind.

As a reference, we also computed the power of a $`\chi ^2`$ test with bins of fixed width. The number of bins $`B2N^{2/5}`$ was chosen according to the prescription proposed in .

The cut-off parameter $`d_{\mathrm{min}}`$ for $`R_\kappa `$ (power law) and $`R_l`$ (logarithmic) was set equal to $`1/(4N)`$.

In Figure 2 we show the results for samples of 100 observations for 70 % contamination with $`f_1`$, 30 % contamination with $`f_2,`$ and 20 % contamination with $`f_3`$. Five different values of the Gaussian width parameter $`s`$, four different power laws and the logarithmic weight function have been studied.

As expected, the linear distribution is best discriminated by slowly varying weight functions like the logarithm, low power laws and the wide Gaussians.

The contamination with the narrow Gaussian $`f_3`$ is better recognized by the narrow weight functions ($`s=1/8,\kappa =0.3`$). Here the two wide Gaussian weight functions fail completely.

Figure 3 illustrates the dependence of the rejection power on the sample size for the three different distortions (10), (11), (12) of the uniform distribution. The amount of contamination was reduced with increasing sample size.

We applied the energy tests with power law $`\kappa =0.3`$, the logarithmic and two Gaussian weight functions with fixed width $`s=1/8`$ (Gfix) and variable width (Gvar). The variable width was chosen such that the full width at half maximum is equal to the $`\chi ^2`$ bin width, chosen according to the $`2N^{2/5}`$ law. This allows for a fair comparison between the two methods.

As expected the linear distribution is best discriminated by the energy test with logarithmic weight function. The power of the $`\chi ^2`$ test is considerably worse. The energy test with variable Gaussian weight function follows the trend of the $`\chi ^2`$ test but performs better in 14 out of the 15 cases.

Comparing the samples with sizes 50 and 100, respectively, we realize one of the caveats of the $`\chi ^2`$ test: For the sample of size 50 there are 9 bins. The central bin coincides favorably with the location of the distortion peak of the background sample and consequently leads to a high rejection power. For the sample of size 100, however, there are 12 bins, thus two bins share the narrow peak and the power is reduced. The Gaussian energy test is insensitive to the location of the distortion.

## 3 Comparison with alternative univariate tests

We have investigated the following goodness-of-fit tests, some of which have been designed for special applications and can be found in : $`\chi ^2`$ test, Kolmogorov test, Kuiper test, Anderson-Darling test, Watson test, Neyman smooth test, Region test, Energy tests with logarithmic and Gaussian weight functions. The region test has been developed to detect localized bumps by partitioning the variate range of the order statistic into three regions with expected probabilities $`p_i,p_j,1p_ip_j`$. For $`N_i`$ and $`N_j`$ observations observed in the first two regions, the test statistic is

$$R=\underset{i,j}{sup}\{(N_iNp_i)^2+(N_jNp_j)^2+[N_i+N_jN(p_i+p_j)]^2$$

where $`N_i`$ is the number of observations with $`xx_i`$.

Samples contaminated by different background sources were tested against $`H_0`$, corresponding to the uniform distribution. The power of each test for a 5 % significance limit was evaluated from $`10^5`$ Monte Carlo simulations containing 100 observations each.

Samples with background contamination by a linear distribution, Gaussians and the following distributions displayed in Figure 4 were generated by

$`b_A(x)`$ $`(1x)^3`$
$`b_B(x)`$ $`\{\begin{array}{c}x^2\\ (1x)^2\end{array}\begin{array}{c}\text{for }0x<0.5\\ \text{for }0.5x<1\end{array}`$
$`b_C(x)`$ $`(x0.5)^2`$
.

The power of the different tests is presented in Figure 5. As expected, none of the tests is optimum for all kind of distortions. The energy tests are quite powerful independent of the background function.

## 4 The energy test for two or more variates

While numerous goodness-of-fit tests exist for univariate distributions in higher dimensions only tests of the $`\chi ^2`$ type have become popular.

The extension of binning-free tests based on linear rank statistics to the multivariate case  is difficult because no obvious ordering scheme of the observations exist. Sequential ordering in the different variates depends on the sequence in which the variates are chosen. In addition to this logical difficulty there is a more fundamental problem: The re-shuffling of the observations destroys the natural metric used to display the data. The situation in goodness-of-fit problems is very similar to that in pattern recognition. The possibility to detect distortions of the expected distribution by background or resolution effects depends on an appropriate selection of the random variables. The natural choice is usually given by the experiment. By mapping a multivariate space onto a space defined by the ordering recipe much of the information contained in the original distribution may be lost.

### 4.1 Test for bivariate Gaussian distribution

It is rather common to deal with samples of observations which are drawn from a multivariate normality. Multivariate goodness-of-fit tests have been developed mostly for multinormality. An overview can be found in . An extension of the moment tests introduced by K. Pearson to two variates is Mardia’s test . Skewness $`\beta _1`$ and kurtosis $`\beta _2`$ are compared to the nominal values for the normal distribution, $`\beta _1=0`$ and $`\beta _2=3`$ respectively. We compared the logarithmic and the Gaussian energy tests to Mardia’s tests and to the two-dimensional Neyman smooth test .

We chose a standard Gaussian

$$f_0=\frac{1}{2\pi }\mathrm{exp}\left((x^2+y^2)/2\right)$$

to represent the null hypothesis.

Background samples added to $`f_0`$ are shown in Figure 6. In addition, a uniform background in the region ($`0<X,Y<1`$) and variations of the Gaussian background were investigated. All Gaussians have the same width in both coordinates $`X`$ and $`Y`$, but differ in the value of the width and the correlation coefficient.

Figure 7 displays the results. We were astonished how well the energy test competes with alternatives especially designed to test normality. We attribute the excellent performance of the energy test to the fact that it is sensitive to all deviations of two distributions whereas the Mardia tests are based only on two specific moments.

### 4.2 Example from particle physics

Figure 8 is a scatter plot of data obtained in the particle physics experiment HERA-B at DESY. The positions of tracks from $`\psi `$ decays are plotted against the momentum of the $`\psi `$ meson. The 20 square entries represent the experimental data, the black dots correspond to a Monte Carlo simulation.

We have computed the energy statistic $`\varphi `$ with the logarithmic weighting function for the experimental data relative to the Monte Carlo prediction. Separately, the energy was determined for 100 independent Monte Carlo samples of 20 events each to estimate the distribution of $`\varphi `$ under $`H_0`$. The experimental value is larger than all Monte Carlo values (see Figure 8), indicating that the data do not follow the prediction.

## 5 Extensions and variations of the energy test

### 5.1 Normalization to probability density

The Gaussian energy test has proven to be similar but superior to the $`\chi ^2`$ test for the uniform distribution. For non-uniform p.d.f.s the similarity can be retained by effectively dividing the elements of the sums in (2) by the p.d.f. $`f_0`$.

$$R(|𝐱_i𝐱_j|)R^{}(|𝐱_i𝐱_j|)=\frac{R(|𝐱_i𝐱_j|)}{\sqrt{f_0(𝐱_i)f_0(𝐱_j)}}$$

If the density is only given by a Monte Carlo simulation, the density can be estimated from the volume of a multidimensional sphere in the variate space containing a fixed number of Monte Carlo observations. In the example of Figure 8 the density at $`𝐱`$ could be estimated from the area $`A(𝐱)`$ of the smallest circle centered at $`𝐱`$ containing 10 Monte Carlo points: $`f_0(𝐱)10/A(𝐱)`$.

### 5.2 Width of the weighting function $`R`$ and metric

Some statisticians recommend to use equal probability bins for the $`\chi ^2`$ test in one dimension. An equivalent procedure for the Gaussian energy test would be to adjust the width of the Gaussian weighting function to the p.d.f.

$`R_s(r_{ij})`$ $`R_s^{}(r_{ij})=\mathrm{exp}\left(r^2/(2s_is_j)\right)`$
$`s_i`$ $`1/f_0(𝐱_i)`$
Such an adjustment can be applied independent of the dimension of the variate space while an equal probability binning is rather difficult in two or higher dimensions.

The resolution of $`R`$ should be adjusted to the expected deviations. The constants $`s`$ can be chosen differently for the different coordinates. Alternatively, a variate transformation may be performed before the test is applied.

The multi-dimensional logarithmic energy test requires a sensible choice of the metric. The distances of the observations should be about equal numerically in all directions of the variate space. This can be achieved by normalizing the variates to the square root of their variance.

### 5.3 Clustering of observations

For very large samples, the computation of $`\varphi `$ may become excessively long even with powerful computers. A possible solution to this problem is to combine observations to clusters. The cluster weight is equal to the sum of the observations included in the cluster and its location is their center of gravity. The details of the clustering are not important. The simplest method consists in combining points in bins formed by a simple grid. Another possibility could be the following: Observations and clusters (once clusters have been formed) are chosen randomly and combined with all points and clusters which are within a fixed maximum distance. The process is terminated when the number of clusters is below an acceptable limit.

### 5.4 Are two samples drawn from the same distribution?

A frequently occurring problem is to check whether two random samples belong to the same unknown parent distribution. In the framework of the “energy” concept we can use bootstrap or permutation methods to deduce a reference energy distribution. The energy computed according to (2) is then compared to this distribution. Quantiles can be used to measure the compatibility of the samples.

## 6 Summary and conclusions

Energy tests represent a powerful alternative to conventional tests. To our knowledge the energy test is the first multi-variate binning-free test which is independent of a subjective ordering system. With a logarithmic weighting function it is very sensitive to long range distortions of the distribution to be tested, a situation which is frequent in physics applications. With a Gaussian weighting function the energy test has similar properties as tests of the $`\chi ^2`$ type, but avoids arbitrary bin boundaries. It is more powerful than the $`\chi ^2`$ test in almost all cases which we have considered. It can cope with arbitrarily small sample sizes. It is astonishing how well the energy test can compete with the Mardia test for testing normality in two dimensions.

The necessity to determine the distribution function of the test statistic for the specific application is not a serious restriction with today’s computing power.

The energy statistic can be used to test whether different samples belong to the same parent distribution. A corresponding study will be presented in a forthcoming article.

## 7 Appendix

### 7.1 Continuous case

We conjecture that $`f(𝐱)f_0(𝐱)`$ corresponds to a minimum of $`\varphi `$ in (1). Substituting $`g(𝐱)=f(𝐱)f_0(𝐱)`$ we obtain

$$\varphi =\frac{1}{2}g(𝐱)g(𝐱^{})R(|𝐱𝐱^{}|)\mathrm{𝐝𝐱𝐝𝐱}^{}0$$
(13)
where $`g`$ fulfils the condition $`g(𝐱)\mathrm{𝐝𝐱}=0`$. We conjecture that $`\varphi `$ is minimum for $`g(𝐱)0`$ for weight functions $`R`$ which decrease continuously with $`|𝐱𝐱^{}|`$.

In the following, we sketch a prove of this conjecture without full mathematical rigor and generality.

The property (13) is invariant under a linear transformation

$$RR^{}=\alpha R+const$$

where $`\alpha `$ is a positive scaling factor. Assuming that $`R`$ is finite, we can restrict its values to the range $`1R1`$ with $`R(0)=1`$. We approximate the function $`g`$ by a histogram with function value $`g_i`$ in bin $`i`$. The size of each bin is $`\mathrm{\Delta }^n`$.

$$\varphi _h=\mathrm{\Delta }^{2n}\underset{i,j}{}g_ig_jR_{ij}$$

with self-explaining notation. The weights can be substituted by cosines, $`\mathrm{cos}\theta _{ij}=R_{ij}`$ with $`\mathrm{cos}\theta _{ii}=1`$ in our approximation. The sum then can be written as a sum of vectors $`𝐠_i`$ in $`^n`$:

$`\varphi _h`$ $`=\mathrm{\Delta }^{2n}{\displaystyle \underset{i,j}{}}g_ig_j\mathrm{cos}\theta _{ij}`$
$`\varphi _h`$ $`=\mathrm{\Delta }^{2n}\left({\displaystyle \underset{i}{}}𝐠_i\right)^20`$
The minimum of $`\varphi _h`$ is realized only if all $`𝐠_i`$ are equal to zero.

### 7.2 Discrete case

When we approximate both densities $`f`$ and $`f_0`$ by distributions of $`N`$ points each, with weights equal $`1/N`$, we require that the energy be minimum if the two samples coincide, $`𝐱_j=𝐲_j`$ $`j`$. To fulfil this condition, we have to replace the factors $`1/(N1)`$ and $`1/(M1)`$ by $`1/N`$ in $`\varphi _{NM}`$.

$`\varphi _{NN}`$ $`={\displaystyle \frac{1}{N^2}}{\displaystyle \underset{j>i}{}}R(|𝐱_i𝐱_j|)+{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{j>i}{}}R(|𝐲_i𝐲_j|)+`$
$`{\displaystyle \frac{1}{N^2}}{\displaystyle \underset{i,j}{}}R(|𝐱_i𝐲_j|)`$
To demonstrate the minimum condition which leads to $`\varphi _{NN}=0`$, we apply an infinitesimal shift to one observation. For $`𝐱_j=𝐲_j`$ for $`jk`$ and $`𝐱_k𝐲_k=\delta 𝐱_k`$, only the pair $`𝐱_k,𝐲_k`$ contributes to the energy, all other terms cancel.

$$\varphi _{NN}=\frac{1}{N^2}[R(|\delta 𝐱_k|R(0)]$$

Since $`R`$ *decreases* with its argument, $`R(0)R(|\delta 𝐱_k|)>0`$, we have found a local *minimum* of the energy.

For the special choice where $`R`$ is the distance function, i. e. *increasing* with the distance, it has been demonstrated by  that $`\varphi `$

$$\varphi =\frac{1}{N^2}\underset{j>i}{}\left[|𝐱_i𝐱_j|+|𝐲_i𝐲_j|\right]\frac{1}{N^2}\underset{i}{}\underset{j}{}|𝐱_i𝐲_j|0$$

is *maximum* for $`𝐱_j=𝐲_j`$ for all $`j`$. Since $`\varphi `$ is constant for constant $`R`$, it is plausible that for discrete distributions $`\varphi `$ is minimum for $`𝐱_j=𝐲_j`$ for all decreasing functions of the distance, and maximum for all increasing functions of the distance, however, we did not succeed in finding a general prove of this conjecture.

Acknowledgments. The authors wish to express their gratitude to Prof. R.-D. Reiss for his constructive comments on earlier drafts of this paper, which led to improvements in the presentation.


format: text

text:
# 1 Finite Temperature Pion Scattering

## 1 Finite Temperature Pion Scattering

There is a growing interest in the Quark-Gluon-Plasma produced in Relativistic Heavy Ion Collisions (RHIC) and its evolution into hadronic matter. After the chiral phase transition, this plasma cools down rapidly, mainly into a pion gas, which requires a description of QCD at temperatures below the chiral phase transition. With this aim we turn to Chiral Perturbation Theory (ChPT)  which is the low-energy effective theory of QCD. It follows from the identification of pions as the Nambu-Goldstone bosons associated to the QCD $`SU(2)_L\times SU(2)_RSU(2)_{L+R}`$ spontaneous chiral symmetry breaking and has been already applied to describe static properties of the low-$`T`$ pion gas . Within ChPT, we will analyze the thermal $`\pi \pi `$ scattering, which is relevant for the gas description in terms of transport equations, for the evolution of the expanding fireball or for estimating the freeze out time (when the scattering mean time is larger than the collective expansion time scale). Pion scattering in the $`\sigma `$ channel has also been related to chiral symmetry restoration . The thermal amplitude in the vector channel has been studied within a modified Nambu-Jona Lasinio (NJL) model  explicitly including the $`\rho \pi \pi `$ vertex and analyzing the complications due to a non-physical constituent quark-antiquark threshold near the $`\rho `$ mass. Within ChPT only the scattering lengths are available at finite $`T`$  (for a NJL calculation see ). In this work we calculate the whole $`\pi \pi `$ elastic scattering amplitude at finite temperature to one loop in ChPT, thus valid at energies and temperatures well below $`\mathrm{\Lambda }_\chi `$ 1.2 GeV, both considered as $`𝒪(p)`$ in the chiral expansion.

A further motivation for our analysis is one of the most striking observations of RHIC: the modification of the dilepton spectrum with respect to the vacuum hadronic emission . This is usually interpreted as an in-medium modification of the $`\rho `$ mass and width . Here we will obtain a completely model independent and systematic description of the low energy tail of this resonance, in terms of $`\pi \pi `$ scattering. Let us remark that we are taking into account the details of the pion interactions, not only their distribution function  whose effect comes from the imaginary part of the amplitude, as we will see below. The description of the full resonance is beyond ChPT. Nevertheless, we show here that our amplitude satisfies a thermal generalization of perturbative unitarity, which is the starting point  of the successful nonperturbative unitarization methods used to study resonances like the $`\rho `$ and $`\sigma `$ at $`T=0`$ .

## 2 The $`T>0`$ Pion Scattering Amplitude

### 2.1 Definition and interpretation of the thermal amplitude

Let us start by discussing briefly what we mean by thermal scattering amplitude in what follows and why it contains interesting physical information. For our purposes, it will be enough to assume that the system is in thermal equilibrium and to neglect finite baryon density effects (as in an ideal RHIC in the central rapidity region).

In Thermal Field Theory, the word ”amplitude” is often used to denote a thermal Green function rather than an $`S`$-matrix element, whose definition in a thermal bath is rather subtle. At $`T=0`$ the $`S`$-matrix element is related to a time-ordered product through the LSZ reduction formula. Such a time-ordered product can be calculated at $`T0`$, where the time components are ordered along a complex contour $`C`$ . In the imaginary-time formalism (ITF) $`C=[0,i\beta ]`$ so that the time arguments become Euclidean $`\tau =itC`$, with $`\beta =T^1`$ and $`t`$ being the Minkowski time. The Feynman rules in the ITF are simply those at $`T=0`$ but replacing all the zeroth components of momenta by a discrete set of frequencies $`q^0i\omega _n=2\pi inT`$ and the loop integrals by Matsubara sums, i.e., $`(2\pi )^1𝑑q^0iT_{n=\mathrm{}}^{\mathrm{}}`$. Since the vertices remain unchanged the set of diagrams is exactly the same as for $`T=0`$. The ITF is appropriate to calculate thermodynamical quantities, but it can also be used to analyze dynamical processes with real energies. In fact, carrying out the Matsubara sums first, the ITF Green functions yield analytic functions of $`i\omega _n`$ everywhere off the real axis . Therefore, they can be continued analytically from discrete external complex to continuous real energies $`E`$, usually with the prescription $`i\omega _nE+iϵ`$. It is worth mentioning that the thermal Green functions thus obtained correspond to the retarded functions in the real-time formalism (RTF)  for which $`C`$ includes the real axis. One could also calculate the time-ordered product directly in the RTF, where the external energies remain real, but the number of fields is doubled . In this sense, let us remark that in the RTF it has been shown  that the $`T0`$ retarded functions have a causal and analytic structure which allows for a description in terms of dispersion relations, unlike the $`T0`$ RTF time-ordered product. Hence, the retarded functions are suitable candidates to generalize important properties of the $`T=0`$ amplitude such as unitarity, as we will see below.

Therefore, back to our present case, consider the $`T=0`$ $`S`$-matrix element corresponding to $`\pi \pi `$ scattering, written through the LSZ formula in terms of the Green function of the time-ordered product of four pion fields. The thermal scattering amplitude we will consider is defined with respect to such $`T=0`$ amplitude by keeping the pion initial and final asymptotic states as those for $`T=0`$ but i) replacing the Green function by a thermal expectation value in the ITF and ii) performing the usual analytical continuation into continuous real energies. As discussed above, this corresponds to the retarded four-point function in the RTF and, as we will see, the thermal amplitude thus defined has the correct analytic structure, allowing for a suitable extension of the perturbative unitarity at $`T0`$. Let us remark that the same definition for the pion scattering thermal amplitude has been used in .

An alternative point of view would have been to consider the RTF thermal reaction rate $``$ for the inclusive reaction $`\pi \pi `$+Heat Bath $``$ $`\pi \pi `$ \+ anything  defined as the thermal average of the transition probability $`S^{}S`$, with the $`T=0`$ $`S`$-matrix element between an initial and final asymptotic state including the heat bath collective states. We remark that this is a different object than the thermal amplitude we consider here, thus providing information on different physical observables. Thus, while $``$ can be interpreted as an ”average of probability” when the free heat bath is taken also as an asymptotic state, the thermal amplitude we consider here can be viewed rather as (the retarded part of) an ”average of amplitude” of the probability that two initial free pions scatter in the interacting thermal bath and two final free pions leave it. In principle, there is no obvious relation between the two, except for the case of two-point functions, where the imaginary part of the retarded RTF Green function is proportional to the difference between the direct decay rate, given by $``$ and the inverse one, related to $``$ by detailed balance . A detailed calculation of the $`T0`$ reaction rate for the pion scattering using the rules given in  would be very interesting but it lies beyond the scope of this work.

### 2.2 The one loop ChPT calculation

The $`N_f=2`$ chiral lagrangian at $`𝒪(p^2)`$ is that of the $`O(4)`$ nonlinear sigma model:

$$_2=\frac{f^2}{2}\left(_\mu U^T^\mu U+2\chi ^TU\right)$$
(1)
where $`f`$ is the $`T=0`$ pion decay constant in the chiral limit, $`U(x)(\sqrt{1\pi ^2/f^2},\stackrel{}{\pi }/f)`$ is an $`O(4)`$ vector satisfying $`U^TU=1`$, with $`\pi _a`$ ($`a=1,2,3`$) the pion fields and $`\chi =(m^2,\stackrel{}{0})`$. Here, $`m`$ is the pion mass to lowest order (in the isospin limit). For simplicity, we do not show external sources, which are nevertheless needed to obtain the $`𝒪(p^4)`$ relations between the physical $`m_\pi `$, $`f_\pi `$ and $`m`$, $`f`$ . The $`𝒪(p^4)`$ lagrangian $`_4`$ includes seven independent terms multiplied by low-energy constants, called $`\overline{l}_i`$ after renormalization.

There is another subtlety to be borne in mind. The loss of Lorentz covariance imposed by the choice of the thermal bath rest frame implies that we have to separate explicitly the time and space-like components of any four-vector. In particular, for an elastic two-body collision with four-momenta $`k_1k_2k_3k_4`$, if we call $`𝐒=k_1+k_2`$, $`𝐓=k_1k_3`$ (not to be confused with the temperature) and $`𝐔=k_1k_4`$, the amplitude will depend differently on $`𝐒_\mathrm{𝟎}`$, $`|\stackrel{}{𝐒}|`$, $`𝐓_\mathrm{𝟎}`$, $`|\stackrel{}{𝐓}|`$, $`𝐔_\mathrm{𝟎}`$ and $`|\stackrel{}{𝐔}|`$. In particular, crossing symmetry on the Mandelstam variables $`s=𝐒^\mathrm{𝟐},t=𝐓^\mathrm{𝟐},u=𝐔^\mathrm{𝟐}`$ no longer holds. Let us recall that at $`T=0`$ any elastic $`\pi \pi `$ amplitude is related to that of $`\pi ^+\pi ^{}\pi ^0\pi ^0`$, called $`A(s,t,u)`$, by isospin and crossing transformations. Nevertheless, since the temperature does not modify the vertices, crossing symmetry still holds in terms of the $`S`$, $`T`$ and $`U`$ four-momenta, and we can still write any $`\pi \pi `$ thermal amplitude in terms of the thermal $`A(𝐒,𝐓,𝐔;\beta )`$:

$`A(𝐒,𝐓,𝐔;\beta )=A_2(𝐒,𝐓,𝐔)+A_4^{pol}(𝐒,𝐓,𝐔)+A_4^{tad}(𝐒,𝐓,𝐔;\beta )+A_4^{uni}(𝐒,𝐓,𝐔;\beta ),`$ (2)
where we have followed the usual $`T=0`$ notation: $`A_2`$ is the $`_2`$ tree level contribution and $`A_4^{pol}`$ contains the $`_4`$ tree level, depending on the $`\overline{l}_i`$, plus those polynomials coming from the renormalization of the loop integrals. Both $`A_2`$ and $`A_4^{pol}`$ are temperature independent. The $`A_4^{uni}`$ term represents those contributions from loops with two propagators proportional to the three independent integrals $`J_0,J_1`$ and $`J_2`$ defined in eq.(14) which yield the correct analytic structure at $`T=0`$ and will ensure perturbative unitarity in all channels. The terms collected in $`A_4^{tad}`$ are proportional to the tadpole integral $`F_\beta `$ in eq.(13) and include all tadpole diagrams as well as terms coming from loop integrals (see Appendix A).

Taking these remarks into account, the thermal amplitude is obtained by adding to the temperature independent part in  the following temperature dependent corrections:

$`\mathrm{\Delta }A_4^{tad}(𝐒,𝐓,𝐔;\beta )={\displaystyle \frac{\mathrm{\Delta }F_\beta }{f^4}}\{s+t+u+{\displaystyle \frac{5m^2}{2}}[{\displaystyle \frac{2𝐓_\mathrm{𝟎}}{|\stackrel{}{𝐓}|^2}}(k_1^0\stackrel{}{k}_2+k_2^0\stackrel{}{k}_1)\stackrel{}{𝐓}+2(k_1k_2)𝐓`$ (3)
$`+{\displaystyle \frac{2(\stackrel{}{k_1}\stackrel{}{𝐓})(\stackrel{}{k_2}\stackrel{}{𝐓})}{|\stackrel{}{𝐓}|^4}}[3({\displaystyle \frac{t}{2}}+𝐓_{\mathrm{𝟎}}^{}{}_{}{}^{2})+|\stackrel{}{𝐓}|^2]+{\displaystyle \frac{t(\stackrel{}{k_1}\stackrel{}{k_2})}{|\stackrel{}{𝐓}|^2}}+(𝐓𝐔)]\}`$
$`\mathrm{\Delta }A_4^{uni}(𝐒,𝐓,𝐔;\beta )=`$ (4)
$`{\displaystyle \frac{1}{f^4}}\{(sm^2)[2k_{1}^{}{}_{}{}^{0}\mathrm{\Delta }J_1(𝐒)+{\displaystyle \frac{\mathrm{\Delta }J_0(𝐒)}{2}}(m^2+3s2k_3𝐒){\displaystyle \frac{\stackrel{}{k}_1\stackrel{}{𝐒}(s\mathrm{\Delta }J_0(𝐒)2𝐒_\mathrm{𝟎}\mathrm{\Delta }J_1(𝐒))}{|\stackrel{}{𝐒}|^2}}]`$
$`+[{\displaystyle \frac{m^2}{|\stackrel{}{𝐓}|^2}}[(\stackrel{}{k}_1\stackrel{}{k}_2)\stackrel{}{𝐓}[t\mathrm{\Delta }J_0(𝐓)2𝐓_\mathrm{𝟎}\mathrm{\Delta }J_1(𝐓)]]+m^4\mathrm{\Delta }J_0(𝐓)2(k_1^0k_2^0)m^2\mathrm{\Delta }J_1(𝐓)`$
$`4k_1^0k_2^0\mathrm{\Delta }J_2(𝐓){\displaystyle \frac{2}{|\stackrel{}{𝐓}|^2}}\left(k_1^0\stackrel{}{k}_2+k_2^0\stackrel{}{k}_1\right)\stackrel{}{𝐓}\left[t\mathrm{\Delta }J_1(𝐓)2𝐓_\mathrm{𝟎}\mathrm{\Delta }J_2(𝐓)\right]`$
$`{\displaystyle \frac{2(\stackrel{}{k_1}\stackrel{}{𝐓})(\stackrel{}{k_2}\stackrel{}{𝐓})}{|\stackrel{}{𝐓}|^4}}\left[\left(3𝐓_{\mathrm{𝟎}}^{}{}_{}{}^{2}|\stackrel{}{𝐓}|^2\right)\mathrm{\Delta }J_2(𝐓)3𝐓_\mathrm{𝟎}t\mathrm{\Delta }J_1(𝐓)+\left({\displaystyle \frac{3t^2}{4}}+m^2|\stackrel{}{𝐓}|^2\right)\mathrm{\Delta }J_0(𝐓)\right]`$
$`+{\displaystyle \frac{2\stackrel{}{k_1}\stackrel{}{k_2}}{|\stackrel{}{𝐓}|^2}}[t\mathrm{\Delta }J_2(𝐓)t𝐓_\mathrm{𝟎}\mathrm{\Delta }J_1(𝐓)+({\displaystyle \frac{t^2}{4}}+m^2|\stackrel{}{𝐓}|^2)\mathrm{\Delta }J_0(𝐓)]+(𝐓𝐔)]\}`$
Note that the above equations in Minkowski space-time are obtained from the Euclidean expressions by replacing $`k^0iE`$ for all Euclidean momenta $`k_0`$, with $`E`$ real, and analytically continuing the $`J_k`$ as showed in eq.(23). Note also that, for clarity, we have not simplified the $`k_i`$ four-momenta in terms of $`𝐒,𝐓`$ and $`𝐔`$. In addition, the $`T`$-corrections to $`m_\pi `$ and $`f_\pi `$  are included in $`\mathrm{\Delta }A_4^{tad}`$, consistently with our definition of the thermal amplitude. Recall that the $`T=0`$ physical $`f_\pi `$ and $`m_\pi `$ are :

$`m_\pi ^2=m^2(1{\displaystyle \frac{\overline{l}_3m^2}{32\pi ^2f^2}}),f_\pi =f(1+{\displaystyle \frac{m^2}{f^2}}{\displaystyle \frac{\overline{l}_4}{16\pi ^2}}).`$ (5)
We have performed three different consistency checks on the thermal amplitude: First, we have recovered the Lorentz invariant $`T0^+`$ limit of . Since all the temperature dependence is on the tadpoles and $`J_k`$ functions, we are thus checking their coefficients, since our calculation was carried out without Lorentz covariance. The other two checks concern partial waves and the c.o.m. frame, i.e, $`\stackrel{}{𝐒}=0`$ or pions at rest with the thermal bath. Customarily, the $`\pi \pi `$ scattering amplitude is projected into partial waves $`a_{IJ}`$ of definite isospin $`I`$ and angular momentum $`J`$, and the standard definitions  are still applicable at $`T0`$, but with crossing symmetry in terms of the $`𝐒,𝐓`$ and $`𝐔`$ four-vectors. This procedure is equivalent to a direct calculation of the isospin amplitudes without using the $`A`$ function, as done in , whose analytic expressions for the thermal scattering lengths we have reobtained, thus checking the real part of our partial waves with an independent calculation. The third check concerns the imaginary parts of the thermal partial waves in the different channels. This requires a generalization of the perturbative unitarity relation that we present next, which, in turn, will provide a natural physical interpretation of the imaginary part of our thermal amplitude in terms of direct (emission) and inverse (absorption) processes in the thermal bath.

### 2.3 Thermal unitarity

For simplicity, we drop the $`IJ`$ indices in the following. At $`T=0`$ unitarity constraints the partial waves, for $`s>4m_\pi ^2`$ and below other inelastic thresholds, to satisfy

$$\text{Im}a(s)=\sigma (s)|a(s)|^2,\sigma (s)=\sqrt{14m_\pi ^2/s},$$
(6)
whereas the ChPT series only satisfies unitarity perturbatively, i.e.,

$$\text{Im}a_2(s)=0,\text{Im}a_4(s)=\sigma (s)\left|a_2(s)\right|^2,\mathrm{}.$$
(7)
We now generalize the above perturbative relation for any one-loop calculation, including those within ChPT. Let us recall that for one-loop elastic amplitudes, any imaginary part should come from loop diagrams with two propagators and in particular from the one-loop integrals $`J_k`$ ($`k=0,1,2`$) given in eq.(14). As an illustration we will analyze $`J_0`$, since the other cases are completely analogous. We start then by performing the frequency sums in the standard fashion :

$`J_0(\omega ,|\stackrel{}{Q}|;T)={\displaystyle \frac{d^{D1}q}{(2\pi )^{D1}}\underset{s_1,s_2=\pm 1}{}\frac{s_1s_2}{4E_qE_{qQ}}\frac{1+n\left(s_1E_q\right)+n\left(s_2E_{qQ}\right)}{\omega s_1E_qs_2E_{qQ}}}`$ (8)
where $`E_k^2=|\stackrel{}{k}|^2+m^2`$, $`n(x)=\left(e^{\beta x}1\right)^1`$ is the Bose-Einstein distribution function and $`\omega iQ_0=2\pi kTi`$ with $`k`$ integer, is the external frequency. Note that eq.(8) provides the analytic continuation of $`J_0(\omega )`$ for $`\omega `$ off the real axis. Taking now $`\text{Im}J_0(E+iϵ)`$ with $`E`$ real (see our previous discussion on the thermal amplitude) the two $`s_1=s_2`$ terms cancel each other. The two cases $`s_1=s_2=\pm 1`$ require $`E\pm E_q\pm E_{qQ}=0`$ which imply $`\pm E\sqrt{|\stackrel{}{Q}|^2+4m^2}>0`$, respectively. Summarizing:

$`\text{Im}J_0(E+iϵ,|\stackrel{}{Q}|;T)`$ $`=`$ $`\pi \text{sgn}(E)\theta (s4m^2)`$
$`\times `$ $`{\displaystyle \frac{d^3\stackrel{}{q}}{(2\pi )^3}\frac{1}{4E_qE_{qQ}}\left[1+n\left(E_q\right)+n\left(E_{qQ}\right)\right]\delta \left(|E|E_qE_{qQ}\right)}`$
where $`s=E^2|\stackrel{}{Q}|^2`$ and $`E`$ real. Hence $`J_0(s)`$ has a cut on the real axis starting at $`s=4m^2`$, as in the $`T=0`$ case. Our analysis is similar to that of a decaying particle in medium , where the integral in eq.(8) is the thermal self-energy. In fact, the interpretation of eq.(2.3) follows from the observation that $`1+n(E_1)+n(E_2)=[1+n(E_1)][1+n(E_2)]n(E_1)n(E_2)`$. Therefore, in our thermal amplitude and for positive energy, the in-medium pions increase the available phase space by $`(1+n_1)(1+n_2)`$ (corresponding to the thermal enhancement of the two pion production) and decrease it by $`n_1n_2`$ (when the two initial pions collide with thermal pions and are absorbed) analogously to the stimulated emission and absorption processes discussed in . Alternatively, we could have calculated in the RTF, arriving to the same result for the retarded self-energy, unlike the RTF time-ordered product, which yields the sum of the emission and absorption terms rather than their difference .

Partial waves are defined in the c.o.m. frame, where we start studying the $`S`$-channel diagrams, setting $`\stackrel{}{Q}=0`$ in eq.(2.3). The generalization to one-loop ChPT is straightforward. Since $`_2`$ contains only two derivatives, one needs loop integrals with up to four momenta in the numerator, whose imaginary part can be related to $`J_0`$, $`J_1`$ and $`J_2`$ in eq.(14) (see Appendix A). The $`J_{1,2}`$ analysis is similar to that for $`J_0`$, and we find

$`\text{Im}J_k(E+iϵ,\stackrel{}{0};T)`$ $`=`$ $`{\displaystyle \frac{\text{sgn}(E)}{16\pi }}\theta (s4m^2)\sigma _T(E)\left({\displaystyle \frac{E}{2}}\right)^k,k=0,1,2`$ (10)
$`\text{where}\sigma _T(E)`$ $`=`$ $`\sigma (E^2)\left[1+{\displaystyle \frac{2}{\mathrm{exp}(\beta |E|/2)1}}\right]`$ (11)
is the thermal two-particle phase space in the c.o.m. frame . Note that in (10) we have continued analytically the $`J_k`$ functions to real energies as indicated in eq.(23) .

Concerning the $`T`$ and $`U`$ channels, in the c.o.m. frame and for $`s>4m_\pi ^2`$ we have to evaluate $`J_k(0,|\stackrel{}{𝐓}|)`$ (see Appendix A) which is real, according to eq.(2.3). Therefore, the same $`T=0`$ derivation of eq.(7) follows for $`T>0`$ by replacing $`\sigma \sigma _T`$, i.e.:

$$\text{Im}a_2(s)=0,\text{Im}a_4(s;\beta )=\sigma _T(𝐒_\mathrm{𝟎})\left|a_2(s)\right|^2,𝐒_\mathrm{𝟎}>2m_\pi $$
(12)
This generalizes the perturbative unitarity relation to $`T0`$, providing a neat physical interpretation of the imaginary part of the thermal amplitude in terms of emission and absorption of pion pairs in the thermal bath. Note that the one-loop thermal amplitude considered here, in terms of retarded propagators, has the same unitarity cut as the $`T=0`$ one. In fact, it can be shown that they have the same analytic structure in the complex plane , which gives support to the idea that the thermal retarded functions provide the natural extension of amplitudes as far as their causal and analytic properties are concerned. We remark that our perturbative result (12) is consistent within ChPT, where the external energies and $`T`$ remain small and so do the Bose-Einstein factors. For higher orders one could have intermediate states in the thermal bath with more than two pions, weighted by higher powers of $`n`$. We have explicitly checked that our final result for the amplitude satisfies eq.(12) for the $`IJ`$=00,11,$`20`$ channels.

Finally, let us note that regarding the $`\pi \pi `$ scattering in the $`11`$ channel as the exchange of a $`\rho `$ resonance , the pion scattering amplitude would be proportional to the $`\rho `$ thermal retarded self-energy, whose imaginary part to one loop is proportional to $`\sigma _T/\sigma `$ and gives the net decay rate of the $`\rho `$ in the thermal bath. Note that, from this point of view, the stimulated emission contribution, proportional to $`\left[1+n(|E|/2)\right]^2`$, would give the reaction rate discussed in  for the decay of a $`\rho `$ at rest into two pions. Hence, we expect to be able to extract physical information about the thermal properties of physical resonances arising in our $`\pi \pi `$ scattering amplitude, which we have calculated without introducing those resonances as explicit degrees of freedom. This is indeed the case, as shown recently in , where the thermal mass and width of the $`\rho `$ and $`\sigma `$ are read off from the complex poles (in the second Riemann sheet) of the thermal nonperturbative amplitude, obtained from the ChPT one by unitarization. At low $`T`$, the $`\rho `$ width grows with the $`\sigma _T/\sigma `$ factor, as our previous discussion suggested. This provides a relevant physical application of our present analysis, which is totally model independent, in contrast with the very recent results in , which maybe more direct to extract information on resonances but relies on the choice of unitarization method.

## 3 Discussion and conclusions

From our previous results, we have also calculated the scattering phase shifts, $`\delta _{IJ}`$ in different channels. This is a simple and easy way of representing the complex amplitude in terms of just one parameter, since elastic unitarity implies $`a=\mathrm{exp}(\delta )\mathrm{sin}(\delta )/\sigma _T`$. Since, for a given energy, the thermal phase space is fixed, the “strength of the interaction”, i.e. the modulus of the amplitude, is given only in terms of $`\mathrm{sin}\delta `$ This is similar to the $`T=0`$ case, and perturbatively corresponds to $`\delta _T\sigma _T\left(\sqrt{s}\right)\left[a_2(s)+\text{Re}a_4(s,\beta )\right]`$ . The results are plotted in Fig.1 up to 500 MeV, roughly where ChPT is reliable. Note that the absolute value of all the phase shifts is increased, but each channel keeps its attractive or repulsive nature (the sign of the phase shifts). The phase shift increase we find is basically due to the Bose enhancement factor in $`\sigma _T`$, the real part and the modulus of the perturbative amplitude changing very slightly with $`T`$ at low energies. In fact, in the $`IJ=20`$ channel, $`|\text{Re}a_4|`$ decreases with $`T`$ for low energies (see ). A strong enhancement of $`|a_{00}|^2`$ near threshold would be a signal of chiral symmetry restoration  but we do not observe such effect here. We also remark that not all the thermal effects in our one-loop calculation can be accounted for with a thermal redefinition of $`f_\pi `$  since in that case, with the decreasing $`f_\pi (T)`$ given in , we would see a larger increase of $`|a|^2`$ in all channels. Note also that our result can be alternatively interpreted as an enhanced effective pion interaction if we were describing them with the $`T=0`$ expressions, although bearing in mind that the main contribution comes from the phase space factor.

The $`IJ=11`$ channel is particularly interesting, since it is dominated by the $`\rho `$ resonance, whose thermal deformation has been related to the anomalous dilepton spectra in RHIC . By analyzing this channel at low energies, we are studying the tail of the resonance. The low-energy temperature enhancement of the phase shift can then be interpreted as a stronger thermal effect of the rho, although with pure ChPT calculations we cannot disentangle whether it is due to a larger width or to a smaller mass. Qualitatively, the chiral combination of parameters that dominate the phase shift in this channel is $`\overline{l}_1\overline{l}_2\mathrm{\Gamma }_\rho f^4/M_\rho ^5`$, following the resonance saturation hypothesis . In order to obtain larger phase shifts at $`T=0`$, we need somewhat larger $`\overline{l}_1\overline{l}_2`$ values. Thus, the low-energy thermal phase shifts would correspond at $`T=0`$ to a larger $`\mathrm{\Gamma }_\rho f^4/M_\rho ^5`$ ratio, although with the $`T=0`$ part only we cannot mimic accurately the low-energy thermal behavior. Therefore, our results are consistent with most models that expect a low-$`T`$ sizable thermal increase of the rho width due to the thermal increase of phase space, with an almost constant thermal $`M_\rho `$ . Our model independent perturbative result is confirmed by our discussion at the end of section 2.3 and by the unitarized results in .

In conclusion, we have calculated the $`\pi \pi `$ scattering to one loop in Chiral Perturbation Theory at finite temperature $`T`$. Our thermal amplitude has been calculated with $`T=0`$ asymptotic states and the thermal ITF Green function, whose analytic continuation corresponds to the retarded RTF Green function, which has a suitable causal and analytic structure. Despite the loss of Lorentz covariance of the thermal formalism, any $`\pi \pi `$ elastic amplitude can still be related to a single one by means of a generalized crossing symmetry in terms of the external four-momenta. We have found the thermal generalization of perturbative partial wave unitarity to one loop, which amounts to replacing the $`T=0`$ two-particle phase space with its thermal version. This accounts for emission and absorption of pion pairs inside the thermal bath, providing a physical interpretation of the imaginary part of the thermal amplitude and giving the expected analytic structure above threshold, consistently with causality and thermal unitarity. We have checked that our amplitudes satisfy this constraint. With these amplitudes, we have studied the low energy behavior of thermal $`\pi \pi `$ elastic scattering in a model independent formalism, finding an enhanced phase shift in all channels, as dictated by the thermal phase space, but with little change in the modulus of the partial waves. Physically, this means that one cannot mimic the low-$`T`$ behaviour of dynamical quantities just by scaling $`f_\pi f_\pi (T)`$.

One of the main physical observables that can be related to our thermal scattering amplitude are the masses and widths of the $`\rho `$ and $`\sigma `$ mesons . We have given here a low-$`T`$ qualitative description for the case of the $`\rho `$. Further quantitative and detailed information can be obtained if these results are extended to higher energies beyond pure ChPT. In particular, by means of unitarization methods  that generate dynamically light resonances like the $`\rho `$ and $`\sigma `$ mesons . Hence, our results are useful as a starting point for the description of relevant physical processes where chiral symmetry and unitarity play a crucial role, such as the properties of those resonances in the plasma formed after a Relativistic Heavy Ion Collision.

### Acknowledgments

The authors thank A. Dobado for useful comments and discussions. F. J. Ll-E thanks P. Bicudo, S. Cotanch, P. Maris, E. Ribeiro and A. Szczepaniak too. Work supported from the Spanish CICYT projects, FPA2000-0956,PB98-0782 and BFM2000-1326.

## Appendix A Loop integrals

The following are the integrals needed for our calculation in an arbitrary reference frame (not necessarily the c.o.m frame, where pions are at rest with the thermal bath). The signature of the metric in Euclidean space is $`()`$ and we will use the notation

$$_\beta \frac{d^Dq}{(2\pi )^D}\frac{1}{\beta }\underset{n=\mathrm{}}{\overset{n=+\mathrm{}}{}}\frac{d^{D1}q}{(2\pi )^{D1}}$$

where $`\beta =1/T`$, $`q_0=2\pi nT`$ and $`Q_0=2\pi kT`$ with $`k`$ integer denotes an external frequency. Also we denote $`\mathrm{\Delta }G(T)G(T)G(T=0)`$. First, we have the tadpole integral:

$$F_\beta _\beta \frac{d^Dq}{(2\pi )^D}\frac{1}{q^2m^2}.$$
(13)
Then, three integrals of pion loops with an external momentum flow:

$`J_k(Q_0,|\stackrel{}{Q}|;T)`$ $``$ $`{\displaystyle _\beta }{\displaystyle \frac{d^Dq}{(2\pi )^D}}{\displaystyle \frac{q_0^k}{(q^2m^2)((qQ)^2m^2)}},k=0,1,2`$ (14)
All other integrals with momenta in the numerator can be reduced to those four as

$`{\displaystyle _\beta }{\displaystyle \frac{d^Dq}{(2\pi )^D}}{\displaystyle \frac{q_i}{(q^2m^2)((qQ)^2m^2)}}`$ $`=`$ $`{\displaystyle \frac{Q_i}{|\stackrel{}{Q}|^2}}\left[Q_0J_1+{\displaystyle \frac{1}{2}}J_0Q^2\right].`$ (15)
$`{\displaystyle _\beta }{\displaystyle \frac{d^Dq}{(2\pi )^D}}{\displaystyle \frac{q_\mu Q^\mu }{(q^2m^2)((qQ)^2m^2)}}`$ $`=`$ $`{\displaystyle \frac{Q^2}{2}}J_0.`$ (16)
$`{\displaystyle _\beta }{\displaystyle \frac{d^Dq}{(2\pi )^D}}{\displaystyle \frac{q_0q_i}{(q^2m^2)((qQ)^2m^2)}}`$ $`=`$ $`{\displaystyle \frac{Q_i}{|\stackrel{}{Q}|^2}}\left[Q_0J_2+{\displaystyle \frac{Q^2}{2}}J_1+{\displaystyle \frac{Q_0}{2}}F_\beta \right].`$ (17)
$`{\displaystyle _\beta }{\displaystyle \frac{d^Dq}{(2\pi )^D}}{\displaystyle \frac{q_iq_j}{(q^2m^2)((qQ)^2m^2)}}`$ $`=`$ $`Q_iQ_jI_{2a}+g_{ij}I_{2b}`$ (18)
$`{\displaystyle _\beta }{\displaystyle \frac{d^Dq}{(2\pi )^D}}{\displaystyle \frac{q_0qq}{(q^2m^2)((qQ)^2m^2)}}`$ $`=`$ $`Q_0F_\beta +m^2J_1`$ (19)
$`{\displaystyle _\beta }{\displaystyle \frac{d^Dq}{(2\pi )^D}}{\displaystyle \frac{q_iqq}{(q^2m^2)((qQ)^2m^2)}}`$ $`=`$ $`Q_iF_\beta m^2{\displaystyle \frac{Q_i}{|\stackrel{}{Q}|^2}}\left(Q_0J_1+{\displaystyle \frac{Q^2}{2}}J_0\right).`$ (20)
$`{\displaystyle _\beta }{\displaystyle \frac{d^Dq}{(2\pi )^D}}{\displaystyle \frac{(qq)^2}{(q^2m^2)((qQ)^2m^2)}}=(Q^2+2m^2)F_\beta +m^4J_0.`$ (21)
where

$`I_{2a}(Q)`$ $`=`$ $`{\displaystyle \frac{1}{(D2)|\stackrel{}{Q}|^4}}[((D1)Q_0^2+|\stackrel{}{Q}|^2)J_2+(D1)Q_0Q^2J_1+({\displaystyle \frac{D1}{4}}Q^4+m^2|\stackrel{}{Q}|^2)J_0`$
$`+`$ $`(({\displaystyle \frac{Q^2}{2}}+Q_0^2)(D1)+|\stackrel{}{Q}|^2)F_\beta ]`$
$`I_{2b}(Q)`$ $`=`$ $`{\displaystyle \frac{1}{(D2)|\stackrel{}{Q}|^2}}\left[Q^2J_2+Q_0Q^2J_1+\left({\displaystyle \frac{Q^4}{4}}+m^2|\stackrel{}{Q}|^2\right)J_0{\displaystyle \frac{Q^2}{2}}F_\beta \right].`$ (22)
The above integrals have the appropriate $`T0`$, $`m0`$ limits and at $`T=0`$ they can be written just in terms of the tadpole $`F`$ and $`J_0`$ . Recall that all the UV divergences of the $`F_\beta `$ and $`J`$ integrals are contained in the $`T=0`$ part, since $`\mathrm{\Delta }F`$ and $`\mathrm{\Delta }J_k`$ always contain Bose-Einstein functions vanishing exponentially for large momentum. Therefore, using dimensional regularization, we always separate the $`T=0`$ part, regularized and renormalized as in  and set $`D=4`$ in $`\mathrm{\Delta }F_\beta `$ and $`\mathrm{\Delta }J_k`$.

The functions $`J_k`$ in eq.(14) can be analytically continued after performing the Matsubara sums. All the Euclidean discrete frequencies are then continued to $`iE`$ with $`E`$ real and continuous, i.e,

$`J_k(Q_0,|\stackrel{}{Q}|;T)`$ $``$ $`(i)^kJ_k(iE,|\stackrel{}{Q}|;T)`$ (23)
In the c.o.m. frame the imaginary part of the $`J`$ integrals, can be calculated exactly, e.g, eq.(10). When the c.o.m frame $`|\stackrel{}{𝐒}|=𝐓_\mathrm{𝟎}=𝐔_\mathrm{𝟎}=0`$, coincides with the thermal bath rest frame, only $`F_\beta `$, $`J_0(𝐒)`$, $`J_0(𝐓)`$, $`J_2(𝐓)`$ have to be evaluated (note that $`J_0(𝐔)`$, $`J_2(𝐔)`$ can be computed from the corresponding $`𝐓`$ integrals). Specifically

$`J_2(𝐒)`$ $`=`$ $`{\displaystyle \frac{𝐒_{\mathrm{𝟎}}^{}{}_{}{}^{2}}{4}}J_0(𝐒)+{\displaystyle \frac{F_\beta }{2}},J_1(𝐒)={\displaystyle \frac{𝐒_\mathrm{𝟎}}{2}}J_0(𝐒),J_1(𝐓)=0`$ (24)
The real part of $`\mathrm{\Delta }J_0`$ analytically continued as in eq.(23) in the $`S`$-channel reads

$$\text{Re}\mathrm{\Delta }J_0(𝐒)=\frac{1}{\pi ^2}PV_m^{\mathrm{}}𝑑E\frac{q(E)n(E)}{s4E^2}$$
(25)
where $`q(E)=\sqrt{E^2m^2}`$ $`n(E)=\left(e^{\beta E}1\right)^1`$ and $`PV`$ denotes Cauchy’s Principal Value. Note that the above expression is valid for $`\omega `$ real ($`s=\omega ^2`$) and the pole of the integrand is only present for $`s>4m^2`$, where $`\text{Im}J_00`$. Finally, in the $`T(U)`$ channels, we get

$`\mathrm{\Delta }J_k(|\stackrel{}{𝐓}|)`$ $`=`$ $`{\displaystyle \frac{1}{8\pi ^2}}{\displaystyle _0^{\mathrm{}}}dq{\displaystyle \frac{qn(E_q)E_q^{k1}}{|\stackrel{}{𝐓}|}}\mathrm{log}\left({\displaystyle \frac{2q+|\stackrel{}{𝐓}|}{2q|\stackrel{}{𝐓}|}}\right)^2,k=0,2`$
where $`E_q^2=q^2+m^2`$. The above result is valid for real $`|\stackrel{}{𝐓}|`$ (or $`s>4m_\pi ^2`$ for $`\pi \pi `$ scattering).