text:
# Optimal Moments for the Analysis of Peculiar Velocity Surveys II: Testing

## 1 INTRODUCTION

Although in principle the measurement of galaxy motions holds great promise as a probe of large–scale structure, in practice there are several obstacles that have prevented there from being robust conclusions made from the analyses of these measurements. First, these measurements are inherently noisy; errors in peculiar velocity determinations are typically of order $`10\%`$ of the redshift of a galaxy or cluster, which for all but nearby objects is comparable or larger than the velocity being measured. From the theoretical side, in order to relate the velocity field to the underlying matter density one must assume that the fields are linear. While this approximation is accurate on large scales, on smaller scales it generally fails due to infall into density concentrations. The difficulty in surmounting these obstacles is illustrated in the fact that attempts to compare different velocity field surveys have shown significant disagreements (Watkins & Feldman,, 1995; Hudson et al. , 1999).

Analyses of catalogs of peculiar velocity measurements have usually taken either of two main approaches. One is to average all the velocities to find the bulk flow, the velocity of the volume occupied by the survey relative to the Universal rest frame defined by the CMBR (e.g. Lauer & Postman (1994); Riess, Press & Kirshner (1995); Branchini, Plionis & Sciama (1996); Colless etal (2001); Aghanim & G rski (2001)). This method has the disadvantage that it discards most of the information contained in a survey and measures only three quantities, the components of the bulk flow vector. The second approach is to use all of the information contained in the survey in a likelihood analysis in order to obtain maximum likelihood estimates of the power spectrum parameters (e.g. Zaroubi et al.  (2001); Yang et al.  (2001); Zehavi & Dekel (2000); Freudling et al.  (1999)). This method potentially suffers from a bias due to small–scale, nonlinear contributions to velocities.

In a recent paper (Watkins et al.  (2002), hereafter Paper I), we introduced a new method for the analysis of peculiar velocity surveys that is a significant improvement over previous methods. In particular, our formalism allows us to separate information about large–scale flows from information about small scales, the latter which can then be discarded in the analysis. By applying specific criteria, we are able to retain the maximum information about large scales needed to place the strongest constraints, while removing the bias that small scale information can introduce into the results.

In paper I we reported on preliminary tests that suggested that our analysis method was effective. Here we present results from more extensive testing that demonstrates conclusively that our method works as advertised. In particular, we show that a likelihood analysis of simulated catalogs that uses all of the information contained in the peculiar velocities results in a bias in the estimation of the power spectrum parameters $`\mathrm{\Gamma }`$ and $`\beta `$, and that our method of analysis effectively removes this bias.

The paper is organized as followed: in sections 2 and 3 we briefly review the analysis formalism introduced in Paper I. In section 4 we discuss the N–body simulations we used and how we constructed synthetic catalogs. In section 5 we discuss results from our synthetic catalogs for both a full likelihood analysis of the catalog as well as our new analysis method. We also present additional evidence that our method is working effectively. In section 6 we conclude.

## 2 The Formalism

Our starting point is the usual statistical model for the line–of–sight peculiar velocities of galaxies. First, we assume that galaxies are tracers of a large–scale linear velocity field $`𝐯(𝐱)`$, a Gaussian random field completely described by the velocity power spectrum $`P_{(v)}(k)`$. With the assumption of linearity, the velocity power spectrum is proportional to the power spectrum of density fluctuations, with $`P_{(v)}(k)=H_o^2\mathrm{\Omega }^{1.2}P(k)`$.

For a set of $`N`$ galaxies with positions $`𝐫_𝐢`$, the observed line–of–sight velocity will be given by

$$v_i=𝐯(𝐫_i)\widehat{𝐫}_i+\delta _i,$$
(1)
the sum of the radial component of the velocity field with a noise term $`\delta _i`$ describing both observational error and any deviation from the linear velocity field due to local gravitational interactions. For simplicity we assume that $`\delta _i`$ is distributed as a Gaussian random variable with variance $`\sigma _i^2+\sigma _{}^2`$, where $`\sigma _i`$ is the observational error associated with that particular galaxy, and $`\sigma _{}`$ describes all other effects and is assumed to be the same for all the galaxies in the set.

With these assumptions, we can construct the probability distribution for the set of measured line–of–sight velocities given a power spectrum $`P(k)`$,

$$L(v_1,\mathrm{},v_N;P(k))=\sqrt{|R^1|}\mathrm{exp}\left(\underset{i,j=1}{\overset{N}{}}v_iR_{ij}^1v_j/2\right),$$
(2)
where $`R_{ij}=v_iv_j`$ is the covariance matrix, which in this case takes the form

$$R_{ij}=R_{ij}^{(v)}+\delta _{ij}(\sigma _i^2+\sigma _{}^2)$$
(3)
where $`R_{ij}^{(v)}=𝐯(𝐫_i)\widehat{𝐫}_i𝐯(𝐫_j)\widehat{𝐫}_j`$ and the second diagonal term is due to the noise. In linear theory, the “signal” part of the covariance matrix $`R^{(v)}`$ can be written as an integral over the density power spectrum

$$R_{ij}^{(v)}=\frac{H^2\mathrm{\Omega }_o^{1.2}}{2\pi ^2}P(k)W_{ij}^2(k)𝑑k,$$
(4)
where $`W_{ij}^2(k)`$ is a tensor window function calculated from the set of positions $`𝐫_i`$ of the galaxies with galaxy velocities weighted by their relative error (for more details see Feldman & Watkins (1994); Watkins & Feldman, (1995); Feldman & Watkins (1998)).

Typically we are given the a catalog of measured velocities $`(v_1,\mathrm{},v_N)`$ and wish to determine $`P(k)`$. Thus we can view $`L(v_1,\mathrm{},v_N;P(k))`$ as a likelihood functional; given $`(v_1,\mathrm{},v_N)`$ we use Eq. (2) to determine the likelihood that they were generated in a Universe with a particular $`P(k)`$.

An analysis of the type described above is susceptible to biases due to small–scale, nonlinear contributions to the galaxy velocities. The method of analysis developed in Paper I eliminates this bias by replacing the full set of $`N`$ line–of–sight velocities $`v_i`$ with moments $`u_i`$, which are designed so that the $`N`$th moment is the linear combination that carries the most information about small scales, the $`(N1)`$th moment is an independent linear combination which carries the second most small–scale information, etc. By using a subset of the first $`N^{}<N`$ moments in our analysis and discarding the rest, we can essentially filter out the small scale information which may carry nonlinear contributions. Our method is based on Karhunen–Loève methods of data compression (Kenney & Keeping, 1954; Kendall & Stuart, 1969); see also Tegmark, Taylor & Heavens (1997), designed to concentrate most of the information in a large set of data into a smaller, more manageable number of moments. However, our method puts a twist on this idea by concentrating unwanted information regarding small scales into a small set of moments which can then be discarded. What follows is a brief review of the mechanics of our method; for details, see Paper I.

Our method is based on linear data compression, so that the moment $`u_n`$ can be written in terms of the line–of–sight velocities $`v_i`$ as

$$u_n=\underset{j=1}{\overset{N}{}}(b_n)_jv_j,$$
(5)
where the $`b_n`$ are a set of vectors of length $`N`$. With this definition, the covariance matrix of the new moments is given as

$$\stackrel{~}{R}_{nm}=u_nu_m=\underset{i,j}{}(b_n)_iv_iv_j(b_m)_j=\underset{i,j}{}(b_n)_iR_{ij}(b_m)_j.$$
(6)
It is convenient to choose the $`b_n`$ so that the $`u_n`$ are linearly independent and of unit variance, so that $`\stackrel{~}{R}_{nm}`$ is the identity matrix. Note that this normalization will hold only for a particular matrix $`R_{ij}`$ and hence a particular power spectrum.

In order to find the vector $`b_n`$ such that the moment $`u_n`$ carries the maximum information about nonlinear scales, we assume a simple model for the power spectrum in which the amount of power on scales below that where density fluctuations have gone nonlinear is specified by a single parameter $`\theta _q`$. Given a single moment $`u_n`$, we can determine the value of $`\theta _q`$ to within a minimum variance given by $`\mathrm{\Delta }\theta _q^2=1/\stackrel{~}{F}_{qq}`$, where $`\stackrel{~}{F}_{qq}`$ is the $`qq`$th element of the Fisher information matrix, which in this case, and with the normalization assumed above, can be shown to take the form

$$\stackrel{~}{F}_{qq}=\underset{ij}{}\frac{1}{2}\left((b_n)_i\frac{R_{ij}}{\theta _q}(b_n)_j\right)^2$$
(7)
We can thus find the single moment $`u_n`$ that carries the maximum information about $`\theta _q`$ by finding the vector $`b_n`$ which minimizes $`\stackrel{~}{F}_{qq}`$ subject to the normalization condition discussed above, which functions as a constraint. After introducing a Lagrange multiplier, the minimization results in an eigenvalue problem

$$\underset{i,j,m}{}\left(L_{ki}^1\frac{R_{ij}}{\theta _q}L_{lj}^1\right)\left(L_{ml}b_m\right)=\underset{j}{}\lambda \left(L_{jk}b_j\right)$$
(8)
where $`L_{ij}`$ is the Cholesky decomposition of the covariance matrix, $`R_{ij}=_{p=1}^NL_{ip}L_{jp}`$.

Solving this eigenvalue problem gives us a set of $`N`$ orthogonal eigenvectors $`_jL_{ji}(b_n)_j`$ with eigenvalues $`\lambda _n`$. Each eigenvector has a corresponding moment $`u_n=_i(b_n)_iv_i`$. The eigenvalue $`\lambda _n`$ of a moment $`u_n`$ is related to the error bar $`\mathrm{\Delta }\theta _q`$ that one could place on $`\theta _q`$ using the single moment $`u_n`$, as can be seen by manipulating the equations above:

$$\frac{1}{\mathrm{\Delta }\theta _q^2}=\stackrel{~}{F}_{qq}=\underset{i,j}{}\frac{1}{2}(b_i\frac{R_{ij}}{\theta _q}b_j)^2=\underset{i,j}{}\frac{1}{2}(b_i\left(\lambda R_{ij}b_j\right))^2=\frac{\lambda ^2}{2}$$
(9)
so that $`\mathrm{\Delta }\theta _q=2/\lambda ^2`$. If we order the moments $`u_n`$ in order of increasing eigenvalue,

$$|\lambda _1||\lambda _2|\mathrm{}|\lambda _N|$$
(10)
then we can interpret each moment as carrying successively more information about $`\theta _q`$, with $`u_N`$ carrying the maximum possible amount of information. Since our goal is to produce a data set that is less sensitive to the value of $`\theta _q`$ than the original data, we should keep moments only up to some $`N^{}`$. The orthogonality of the eigenvectors ensures that the moments are statistically independent, as we assumed above. Thus if we compress the data by discarding the moments with a large value of $`|\lambda |`$, the information contained in those moments will be completely removed from the data. However, we would also like to keep as many moments as possible in order to retain the maximum information about large scales.

In order to choose a value of $`N^{}`$, we need to examine what error bar $`\mathrm{\Delta }\theta _q`$ we can put on the parameter $`\theta _q`$ using the compressed data. Since the moments are independent, we can write the Fisher matrix for the $`N^{}`$ moments that were not discarded as

$$\stackrel{~}{F}_{qq}=\underset{n=1}{\overset{N^{}}{}}\frac{1}{2}\lambda _n^2$$
(11)
so that the error bar that can be put on $`\theta _q`$ using the compressed data is given by

$$\mathrm{\Delta }\theta _q=\frac{1}{\sqrt{\stackrel{~}{F}_{qq}}}=\left[\frac{1}{2}\underset{n=1}{\overset{N^{}}{}}\lambda _n^2\right]^{1/2}$$
(12)
This suggests that $`N^{}`$ should be chosen by adding up the sum of the squares of the smallest eigenvalues until the desired sensitivity is reached. The criterion that we use is as follows: First, we estimate the actual size of the parameter $`\theta _q=\theta _{qo}`$ from peculiar velocity data. Then, we keep the largest number $`N^{}`$ moments that is still consistent with the requirement that $`\mathrm{\Delta }\theta _q\theta _{qo}`$. With this requirement, as long as our estimate of the true value of $`\theta _q`$ is correct, our final set of moments $`u_1\mathrm{}u_N^{}`$ will not contain enough information to distinguish the value of $`\theta _q`$ from zero.

## 3 Other Selection Criteria

Our method of selecting moments by their lack of information about small scales has the disadvantage of not discarding moments which have little information about any scale; this issue was briefly touched on in Paper I. Thus we have developed a second criterion for moment selection whereby we discard moments that are dominated by noise. That is, they have no cosmologically useful information.

Recall from Eq. 3 that the covariance matrix for the line–of–sight velocities is the sum of a “signal” part and a noise part. Since the covariance matrix for the moments $`u_n`$ is essentially a “rotation” of the velocity covariance matrix, this matrix can separated in a similar fashion,

$`\stackrel{~}{R}_{nm}`$ $`=`$ $`{\displaystyle \underset{ij}{}}(b_n)_iR_{ij}(b_m)_j={\displaystyle \underset{ij}{}}(b_n)_iR_{ij}^{(v)}(b_m)_j+{\displaystyle \underset{ij}{}}(b_n)_i\delta _{ij}(\sigma _i^2+\sigma _{}^2)(b_m)_j`$ (13)
$`=`$ $`\stackrel{~}{R}_{nm}^{(v)}+{\displaystyle \underset{i}{}}(b_n)_i(b_m)_i(\sigma _i^2+\sigma _{}^2)`$
where the second term is the noise contribution to the variance of the moment. Given that the $`b_n`$ are normalized such that the moments are independent and have unit variance, i.e. that $`\stackrel{~}{R}`$ is the identity matrix, we see that the quantity

$$\xi _n=\underset{i}{}(b_n)_i^2(\sigma _i^2+\sigma _{}^2)$$
(14)
is a measure of the fraction of the variance of the moment $`u_n`$ that is due to noise. If $`\xi _n1`$, then the moment has very little noise and should be retained. If, on the other hand, $`\xi _n1`$, then the value of the moment is mostly determined by the errors in the data and should be discarded.

Generally there is a correlation between a moment’s $`\xi _n`$ and its eigenvalue $`|\lambda _n|`$; moments that are most sensitive to small scales tend to be very noisy. Low noise moments tend to be those that probe large–scale power. This is due to the fact that the measurement errors in the velocities, which vary independently from galaxy to galaxy, are much more effective at masking small scale modes than large scale modes. However, we have found that some moments with small eigenvalues also have large noise; these moments carry little information about any scale.

The correlation between $`\xi _n`$ and $`|\lambda _n|`$ suggests that our two selection criterion can be somewhat redundant; eliminating noisy moments often can also accomplish the goal of removing moments with large $`|\lambda _n|`$. Similarly, eliminating moments with large $`|\lambda _n|`$ leaves one with moments which generally have smaller noise. Thus we will see below that applying the second criteria after we have already applied the first typically does not change the results of the analysis significantly.

Once one has a set of moments which have both small $`\xi _n`$ and $`|\lambda _n|`$, it is desirable to have a way of determining which scales each moment is most sensitive to. From Eq. 4 we recall that the “signal” part of the covariance matrix for the velocities $`R^{(v)}`$ is given by an integral of a tensor window function $`W_{ij}(k)`$ with the power spectrum. The “signal” part of the covariance matrix $`\stackrel{~}{R}`$ of the moments is given by a “rotation” of this tensor window function. Since $`\stackrel{~}{R}`$ is diagonal, this results in a scalar window function for each moment $`u_n`$,

$$W_n^2(k)=(b_n)_iW_{ij}^2(k)(b_n)_j$$
(15)
By examining the window functions for each moment, we can determine which scales the moments are sensitive to and confirm that our method is working. In principle, examination of the window functions could also provide a further criterion for the discarding of moments.

## 4 Synthetic Catalogs

The simulations used here are numerical models for the gravitational dynamics of collisionless particles in an expanding background. We are studying evolution of initial Gaussian perturbations in a matter–dominated universe. All the simulations are done with a particle–mesh (PM) code with $`256^3`$ particles in an equal number of grid points (Melott, 1986; Melott, Weinberg & Gott, 1988). More details about the peculiarities of the simulations used here can be seen in Melott & Shandarin (1993). Although the parameter $`\mathrm{\Omega }`$ appears in both the dynamics of the expanding background and as a parameter in the fit to the power spectrum shape in the CDM family of models (Bardeen et al. , 1986), they serve two different functions.

We ran simulations with a dynamical background $`\mathrm{\Omega }_0`$ = 1. or 0.34. These were normalized to an amplitude $`\sigma _8`$=0.93 at redshift moment $`z`$=0, but we also took data at $`z`$=1, which for our purposes can be described as studying a Universe with a lower perturbation amplitude normalization, and possibly a higher $`\mathrm{\Omega }_0`$. All models were interpreted with an assumed Hubble Constant H$`=100h`$km s-1Mpc-1 where $`h`$ = 2/3. When power spectra are parameterized in $`Mpc`$, the shape is dependent upon $`\mathrm{\Omega }_0h^2`$, which we set equal to 0.15, 0.35, and 1.0 for our low–$`\mathrm{\Omega }`$ tests. We have used some values which are inconsistent with other constraints in CDM linear theory in order to test our method over a wide range of values. We also ran another simulation with $`\mathrm{\Omega }_0`$ = 1, and $`\mathrm{\Omega }_0h^2`$ = 0.15; such models have been called $`\tau `$CDM in the past. The set with $`\mathrm{\Omega }_0`$ = 0.34, $`\mathrm{\Omega }_0h^2`$ = 0.15 is most consistent with a variety of findings at this time, but we do not wish to test our method only against currently favored cosmologies. There are a variety of alternative models in addition to the cosmological constant $`\lambda `$; all have very small and totally linear effects on large–scale velocities. We omit this in favor of wider exploration of parameter shifts which have large effects.

For testing our method we created synthetic redshift–distance catalogs from the 2563 $`N`$–body PM simulations. In these simulations, the box size was taken to equal 512 Mpc, or 34,133 km s-1 in redshift space for $`h`$=2/3.

Each of the points defined by the mesh represented a galaxy with a corresponding location and velocity. Testing the optimal moments method accurately requires proper modeling of statistical errors so we wanted to be sure to include the effects of cosmic variance and scatter in distance indicators. We chose three coordinates of scale 1/6, 1/2 and 5/6 the box width. An exhaustive permutation of combining these coordinates results in 27 center locations within the box, each corresponding to an individual synthetic catalog as described next. About each central location, an annular volume was defined by the redshift range 500 km s$`{}_{}{}^{1}<`$ z $`<`$ 10,000 km s-1. If the volume intercepted the boundaries of the box then we appealed to the periodic boundary conditions of the simulation and included galaxies from the opposite box side. Within each of the defined regions, we selected $``$ 1000 galaxies under the assumption of a radial selection function and the additional requirement that there is a zone of avoidance below the galactic latitude $`|b|<10^{}`$. The selection function was chosen generically (we ignored galactic properties) to mimic existing popular redshift–distance surveys (i.e. the SFI survey; da Costa et al. (1996)). The effect of scatter in distance indicators was replicated by adding a random error to each peculiar velocity drawn from a Gaussian distribution of width 10% of the galaxy redshift distance.

In the end, for each simulation box we have 27 surveys sampling the simulation that give information about the positions and radial velocities of the galaxy distribution in some volume. We analyzed these surveys by using the actual positions (i.e. no errors) and by perturbing the velocities with a 10% Gaussian error. The value of sampling the simulation in this way is that we are able to model the effects of cosmic variance. The errorbars of the unperturbed catalogs are predominantly from the cosmic variance, whereas those of the perturbed, 10% catalogs include the effects of both cosmic variance and the inaccuracies of the distance indicators.

## 5 Results

The main purpose of the formalism we presented here and in Paper I was to allow the removal or filtering of small scale noise while keeping the large scale signal. To test the success of the formalism we have created synthetic surveys from simulations with known parameters, specifically, $`\mathrm{\Gamma }`$, the CDM power spectrum shape parameter, and $`\beta `$, its amplitude (the values of the parameters we simulated is given in the table 1 below.) To compare our method with the full analysis method, we reemphasize that the optimal moment analysis presented here allows for two semi–independent methods of cleaning up a survey: 1) Ordering the moments by their eigenvalues (Eq. 10); and 2) Removing the noisiest moments (Eq. 14). In Figs. 1 and 36 we show the comparison between choosing the modes least susceptible to small scale signal (open triangles); those that are least susceptible to small scale signal and are not noisy (open squares); and the full analysis (that is, keeping all moments, the usual analysis, solid triangles). We see that the full analysis fails to recover the “true” parameters by a significant amount ($`4\sigma `$ for no errors and $`>2\sigma `$ for 10% errors). In contrast, the mode analysis recovers the values of the parameters very well, with or without the removal of the noisy moments.

For each one of the models we simulated we have extracted 27 catalogs from each simulation box, as described in Sec. 4. The points and errorbars in the figures are the maximum likelihood mean and standard deviation of the mean for the analysis of all catalogs. Each one of the catalogs were analyzed using the full maximum likelihood analysis, keeping all moments; the maximum likelihood analysis discarding large eigenvalue modes; and the maximum likelihood mode analysis without the noisiest moments ($`\xi >0.95`$).

In Fig. 2 we show the value of the estimated parameters as a function of the $`\mathrm{\Sigma }\lambda ^2`$ (see Eq. 12) where we see that as the number of modes is increased, we get closer and closer to the “true” value. When we keep more than the number of moments that corresponds to the fulfillment of our criterion (Eq. 12), the values start diverging from the “true” results. This is due to the fact that small–scale modes that have become nonlinear are introducing a bias. This tendency of the full analysis to systematically overestimate the parameter values can be seen in the analyses done for simulations with various cosmological parameters, Figs. 1 and 36.

We have experimented with the choice of $`N^{}`$, the number of modes to keep, as discussed in the text after Eq. 12. This choice depends on our power spectrum and more specifically on $`k_{nl}`$, the wavenumber of the largest scale for which density perturbations have become nonlinear (see Paper I). We chose $`k_{nl}`$ by comparing the linear power spectrum from the initial conditions of the simulation to the power spectrum at the end of the run. $`k_{nl}`$ is where the power spectra started to diverge. In general we found that $`k_{nl}0.2`$, though choosing $`0.15<k_{nl}<0.4`$ did not affect our results significantly. Further, as can be seen in Fig. 2, the $`N^{}`$ choice need not be finely tuned.

As was discussed in the text, the reason for the full analysis failure to recover the “true” parameters when the mode analysis succeeds so well can be shown by looking at the window functions themselves. In fig. 7 we show the window functions corresponding to the five lowest eigenvalues and lowest noise (lower left panel). Clearly, these probe only large scales. As we move up the panels we see the window functions with larger noise components not removed, whereas when we move to the right we see window functions corresponding to larger eigenvalues. Here the reasons for the particular choices for our criteria Eqs. (10) and (12) become clear. As the eigenvalues or the noise level become large, the window functions generally probe more small scale and less of large scale modes. Since we are primarily interested in large scale information, discarding the noisy, high $`\lambda `$ modes allows us to discard small–scale signal that might interfere with with our analysis.

One more advantage the formalism provides is efficiency. For a catalog of $`1,000`$ galaxies it takes the mode analysis about one hour CPU time, whereas it takes the full analysis about seven hours to complete. The differences are more dramatic for larger surveys: A 5,000 galaxy catalog completes in about 30 hours with the mode analysis and about 1,300 hours of CPU time for the full analysis. All runs were done on the Origin2000 at the NCSA, University of Illinois, Urbana–Champaign.

As was mentioned in the analysis (Sec. 3) there is a general correlation between the $`n^{\mathrm{th}}`$ moment’s eigenvalue $`\lambda _n`$ and the noise $`\xi _n`$ associated with it (dots in Fig. 8). As can be seen in the figure, the correlation is not perfect, that is, there are low eigenvalues with large noise component, but in general as we move to larger eigenvalues (large $`n`$) the noise component is larger. The line in the figure represents the running mean of the noise which shows clearly the correlations between the noise and the eigenvalue.

In Fig. 9 we show the contours that contain 68% and 94% of the total likelihood for six typical catalogs. The diamond shows the maximum likelihood results, whereas the asterix in each panel shows the “true” values of the parameters. These contours allow us to estimate the uncertainty in the maximum likelihood values obtained from the analysis of a single catalog, as is the case when analyzing observational data. From the figures it is clear that the uncertainties obtained in this way are comparable to those we get from the Monte–Carlo simulations. In general, when we try to test the reliability of results from an observational data set, we apply our formalism to mock catalogs extracted from N–body simulations as was done here. This compatibility between the uncertainties obtained in two different ways gives us confidence that using the likelihood contours will give us an accurate assessment of the uncertainties of our maximum likelihood values when we apply our method to real catalogs.

## 6 Conclusions

We have described the power and elegance of a new statistic that was designed and formulated in order to address a crisis in the analysis of proper distance cosmological surveys. We have shown that our formalism mostly overcomes the problems with the traditional analysis of the data. Whereas the full maximum likelihood analysis tends to overestimate the values of the parameters that describe the power distribution on large scale, our mode analysis makes very accurate estimates of these parameters.

The formalism presented here assumes Gaussian statistics. The natural question should be: Can the deviations from Gaussianity caused by the collapse of perturbations interfere with the removal of small scale power and introduce additional unpredictable biases? As the results in Sec. 5 indicate, deviations from Gaussianity do not have a measurable effect and the effectiveness of filtering small–scale power is unbiased. Further, we have explored in detail such issues as moment and noise selection, the window functions’ effectiveness and criteria for which modes to keep.

As was shown in Paper I and in more detail here, the formalism we presented is highly adaptive and versatile. It can be applied surveys with any geometry and density, and since it retains maximum information should be particularly useful for sparse data such as that obtained in cluster peculiar velocity surveys. Overall, we consider this method to be a significant improvement over previous methods used for the analysis of peculiar velocity data.

We wish to thank Jim Fry, Roman Juszkiewicz and Avishai Dekel for illuminating conversations. HAF and ALM wish to acknowledge support from the National Science Foundation under grant number AST–0070702, the University of Kansas General Research Fund and the National Center for Supercomputing Applications for allocation of computer time. This research has been partially supported by the Lady Davis Foundation at the Hebrew University, Jerusalem, Israel and by the Institute of Theoretical Physics at the Technion, Haifa, Israel.


format: text

text:
# Searching for multi-boson effects and space-time asymmetries in particle production

## 1 Introduction

The correlations of particles at small relative velocities are widely used to study space-time characteristics of the production processes. Particularly, for non-interacting identical particles, like photons, this technique is called intensity or particle interferometry. In this case the correlations appear solely due to the effect of quantum statistics (QS) . Similar effect was first used in astronomy to measure the angular radii of stars by studying the dependence of the two-photon coincidence rate on the distance between the detectors (HBT effect ). In particle physics the QS interference was first observed as an enhanced production of the pairs of identical pions with small opening angles (GGLP effect ). Later on, similar to astronomy, Kopylov and Podgoretsky  suggested to study the interference effect in terms of the correlation function.111Note that though both the KP and HBT methods are based on the QS interference, they represent just orthogonal measurements . The former, being the momentum-energy measurement, yields the space-time picture of the source, while the latter does the opposite. In particular, the HBT method provides the information about the angular size of a star but, of course, - no information about the star radius or its lifetime.

The effect of QS is usually considered in the limit of a low phase-space density such that the possible multi-particle effects can be neglected. This approximation seems to be justified by present experimental data which does not point to any spectacular multi-boson effects neither in single-boson spectra nor in two-boson correlations. These effects can however clearly manifest themselves in some rare events (e.g., those with large pion multiplicities) or in the eventually overpopulated regions of momentum space; see, e.g., , references therein and Section 3.

The particle correlations are also influenced by the effect of particle interaction in the final state (FSI) . Thus the effect of the Coulomb interaction dominates the correlations of charged particles at very small relative momenta (of the order of the inverse Bohr radius of the two-particle system), respectively suppressing or enhancing the production of particles with like or unlike charges. Regarding the effect of the strong FSI, it is quite small for pions, while for nucleons it is often a dominant one due to the very large magnitude of the s-wave singlet scattering length of about 20 fm.

Though the FSI effect complicates the correlation analysis, it is an important source of information allowing one to measure the space-time characteristics of the production process even with the help of non-identical particles . Moreover, the unlike particle correlations, contrary to those of identical particles, are sensitive to the relative space-time asymmetries in their production, e.g. \- to the relative time delays, thus giving an important information not accessible in the standard interferometry measurements . In Section 4 we briefly formulate the theory of these correlations and demonstrate the possibilities of the corresponding correlation technique.

## 2 Formalism

As usual, we will assume sufficiently small phase-space density of the produced multi-particle system, such that the correlation of two particles emitted with a small relative velocity in nearby space-time points is influenced by the effects of their mutual QS and FSI only.222 This assumption may be not valid in the case of low energy heavy ion reactions when the particles are produced in a strong Coulomb field of residual nuclei. To deal with this field a quantum $`adiabatic`$ approach can be used .  We define the ideal two-particle correlation function $`R(p_1,p_2)`$ as a ratio of the differential two-particle production cross section to the reference one which would be observed in the absence of the effects of QS and FSI. In heavy ion or high energy hadronic collisions we can neglect kinematic constraints and most of the dynamical correlations and construct the reference distribution by mixing the particles from different events.

Assuming the momentum dependence of the one-particle emission probabilities in- essential when varying the particle 4-momenta $`p_1`$ and $`p_2`$ by the amount characteristic for the correlation due to QS and FSI (smoothness assumption), i.e. assuming that the components of the mean space-time distance between particle emitters are much larger than those of the space-time extent of the emitters, we get the well-known result of Kopylov and Podgoretsky for identical particles, modified by the substitution of the plane wave $`e^{ip_1x_1+ip_2x_2}`$ by the nonsymmetrized Bethe-Salpeter amplitudes in the continuous spectrum of the two-particle states $`\psi _{p_1p_2}^{S(+)}(x_1,x_2)`$, where $`x_i=\{t_i,𝐫_i\}`$ are the 4-coordinates of the emission points of the two particles and $`S`$ is their total spin . At equal emission times in the two-particle c.m.s. ($`t^{}=t_1^{}t_2^{}=0`$) this amplitude coincides (up to an unimportant phase factor due to the c.m.s. motion) with a stationary solution of the scattering problem $`\psi _𝐤^{}^{S(+)}(𝐫^{})`$, where $`𝐤^{}=𝐩_1^{}=𝐩_2^{}`$ and $`𝐫^{}=𝐫_1^{}𝐫_2^{}`$ (the minus sign of the vector $`𝐤^{}`$ corresponds to the reverse in time direction of the emission process). The Bethe-Salpeter amplitude can be usually substituted by this solution (equal time approximation).333The equal time approximation is valid on condition  $`|t^{}|m_{2,1}r^2`$ for $`\mathrm{sign}(t^{})=\pm 1`$ respectively. This condition is usually satisfied for heavy particles like kaons or nucleons. But even for pions, the $`t^{}=0`$ approximation merely leads to a slight overestimation (typically $`<5\%`$) of the strong FSI effect and, it doesn’t influence the leading zero–distance ($`r^{}|a|`$) effect of the Coulomb FSI. Then, for nonidentical particles,

$$R(p_1,p_2)=\underset{S}{}\rho _S|\psi _𝐤^{}^{S(+)}(𝐫^{})|^2_S.$$
(1)
Here the averaging is done over the emission points of the two particles in a state with total spin $`S`$ populated with the probability $`\rho _S`$, $`_S\rho _S=1`$.444 For unpolarized particles with spins $`s_1`$ and $`s_2`$ the probability $`\rho _S=(2S+1)/[(2s_1+1)(2s_2+1)]`$. Generally, the correlation function is sensitive to particle polarization. For example, if two spin-1/2 particles are emitted with polarizations $`𝐏_1`$ and $`𝐏_2`$ then  $`\rho _0=(1𝐏_1𝐏_2)/4`$ and $`\rho _1=(3+𝐏_1𝐏_2)/4`$.  For identical particles, the amplitude in Eq. (1) should be properly symmetrized:

$$\psi _𝐤^{}^{S(+)}(𝐫^{})[\psi _𝐤^{}^{S(+)}(𝐫^{})+(1)^S\psi _𝐤^{}^{S(+)}(𝐫^{})]/\sqrt{2}.$$
(2)
Particularly, for non-interacting identical particles the characteristic feature of the correlation function is the presence of the interference maximum or minimum at small $`|𝐪|`$ changing to a horizontal plateau at sufficiently large $`|𝐪|`$. For example, assuming that particles are produced independently and that the space and time limitation of the production process is effectively described by the one-particle probability of a Gaussian form with the corresponding dispersions $`r_0^2`$ and $`\tau _0^2`$, Eqs. (1) and (2) yield:

$$R(p_1,p_2)=1+\underset{S}{}(1)^S\rho _S\mathrm{exp}[r_0^2q_T{}_{}{}^{2}(r_0{}_{}{}^{2}+v^2\tau _0{}_{}{}^{2})q_L{}_{}{}^{2}].$$
(3)
Here $`q_T`$ and $`q_L`$ are the transverse and longitudinal components of the relative 3-momentum $`𝐪`$ with respect to the pair velocity vector $`𝐯`$ ($`q_T=2k_T^{}`$ and $`q_L=q_0/v=2\gamma k_L^{}`$, where $`\gamma `$ is the pair Lorentz factor). We see that, due to the relation $`q_0=\mathrm{𝐯𝐪}vq_L`$, the correlation function at $`v\tau _0>r_0`$ substantially depends on the direction of the vector $`𝐪`$ even in the case of spherically symmetric spatial form of the production region, the interferometric radii squared in the $`T`$ and $`L`$directions being $`r_T^2=r_0^2`$ and $`r_L^2=r_0^2+v^2\tau _0^2`$. Generally, the directional dependence of the correlation function can be used to determine both the characteristic emission time and the form of the production region .

It should be noted that particle correlations at high energies usually measure only a small part of the space-time emission volume since, due to substantially limited decay momenta of few hundred MeV/c, the sources, despite their fast longitudinal motion, emit the correlated particles with nearby velocities mainly at nearby points in the c.m.s. of particle pair. The dynamical examples are sources-resonances, colour strings or hydrodynamical expansion. To get rid of a fast longitudinal motion, sometimes the longitudinally comoving system (LCMS) is introduced, in which each pair is emitted transverse to the reaction axis.

## 3 Multi-boson effects

The two-body formalism of Section 2 assumes that the phase-space density $`f`$ of the produced multi-particle system is small compared to unity so that the multi-particle correlations can be neglected. The mean phase-space density of the pions of a given type at a given momentum p (rapidity $`y`$ and transverse momentum $`𝐩_t`$) can be estimated, in the low density limit, as a number of the pions interfering with a pion of momentum p and building the Bose-Einstein (BE) enhancement . Using the usual Gaussian parametrization for the correlation function in the LCMS: $`R(p_1,p_2)=1+\lambda \mathrm{exp}(r_x^2q_x^2r_y^2q_y^2r_z^2q_z^2)`$, and parametrizing the cut-off of transverse momenta by an exponential law in transverse mass $`m_t`$, we can estimate the mean pion phase-space density as

$$f_𝐩\lambda \pi ^{3/2}\frac{coshy}{V}\frac{d^3n}{d^3𝐩}=\lambda \frac{\sqrt{\pi }}{2}\frac{\mathrm{exp}[(m_tm)/T]}{VT(T+m)m_{}}\frac{dn}{dy},$$
(4)
where $`V=r_xr_yr_z`$ is the interference volume. The mean density is maximal for soft pions ($`p_ty0`$). At present energies it is typically less than $`0.2`$ thus indicating rather small multi-boson effects. These effects can be expected small also in future heavy ion exeperiments at RHIC and CERN since presently the LCMS interference volume $`V`$ seems to scale with the density $`dn/dy`$ pointing to the freeze-out of the particles at a constant phase-space density. Nevertheless, the multi-boson effects can show up in certain classes of events or in some regions of momentum space.

The multi-boson effects can be practically treated provided that we can neglect particle interaction in the final state and assume independent particle emission, supplemented by the requirement of a universal single-particle emission function, independent of the origin of single-particle sources. The multi-particle emission function then being a product of the single-particle ones. The latter allows to calculate the Bose-Einstein effect on the original multiplicity distribution or boson spectra with the help of the BE weights $`\omega _n`$  which can be expressed through the so called cumulants or differential cumulants . The calculation of all the cumulants is generally a difficult task. Thus, for realistic transport models used to predict particle production in ultra-relativistic heavy-ion collisions, the numerical limitations allow to determine only a few lowest order cumulants (up to about the fourth order) . Fortunately, they are sufficient for typical rather moderate pion freeze-out phase-space densities .

The analytical calculation of the cumulants of all orders is possible in a simple model in which the pions are instantaneously emitted according to the Gaussian emission function (in Wigner-like phase-space) characterized by the dispersions $`\mathrm{\Delta }^2`$ and $`r_0^2`$ in the momentum and ordinary space respectively .555 Note that in this model the two-pion correlation function in the low density limit is given by Eq. (3) with $`S=0`$, $`\tau _0=0`$ and $`r_0^2`$ substituted by the dispersion $`r_0^21/(4\mathrm{\Delta }^2)`$ of the centers of the elementary emitters, where $`1/(4\mathrm{\Delta }^2)`$ represents the minimal possible dispersion of the Wigner-like spatial coordinates due to emitter finite sizes.

Assuming a Poissonian distribution of the original boson multiplicities with the mean multiplicity $`\eta `$, the resulting one, at large $`n`$, tends to the BE distribution with the mean multiplicity $`\xi /(1\xi )`$, where $`\xi =\eta /(r_0\mathrm{\Delta }+1/2)^3`$ can be considered as an intrinsic phase space parameter. Its maximal value of 1 corresponds to the explosion of the multiplicity distribution. Comparing the model phase-space density $`\stackrel{~}{f}_{𝐩=0}\eta /(\sqrt{2}r_0\mathrm{\Delta })^3`$ with the experimental estimate of about 0.2, we get for the density parameter $`\xi 0.40.5`$. It was shown  that at such densities an approximate $`\xi `$-scaling of the ratio $`n/\eta `$ takes place; $`n`$ becomes close to the asymptotic scaling value of $`\xi /(1\xi )`$ only at $`\xi >0.9`$.

Regarding the influence of the BE effect on the single-boson spectrum, at sufficiently large momenta, when the local density $`f_𝐩`$ remains small even at large multiplicities, it is dominated by the contribution of the original spectrum, otherwise, at large local densities, it is determined by the asymptotic large-density spectrum with the original dispersion $`\mathrm{\Delta }^2`$ substituted by $`\mathrm{\Delta }/(2r_0)`$ . Experimentally the effect of BE ”condensate” was searched for at SPS CERN as a low-$`p_t`$ enhancement, however, with rather uncertain results.

As for the two-boson correlation function, it is well known that for a given multiplicity $`n`$, the intercept $`R_n(0)`$ decreases and the correlation function width increases with the increasing $`n`$ or decreasing momentum $`p`$, both corresponding to the increasing local density. In Fig. 1 we show  the intercept as a function of the multiplicity $`n`$ and the density parameter $`\xi _{n,𝐩}=\xi _n\mathrm{exp}(𝐩^2/2\mathrm{\Delta }^2)`$, $`\xi _n=n/(r_0\mathrm{\Delta }+1/2)^3`$, for several values of the mean momentum: $`p=0`$, 0.1, 0.2 and 0.4 GeV/c. As expected, the intercept is practically constant at low local densities (small $`n`$ or high $`p`$). As condensate develops, the intercept sharply falls down. The sharpness of this drop is however less pronounced at higher momenta even if plotted as a function of the local density. Clearly, this lack of density scaling is related to a strong decrease of the condensate contribution with the increasing momentum. To demonstrate the possibility of the observation of the condensate effect for soft pions produced in certain high multiplicity events not following the ordinary proportionality rule between the freeze-out phase space density and pion multiplicity, in Fig. 1 we indicate by the arrows the intercept values corresponding to $`\xi _n=3\xi =1.5`$ ($`n3n`$).

In the inclusive case corresponding to the original Poissonian multiplicity distribution, the correlation function intercept is equal to 2 for any pion momenta or the local densities  and its width logarithmically increases with the increasing condensate contribution .

The results of the considered model should not be taken, however, too literally since:

a) in contradiction with the experimental indications on a constant freeze-out phase space density, in the model there is no correlation between the emission volume and pion multiplicity;

b) the static character of the model can be justified (neglecting the transverse expansion) in a limited rapidity region only. Thus the pions with a rapidity difference greater than about unity have to be considered as originating from different static sources;

c) since the inclusive spectra at sufficiently high mean multiplicity are dominated by the condensate contributions in both extreme cases of very narrow or very wide multiplicity distributions of the originally uncorrelated bosons, the intercept of the inclusive correlation function is likely to be less than the value of 2 corresponding to the Poissonian case;

f) for identical charged pions, the BE effects are suppressed due to the Coulomb repulsion. Since this repulsion is important only in a weakly populated region of very small relative momenta, the suppression of the global BE weights $`\omega _n`$ is rather small. For example, for $`\omega _2`$ this suppression, being about $`(ar_0\mathrm{\Delta }^2)^1`$, is usually less than one per mill. The Coulomb distortion of the global multi-boson effects is therefore negligible in the rare gas limit. Nevertheless, since the Coulomb repulsion destroys the formation of the condensates made up from positive and negative pions in the disjoint phase space regions, it can lead to noticeable differences between charged and neutral pions in dense systems. Particularly, we can expect a decrease of the charge-to-neutral multiplicity ratio with the increasing phase space density.

## 4 Measuring the relative space-time asymmetries

The correlation function of two nonidentical particles, compared with the identical ones, contains a principally new piece of information on the relative space-time asymmetries in particle emission . This is clearly seen in the case of neutral particles when the two-particle amplitude $`\psi _𝐤^{}^{S(+)}(𝐫^{})`$ takes on the form

$$\psi _𝐤^{}^{S(+)}(𝐫^{})=e^{i𝐤^{}𝐫^{}}+\varphi _k^{}^S(r^{}),$$
(5)
where the scattered wave $`\varphi _k^{}^S(r^{})`$, in the considered region of small relative momenta, is independent of the directions of the vectors $`𝐤^{}`$ and $`𝐫^{}`$. Inserting Eq. (5) into the formula (1) for the correlation function, we can see that the latter is sensitive to the relative space-time asymmetry due to the odd term $`\mathrm{sin}𝐤^{}𝐫^{}`$. Particularly, it allows for a measurement of the mean relative delays $`tt_1t_2`$ in particle emission. To see this, let us make the Lorentz transformation from the rest frame of the source to the c.m.s. of the two particles: $`r_L^{}=\gamma (r_Lvt)`$, $`r_T^{}=r_T`$. Considering, for simplicity, the behavior of the vector $`𝐫^{}`$ in the limit $`|vt|r`$, we see that this vector is only slightly affected by averaging over the spatial distance $`r|vt|`$ of the emission points in the rest frame of the source so that $`𝐫^{}\gamma 𝐯t`$. Therefore, the vector $`𝐫^{}`$ is nearly parallel or antiparallel to the velocity vector $`𝐯`$ of the pair, depending on the sign of the time difference $`t\mathrm{\Delta }t=t_1t_2`$. The sensitivity to this sign is transferred to the correlation function through the odd in $`𝐤^{}𝐫^{}\gamma 𝐤^{}𝐯t`$ term provided the sign of the scalar product $`𝐤^{}𝐯`$ is fixed.

For charged particles there arise additional odd terms due to the confluent hypergeometrical function $`F(\alpha ,1,z)=1+\alpha z/1!^2+\alpha (\alpha +1)z^2/2!^2+\mathrm{}`$, modifying the plane wave in Eq. (5):

$$\psi _𝐤^{}^{S(+)}(𝐫^{})=\mathrm{e}^{i\delta }\sqrt{A_c(\eta )}\left[\mathrm{e}^{i𝐤^{}𝐫^{}}F(i\eta ,1,i\rho )+\varphi _{ck^{}}^S(r^{})\right],$$
(6)
where $`\rho =𝐤^{}𝐫^{}+\mathrm{k}^{}r^{}`$, $`\eta =(k^{}a)^1`$, $`a=(\mu z_1z_2e^2)^1`$ is the Bohr radius of the two-particle system taking into account the sign of the interaction ($`z_ie`$ are the particle electric charges, $`\mu `$ is their reduced mass), $`\delta =\text{arg}\mathrm{\Gamma }(1+i\eta )`$ is the Coulomb s-wave shift and $`A_c(\eta )=2\pi \eta /[\mathrm{exp}(2\pi \eta )1]`$ is the Coulomb penetration factor.666This factor substantially deviates from unity only at $`k^{}<2\pi /|a|`$ (e.g., at $`k^{}<22`$ MeV/c for two protons). Note that for the distances $`r^{}>|a|`$ the confluent hypergeometrical function becomes important and compensates the deviation of the Coulomb factor from unity except for the classically forbidden region of $`k^{}<(|a|r^{}/2)^{1/2}`$, narrowing with the increasing $`r^{}`$. Thus the FSI practically vanishes if at least one of the two particles comes from a long lived source ($`\eta ,\eta ^{},\mathrm{\Lambda },K_s^0,\mathrm{}`$).  Clearly, at a given distance $`r^{}`$, the effect of the odd component in the Coulomb wave function is of increasing importance with a decreasing Bohr radius of the particle pair, i.e. for particles of greater masses or electric charges. At low energies, the sensitivity of the correlation to the odd component can be somewhat modified due to Coulomb interaction with the residual charge .

It is clear that in the case of a dominant time asymmetry, $`v|t||r_L|`$, a straightforward way to determine the mean time delay $`t`$ is to measure the correlation functions $`R_+(𝐤^{}𝐯0)`$ and $`R_{}(𝐤^{}𝐯<0)`$. Depending on the sign of $`t`$, their ratio $`R_+/R_{}`$ should show a peak or a dip in the region of small $`k^{}`$ and approach 1 at large values of $`k^{}`$. As the sign of the scalar product $`𝐤^{}𝐯`$ is practically equal to that of the difference of particle velocities $`v_1v_2`$ (this equality is always valid for particles of equal masses), the sensitivity of the correlation functions $`R_+`$ and $`R_{}`$ to the sign of the difference of particle emission times has a simple classical explanation. Clearly, the interaction between the particles in the case of an earlier emission of the faster particle will be weaker compared with the case of its later emission (the interaction time being longer in the latter case leading to a stronger correlation). This expectation is in accordance with Eqs. (1) and (6) at $`k^{}0`$, $`r^{}|a|`$ and $`|\varphi _{ck^{}}^S(r^{})|1`$, when (the arrow indicates the limit $`v|t||r_L|`$):

$$R_+/R_{}1+2\frac{r_L^{}}{a}12\frac{\gamma v(t_1t_2)}{a}.$$
(7)
The sensitivity of the $`R_+/R_{}`$ correlation method to the mean relative time shifts (introduced ad hoc) was studied  for various two-particle systems simulated in $`Pb+Pb`$ collisions at SPS energy using the event generator VENUS 5.14 . The scaling of the effect with the space-time asymmetry and with the inverse Bohr radius $`a`$, indicated by Eq. (7), was clearly illustrated for the $`K^+K^{}`$ system ($`a=110`$ fm) and for the like- and unlike-sign $`\pi K`$, $`\pi p`$ and $`Kp`$ systems ($`a=\pm 249`$, $`\pm 226`$ and $`\pm 84`$ fm respectively). It was concluded that for sufficiently relativistic pairs ($`\gamma v>0.5`$) the $`R_+/R_{}`$ ratio can be sensitive to the shifts in the particle emission times of the order of a few fm/c. Motivated by this result the $`R_+/R_{}`$ method was recently applied to the $`K^+K^{}`$ system simulated in a two-phase thermodynamical evolution model and the sensitivity was demonstrated to the production of the transient strange quark matter state even if it decays on strong interaction time scales .

The method sensitivity was also studied for AGS and SPS energies using the transport code RQMD v2.3 . Thus at SPS energy the central $`Pb+Pb`$ collisions have been simulated and the $`\pi ^+K^+`$, $`\pi ^+p`$ and $`K^+p`$ correlations have been studied . To get rid of the effect of a fast longitudinal motion, the study was done in the longitudinally co-moving system (LCMS) in which the pair is emitted transverse to the reaction axis so that $`v=v_{}`$, $`r_L=\mathrm{\Delta }x`$ and

$$\mathrm{\Delta }x^{}=\gamma _{}(\mathrm{\Delta }xv_{}\mathrm{\Delta }t),\mathrm{\Delta }y^{}=\mathrm{\Delta }y,\mathrm{\Delta }z^{}=\mathrm{\Delta }z.$$
(8)
The simulated correlation functions $`R_+`$, $`R_{}`$ and their ratios are plotted in Fig. 2. We can see that for $`\pi ^+p`$ and $`\pi ^+K^+`$ systems these ratios are less than unity at small values of $`qk^{}`$, while for $`K^+p`$ system the ratio $`R_+/R_{}`$ practically coincides with unity. These results well agree with the mean values of $`\mathrm{\Delta }t`$, $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }x^{}`$ presented in Table 1 ($`\mathrm{\Delta }y\mathrm{\Delta }z0`$). It can be seen from Eqs. (7) and (8) that the absence of the effect in the $`R_+/R_{}`$ ratio for the $`K^+p`$ system is due to practically the complete compensation of the space and time asymmetries leading to $`\mathrm{\Delta }x^{}0`$. For $`\pi ^+p`$ system the effect is determined mainly by the x-asymmetry. For $`\pi ^+K^+`$ system both the x- and time-asymmetries contribute in the same direction, the latter contribution being somewhat larger. The separation of the relative time delays from the spatial asymmetry is, in principle, possible (see Eq. (8)) by studying the ratio $`R_+/R_{}`$ in different intervals of the pair velocity.

At AGS energy the $`Au+Au`$ collisions have been simulated and the $`\pi ^+p`$ correlations have been studied in the projectile fragmentation region where proton directed flow is most pronounced and where the proton and pion sources are expected to be shifted

relative to each other both in the longitudinal and in the transverse (flow) directions in the reaction plane. It was demonstrated that a modification of the $`R_+/R_{}`$ method ($`\pm `$ corresponding now to the signs of the respective components $`k_i^{}`$) is sufficiently sensitive to reveal these shifts . Recently, the shifts predicted by RQMD were confirmed by applying the $`R_+/R_{}`$ method to the experimental AGS data .

At low energies, the particles in heavy ion collisions are emitted with the characteristic emission times of tens to hundreds fm/c so that the observable time shifts should be of the same order . In fact the $`R_+/R_{}`$ method has been successfully applied to study proton-deuteron correlations in several heavy ion experiments at GANIL . It was observed, in agreement with the coalescence model, that deuterons are on average emitted earlier than protons.

## 5 Conclusion

Using the analytically solvable Gaussian model, we have discussed the influence of the multi-boson effects on boson multiplicities, single-boson spectra and two-boson correlations, including an approximate scaling behavior of some of their characteristics with the phase space density. Though these effects are hardly to be observable in typical events of heavy ion collisions in present and perhaps also in future heavy ion experiments, they can show up in certain classes of events, e.g., in those with high pion multiplicities.

We have shown that unlike particle correlations, compared with those of identical particles, contain a principally new piece of information on the relative space-time asymmetries in particle emission, thus allowing, in particular, a measurement of the mean relative delays in particle emission at time scales as small as $`10^{23}`$ s. To determine these asymmetries, the unlike particle correlation functions $`R_+`$ and $`R_{}`$ have to be studied separately for positive and negative values of the projection of the relative momentum vector in pair c.m.s. on the pair velocity vector or, generally, - on any direction of interest. We have discussed here the results of recent studies of these correlation functions for a number of two-particle systems simulated with various event generators. It was shown that the $`R_+/R_{}`$ ratio is sufficiently sensitive to the relative space-time asymmetries arising due to the formation of the quark-gluon plasma and strangeness distillation and even to those expected in the usual dynamical scenarios at AGS and SPS energies. As to the detection of the unlike particles with close velocities ($`p_1/m_1p_2/m_2`$), there is no problem with the two-track resolution since these particles, having either different momenta or different charge-to mass ratios, have well separated trajectories in the detector magnetic field. For the same reason, however, a large momentum acceptance of the detector is required.

This work was supported by GA AV Czech Republic, Grant No. A1010601 and by GA Czech Republic, Grant No. 202/98/1283.


text:
# Possible spin triplet superconductivity in NaxCoO₂⋅𝑦H20

## Abstract

Combining symmetry based considerations with inputs from available experimental results, we make the case that a novel spin-triplet superconductivity triggered by antiferromagnetic fluctuations may be realized in the newly discovered layered cobaltide NaxCoO$`{}_{2}{}^{}y`$H2O. In the proposed picture, unaccessable via resonating-valence-bond physics extrapolated from half-filling, the pairing process is similar to that advanced for Sr2RuO4, but enjoys a further advantage coming from the hexagonal structure of the Fermi-surface which gives a stronger pairing tendency.

Anderson’s suggestion that a doped antiferromagnet on a 2d square lattice be best described in terms of resonating valence bond (RVB) states, set the direction for a huge subsequent effort pursuing this picture in an attempt to understand superconductivity in the cuprate oxide compounds. While this has lead to the development of new avenues in the physics of low dimensional condensed matters, conclusive evidence of such novel states has not turned up. Meanwhile, renewed interests in RVB physics in systems on triangular lattices (TL) have been emerging recently - though the S=1/2 antiferromagnetic (AF) Heisenberg model on a TL is generally regarded as having a ground state with long range order, closely related systems such as the undoped quantum dimer model and the $`tJ`$ model on a TL with the hopping sign $`t>0`$ have been shown to exhibit a certain amount of stability of RVB-like states.

The discovery by Takada et al of superconductivity in NaxCoO2 $`y`$ H2O ($`x0.35,y1.3`$) , a system made of stacked layers of cobalt ions sitting on a 2d TL and surrounded by octahedrally displaced oxygen atoms will no doubt spark new interest in this line of investigation. Superconductivity was achieved at Tc=5K by oxidation of the sodium ions (changing the effective charge density) followed by intercalation of the CoO2 layers in the mother compound NaCo2O4 with H2O, thus effectively enhancing the two dimensionality of the system. Guided by analogy with the cuprates, the authors of ref were lead to suggest that this system maybe understood in terms of doped AFs on TLs, which may well be related to RVB states.

As much as such directions are tempting, one must proceed with due caution since the experimental situation is rather complex. As discussed shortly, the system is actually away from the vicinity of half-filling, the regime where RVB physics is considered to be valid. It is reasonable then to resort to generic methods no more specialized than a combination of symmetry analysis and simple “fermiology”, reinforced by experimental inputs. While experiments conducted on NaxCo2 $`y`$ H2O are at present basically limited to resistivity and preliminary susceptibility data, we have noticed that several conclusions -which are in fact remote from what an RVB-based investigation would lead to- can already be drawn. We therefore believe it useful at this stage to raise notice on this point with the hope of motivating further activity. The seminal work of Rice and Sigrist, who adopted a similar approach in the case of Sr2RuO4, had played an instrumental role in the subsequent identification of the triplet superconductivity in that compound.

During the final stage of this work, we came across a preprint by Baskaran addressing the same system. Although his reduction of the problem to a one-band model is similar, the discussion there is based on the RVB picture. The present work is founded on a complimentary approach as explained above.

Deduction of effective theory

The first point to make is that the present system cannot be taken as a realization of the dilutely doped AF on a TL. The 3d cobalt ions in the parent compound NaCo2O4 with a formal valence of +3.5, are actually a mixture of the low-spin states Co3+(t$`{}_{}{}^{6}{}_{2\mathrm{g}}{}^{}`$e$`{}_{}{}^{0}{}_{\mathrm{g}}{}^{}`$ with S=0) and Co4+(t$`{}_{}{}^{5}{}_{2\mathrm{g}}{}^{}`$e$`{}_{}{}^{0}{}_{\mathrm{g}}{}^{}`$ with S=1/2). With no experimental indications of charge ordering in NaCo2O4, the most natural starting point for studying these systems is as a random mixture (1:1) of S=0 and S=1/2 states distributed on the TL. Hence instead of starting near half-filling, we must consider the case of three-quarter filling, which is a totally new situation. Superconductivity occurs after oxidation of the sodium ion so that the compound becomes Na0.35CoO2 $`y`$ H2O, in which case the system is approximately 1/3 filled. We therfore stress that RVB physics, presumably valid close to half-filling, by no means has a trivial extrapolation to this regime.

The next issue to consider is the magnetic properties of the system. The susceptibility of NaCo2O4 contains a Curie-Weiss component with large negative Curie temperature T$`{}_{\mathrm{Curie}}{}^{}285K`$, a trend also confirmed in measurements on Na0.35CoO2 $`y`$ H2O, though with a reduced value of T$`{}_{\mathrm{Curie}}{}^{}18.8K`$. Taken together with the finding of spin density wave (SDW) like tendencies upon Co to Cu subsitution in NaCo2O4, it is natural to assume a fair amount of AF fluctuation to be present. LDA calculations show that the Fermi surface for NaCo2O4 sits near the top of the a1g like state, split off from the t2g manifold due to oxygen distortions. These circumstances suggest a one-band Hubbard model ($`U58eV,t1eV`$) on a TL around quarter to third filling as the simplest model capturing the essential physics of this system. At this point one may already expect nontrivial pairing tendencies. For instance, there are perturbative studies that claim that triplet pairing is favored over singlet pairing. With direct quantum Monte-Carlo methods unavailable due to severe negative signs, we proceed with the generic approach mentioned earlier.

Symmetry considerations

Group theoretical studies, valid irrespective of the pairing mechanism, have played essential roles in classifying superconductivity in heavy fermions, organics, and e.g. the ruthenate compound. For the present case, the relevant symmetry group is $`𝒢=D_{6h}U(1)T`$, with $`D_{6h}`$ the hexagonal group together with the inversion symmetry about the plane, $`U(1)`$ the gauge symmetry broken by the onset of superconductivity, and $`T`$ the time reversal symmetry. For Na0.35CoO2 $`y`$ H2O, magnetization measurements suggest the presence of a magnetic anisotropy. Hence the absence of $`SU(2)`$ from $`𝒢`$. The irreducible representation for this situation can be classified into 12 classes, $`\mathrm{\Gamma }_a^\pm ,1a6`$, with +(-) standing for spin singlet (triplet) states. States with subscripts $`1a4`$, and those with $`a=5,6`$ are 1d and 2d representations respectively. The conventional point of view that AF fluctuations are most compatible with singlet pairing will lead us to one of these (+) representations.

A more intriguing possibility however is the case of triplet pairing. Not easily accessed from RVB treatments which deal primarily with spin singlets and their fluctuations, it is perhaps worthwhile to highlight their exotic properties. Later on we will mention that experiments so far seem to spell out such a state as the one most probable, and a novel scenario from which we recover this pairing choice will be described. As mentioned there, we will be interested in states with the d-vector perpendicular to the cobalt plane, i.e. $`\stackrel{}{d}//\widehat{z}`$, which ensures that the electron spins contributing to the wavefunctions have a vanishing z-component. It turns out that the $`\mathrm{\Gamma }_5^{}`$ representation realizes such states. The fourth order terms of a Ginzburg Landau expansion of the free energy is of the form

$$_{quad}=d^2\stackrel{}{r}[\beta _1(|\eta _1|^2+|\eta _2|^2)^2+\beta _2(\eta _1^{}\eta _2\eta _1^{}\eta _2)^2],$$
(1)
where $`\eta _1`$ and $`\eta _2`$ are coefficients of the two basis functions $`\widehat{z}k_x`$ and $`\widehat{z}k_y`$. When the coefficient $`\beta _2>0`$, we have $`\stackrel{}{d}=\widehat{z}(k_x+ik_y)`$, a state violating $`T`$. The same pairing state also arises in the tetragonal ($`D_{4h}`$) case, though additional conditions on the GL coefficients are required there. The $`T`$-violation is most dramatically manifested in electric/magnetic properties. The derivative terms of the GL theory for representations $`\mathrm{\Gamma }_5^\pm `$ are

$`={\displaystyle d}`$ $`\stackrel{}{r}[K_1|D_x\eta _1+D_y\eta _2|^2+K_2|D_x\eta _2D_y\eta _1|^2`$ (3)
$`+K_3(|D_x\eta _1D_y\eta _2|^2+|D_x\eta _2+D_y\eta _1|^2)]`$
with $`D_\mu =_\mu 2ieA_\mu `$ the covariant derivatives. All z-derivatives are omitted on account of the two-dimensionality. Adapting a method due to Furusaki et al who discussed the tetragonal case $`D_{4h}`$, we find that the terms in the GL Lagrangian density coupling linearly to the electric field $`\stackrel{}{E}`$ read

$`_E`$ $`=`$ $`c_1(E_x\eta _1^{}+E_y\eta _2^{})(D_x\eta _1+D_y\eta _2)+c.c.`$ (4)
$`+`$ $`c_2(E_x\eta _2^{}E_y\eta _1^{})(D_x\eta _2+D_y\eta _1)+c.c.`$ (5)
$`+`$ $`c_3\{(E_x\eta _1^{}E_y\eta _2^{})(D_x\eta _1D_y\eta _2)`$ (7)
$`+(E_x\eta _2^{}+E_y\eta _1^{})(D_x\eta _2+D_y\eta _1)+c.c.\}`$
where the coefficients $`c_j`$ (1$``$j$``$3) are proportional to the $`K_j`$’s in the free energy. From these one extracts the following contribution;

$`_{\mathrm{CS}\mathrm{like}}`$ $`=`$ $`i(c_3c_1c_2)(\eta _1^{}\eta _2\eta _2^{}\eta _1)`$ (9)
$`\times (A_y_xA_0A_x_yA_0),`$
which is of the Chern-Simons (CS) form. Several consequences follow, which are basically similar to what has been discussed in the context of SrRu2O4. First, a zero field Hall effect is expected, due to a spontaneous (orbital) magnetization $`\stackrel{}{d}(\stackrel{}{k}\times _\stackrel{}{k})\stackrel{}{d}\stackrel{}{\eta }^{}\times \stackrel{}{\eta }`$, where $`\stackrel{}{\eta }(\eta _1,\eta _2)`$. Futhermore one sees from the nonuniversality of the coefficient of the CS term that the Hall conductivity $`\sigma _{xy}`$ is not quantized. The physical reason for the latter feature can be identified microscopically. (Readers are referred to refs. for the tetragonal case.) Inserting the d-vector $`\stackrel{}{d}=\mathrm{\Delta }_0(f_1(\stackrel{}{k})+if_2(\stackrel{}{k}))\widehat{z}`$, with $`f_1(\stackrel{}{k})=2\mathrm{sin}(\frac{\sqrt{3}}{2}k_x)\mathrm{cos}(\frac{1}{2}k_y)`$ and $`f_2(\stackrel{}{k})=2\mathrm{cos}(\frac{\sqrt{3}}{2}k_x)\mathrm{sin}(\frac{1}{2}k_y)+2\mathrm{sin}k_y`$ into the Bogoluibov-de Gennes Hamiltonian

$$H=\mathrm{\Psi }^{}[(ϵ(\stackrel{}{k})\mu )\tau _3\mathrm{𝟏}+\tau _1i(\stackrel{}{d}(\stackrel{}{k})\stackrel{}{\tau })\tau _2]\mathrm{\Psi },$$
(10)
where $`ϵ(\stackrel{}{k})=4\mathrm{cos}(\frac{\sqrt{3}}{2}k_x)\mathrm{cos}(\frac{1}{2}k_y)2\mathrm{cos}(k_y)`$ and performing a gradient expansion to one-loop order, one sees that the correction to the universal contribution $`\sigma _{xy}=e^2/4\pi `$ comes from the normalization of the generating functional due to the Goldstone mode associated with the breaking of the gauge symmetry. This simply reflects the fact that the charge is not a well-defined quantum number in a BCS condensate. On the other hand, the residual U(1) degree of freedom left of the SU(2) spin rotational symmetry suggests that the transport of $`S_z`$ has universal feautures. This expectation is verified by gauging the system via the coupling to a spin-gauge field, from which one confirms that the spin Hall conductivity $`\sigma _{xy}^s`$ is quantised in integer multiples of $`1/2\pi \mathrm{}`$. The integer factor is a topological invariant related to the chirality of the superconducting state, characterized in terms of the quantity $`𝐠={}_{}{}^{t}(\mathrm{Re}d_z\tau _1,\mathrm{Im}d_z\tau _1,ϵ\mathrm{𝟏})`$ as

$$N_{chiral}=\frac{d^2\stackrel{}{k}}{(2\pi )^2}\frac{\mathrm{tr}\left(𝐠_{k_x}𝐠\times _{k_y}𝐠\right)}{\mathrm{tr}|𝐠|^{\frac{3}{2}}},$$
(11)
which is +1 for the present choice of d-vector (-1 for the opposite chirality). To complete the discussion on symmetries and exotic topological properties, we note that vortices obeying nonabelian statistics can be realized in this chiral p-wave state .

Relation to experimental inputs

We list the salient features of the inputs provided by experiments conducted so far.

1. The normal state resistivity as a function of temperature fits rather well to a $`T^2`$-dependence.
2. Preliminary NMR measurement of $`T_1^1`$ at the cobalt site shows a prominent peak just below Tc, reminiscent of the Hebel-Schlichter peak for conventional s-wave superconductors.
3. Spins direct their moment in-plane under an applied magnetic field $``$ 8T.

Deferring discussions on feature 3 for the moment, we focus on the first two. We take feature 1 to be an indication that the system lies in the weak coupling regime. It then becomes plausible to interpret feature 2 in terms of a fully-gapped pairing state without nodes. Although this is usually considered a hallmark of s-wave pairing, it is worth recalling that the Balian-Werthamer state, with $`\stackrel{}{d}_{BW}=\mathrm{\Delta }_0\stackrel{}{k}`$ has a sufficiently strong singularity of the density of states to exhibit a coherence peak as well. The chiral p-wave state (with the cylindrical Fermi surface) discussed above is the two-dimensional analogue of this situation. An in-plane measurement of the phase-sensitive thermal conductivity on single crystals, when they become available, is expected to give an isotropic result. (The absence of a peak in the case of Sr2RuO4, which had also been considered to be a chiral-p wave superconductor, is consistent with the existence of horizontal nodes related to interlayer pairing.) States such as $`\stackrel{}{d}=p_x\widehat{z}`$ and $`\stackrel{}{d}=p_y\widehat{z}`$, permissable on symmetry grounds, can be dismissed on bases of these observations. Alternative fully-gapped states include the singlet $`d_{x^2y^2}+id_{xy}`$ and $`d_{x^2y^2}+is`$ states. One can directly distinguish between singlet and triplet states by performing a Knight-shift measurement. For singlet pairings, the shift $`K`$ will drop to zero on approaching $`T=0`$ for all magnetic field directions, with $`KY(\stackrel{}{k},T)e^{\frac{\mathrm{\Delta }_0}{T}}`$, where $`Y(\stackrel{}{k},T)`$ is the Yoshida function. Triplets on the other hand will exhibit an anisotropic magnetic susceptibility in the superconducting state; an in-plane field $`\stackrel{}{H}\widehat{z}`$ will polarize both the Cooper pairs and the quasiparticles, so that $`K_{//}(T)=\mathrm{const}.`$ is expected, whereas the shift for $`\stackrel{}{H}//\stackrel{}{z}`$ will lead to a behavior similar to the singlet case, $`K_{}Y(\stackrel{}{k},T)`$. Remarkably, preliminary Knight shift measurements do show a temperature-independent behavior. (A large $`H_{c2}61`$T also is suggestive of triplet pairing.) These results, taken at face value single out the nodeless triplet pairing, of which the chiral-p state is the simplest.

Possible scenario for triplet pairing in Na0.35CoO2 $`y`$ H2O

According to LDA calculations, the central cylindrical Fermi surface of NaCo2O4, with dominant $`a_{1g}`$ character bears a shape close to a hexagon when viewed in the $`k_z=\mathrm{const}.`$ plane. Hence there are three (approximate) nesting vectors, given by $`𝐐_i=\frac{4}{5}\mathrm{\Gamma }K_i,i=1,2,3`$. Here we define the $`K`$-points as $`K_i=C_6^{i1}(\frac{1}{3}𝐆_1+\frac{1}{3}𝐆_2)`$, $`𝐆_1=2\pi \frac{𝐛\times 𝐜}{𝐚𝐛\times 𝐜}`$, $`𝐆_2=2\pi \frac{𝐜\times 𝐚}{𝐚𝐛\times 𝐜}`$, in which $`𝐚=(1,0,0),𝐛=(1/2,\sqrt{3}/2,0),𝐜=(0,0,1)`$ are the primitive vectors of the layered TL. $`C_6`$ denotes the anticlockwise rotation about the z-axis by an angle of $`\pi /6`$. The presence of such nesting instabilities in the $`\mathrm{\Gamma }K`$ direction has been suggested by Terasaki et al based on their observation of an order-from-disorder type emergence of a spin density wave (SDW) in Cu-substituted NaCo1.8Cu0.2O4. We are assuming below that these qualitative Fermi surface features remain intact for Na0.35CoO2 $`y`$ H2O, which appears to be in accord with the observed Curie-Weiss contribution to the susceptibility.

Fermiology, the study of many-body effects with an emphasis on Fermi surface properties has been playing a central role in advancing our knowledge on superconductivity in various correlated electron systems. One such scenario, which takes advantage of nesting tendencies and leading to triplet pairing, has been advanced for the tetragonal lattice case, with the aim of explaining the pairing in Sr2RuO4. We would like to point out that a very natural extension to the present TL case is possible. This argument starts by recasting of the linearized BCS gap equation at Tc into the form, $`T_ce^{\frac{1}{N(O)<<V_\varphi >>_{FS}}}`$, where we are following the notations of ref.: $`<<V_\varphi >>_{FS}`$ is the pairing interaction $`V(\stackrel{}{k}\stackrel{}{k}^{})`$ averaged over the Fermi surface, i.e.

$$<<V_\varphi >>_{FS}=\frac{_{FS}𝑑\stackrel{}{k}_{FS}𝑑\stackrel{}{k}^{}V_\varphi (\stackrel{}{k}\stackrel{}{k}^{})\varphi (\stackrel{}{k})\varphi (\stackrel{}{k}^{})}{[_{FS}𝑑\stackrel{}{k}^{}]_{FS}𝑑\stackrel{}{k}\varphi ^2(\stackrel{}{k})},$$
(12)
with $`\varphi (\stackrel{}{k})`$ the $`\stackrel{}{k}`$-dependent part of the order parameter, i.e. $`\mathrm{\Delta }(\stackrel{}{k})=\mathrm{\Delta }_0\varphi (\stackrel{}{k})`$. We may view this as a sort of variational principle in which the actual pairing occurs for cases which yield a (1)positive and (2)large $`<<V_\varphi >>_{FS}`$. In the present situation, where we are dealing with AF fluctuations, the pairing potential $`V(𝐐_i)>0`$. For this to lead to pairing requires that $`\varphi (\stackrel{}{k})\varphi (\stackrel{}{k}+𝐐_i)<0`$ for some $`i`$. The simplest way to realize this is to introduce a nodal direction which coincides e.g. with the $`k_x`$ axis, i.e. a $`p_y`$-pairing state. One immediately sees that this choice is highly beneficial as all three nesting vectors satisfy the required condition (see Fig. 1(a)), and therefore the entire Fermi surface is available for the pairing. (For this $`p_y`$-pairing, the $`𝐐_1`$-channel should be the most dominant because it coincides with the direction in which the magnitude of the gap becomes maximal.) In this respect, the advantage of such node formation (vertical nodes on the central Fermi surface) is bigger than in the ruthenates where only one pair of nested segments is involved in the pairing. This state is degenerate (i.e. has the same Tc) with those in which the nodes are located at the intersections between the hexagon Fermi surface and the lines $`k_y=\pm \sqrt{3}k_x`$ (Fig. 1(b) and (c)). Further gain in condensation energy is achieved by constructing a linear combination of these three states in such a way that the nodes will cancel one other, opening a full gap at all points on the Fermi surface. This becomes possible through the construction

$$\varphi (\stackrel{}{k})_{p_y+ip_x}\varphi _{p_y}(\stackrel{}{k})+i\left(\varphi _{p_y}(C_6\stackrel{}{k})+\varphi _{p_y}(C_6^2\stackrel{}{k})\right).$$
(13)
As will be discussed elsewhere, one can check that a similar construction cannot be made for an f-wave state. The appearance of the chiral p-wave state above relied on the Fermi surface geometry (reflecting the hexagonal symmetry), and the assumption of a nesting tendency. It is desirable that fine quality samples, amenable to detailed neutron scattering measurements be made in order to make a direct verification on this point possible. Further support for triplet pairing comes from the aformentioned magnetic anisotropy; when we interpert this in terms of an easy-axis along the c-axis, evaluation along the lines of ref. can be used to demonstrate that the anisotropy ratio of the in-plane and longitudinal susceptibilities will directly enter into expressions for the triplet and singlet pairing potentials, and a sufficient anisotropy will eventually favor triplets (with $`\stackrel{}{d}//\widehat{z}`$) over singlets.

We are indebted to H. Sakurai, K. Takada, T. Sasaki, M. Arai and E. Takayama-Muromachi for access to their results prior to publication, and for extensive discussions. We would also like to acknowledge K. Kindo, K. Yoshimura and their group for generously allowing us to mention their results on magnetization and NMR. The similarity of the present system with the ruthenates was pointed out to the authors by I. Terasaki. We have also meritted from discussions with T. Hikihara, T. Miyazaki, T. Sasaki, Y. Tateyama, K. Izawa, Y. Matsuda and J. Goryo.