text:
# References

Remarks on: ‘Theory of Neutrino Oscillations’ (hep-ph/0311241) by C.Giunti, the comments by L.B.Okun et al in hep-ph/0312151, and Giunti’s reply in hep-ph/0312180

J.H.Field

Département de Physique Nucléaire et Corpusculaire, Université de Genève

24, quai Ernest-Ansermet CH-1211Genève 4.

E-mail: john.field@cern.ch

## Abstract

As pointed out previously, some of the hypotheses used in the standard discussion of the quantum mechanics neutrino oscillations, recently summarised in the paper ‘Theory of Neutrino Oscillations’ by C.Giunti, lack any physical foundation and lead to neglect of both a factor of two correction in the contribution of neutrino propagation to the oscillation phase and an important contribution to the latter from the propagator of the source particle.

PACS 03.65.Bz, 14.60.Pq, 14.60.Lm, 13.20.Cz
Keywords ; Quantum Mechanics, Neutrino Oscillations.

C.Giunti has performed a valuable service in making a clear and concise statement in Ref.  (TNO) of the basic physical assumptions underlying the standard quantum mechanical treatment of neutrino oscillations. Of particular importance are the assumptions A2 and A3 of TNO which are the essential ones in order to derive the standard result for the oscillation phase. However, A2, as written, is either factually incorrect, or, on the most generous interpretation highly likely to be misinterpreted. Indeed, it is the incorrect interpretation that is used in the standard derivation. In fact the initial and final states in the amplitudes of any standard model process are, uniquely, the mass eigenstates $`|\nu _k`$ not the ‘neutrino flavour eigenstate’ $`|\nu _\alpha `$ defined in Eqn(1) of TNO. Since the MNS matrix element, $`U_{\alpha k}`$, necessarily appears in the amplitude when the state $`|\nu _k`$ is created or destroyed, in association with a charged lepton with label $`\alpha =e,\mu ,\tau `$, one can introduce, formally, the state $`|\nu _\alpha `$ in order to write, in a compact way, a general expression for the charged weak current or the corresponding Lagrangian, always bearing in mind, however, that it has no actual physical significance, as the state does not appear in the amplitude of any standard model process. Now what is done in the standard derivation is to assume that this formal, physically non-existent state is actually created by the neutrino production process at some fixed time. Since the detection process evidently occurs at another fixed time, the implicit assumption is made that all the mass eigenstates have the same velocity, which is just what is explicitly stated in A3, so that two assumptions are certainly consistent with each other. Two remarks follow. Firstly the state $`|\nu _\alpha `$ is  not produced, because there is no physical reason why it should be. As already stated above only the mass eigenstates $`|\nu _k`$ appear in physical amplitudes as initial or final states.. As I have shown in a recent paper , assuming that such states are created in pion decay gives a prediction for the ratio $`\mathrm{\Gamma }(\pi e\nu )/\mathrm{\Gamma }(\pi \mu \nu )`$ completely excluded by experiment. Interestingly enough, Giunti has himself recently given arguments based on quantum field theory why such states are not produced in pion decay , wthout, apparently, realising the contradiction with the assumption A2 of TNO! Secondly, which is essentially what is pointed out in the comment of Okun et al , the assumptions A2 and A3 of TNO are in clear contradiction with A4. On the one hand, it is assumed (A2 and A3) that the space-time velocities, $`L/T`$, are the same for all mass eigenstates, on the other hand (A4) they must have different kinematical velocities, $`p/E_k`$. Aristotle made a wise remark that seems quite appropriate to this situationaaa On Giunti’s website ‘Neutrino Unbound’: http://www.nu.to.infn.it/ a fair number of wise quotations are to be found, even one in connection with a paper of mine. I therefore consider that I have the right to include at least one here

In inventing a model we may assume what we wish, but should avoid
impossiblities.

Since, however, the different mass eigenstates do not have to be produced at the same time in the different path amplitudes whose interference results in the ‘neutrino oscillation’ phenomenon, the hypotheses A2 and A3 are generally incorrect so that there in no contradiction with A4. I am, of course, using here Feynman’s space-time formulation of QM which is the one best adapted to flavour oscillation problems, for the theoretical description of the problem.

Thus Assumptions A2 and A3 are usually false. A coherent state $`|\nu _\alpha `$ is  not produced , and although the source-detector distance, $`L`$, is the same for all states $`|\nu _k`$, the times-of-flight, $`T`$, (and hence the times of decay of the source) are different in the different amplitudes: $`T_k/L=p_k/E_k`$. Taking into account the different decay times of the source increases the contribution of the neutrino propagators to the interference phase by a factor of two and results in a numerically important contibution to the interference phase from the propagator of the source particle. There is therefore in fact no incompatiblity between the space-time and kinematical velocities, and no need, as asserted by Giunti in TNO and Ref. , to introduce an ad hoc gaussian ‘spatial wave-packet’ to circumvent the problem.

There should be nothing shocking in the idea that a source can decay at different (but unobserved) times in quantum interference phenomena. Consider a Michelson interferometer with a 1m path difference between the two arms. The times of decay of the source atom that produces the photon that ‘interferes with itself’ must obviously differ by 2m/c = 6.7 nsec when interference fringes are observed.

Unfortunately, practitioners of the theory of neutrino oscillations appear to have little familiarity with Feynman’s space-time approach to QM, being instead usually experts in text-book second-quantised field theory, which has limited relevance to space-time flavour oscillations. This is not the case, however, for atomic physics theorists. In concluding below I will briefly describe an atomic physics experiment that is an almost perfect analogue of the two-flavour neutrino oscillation problem. Feynman’s path amplitude method has been used to predict quantum interference effects in the experiment in very good agreement with observation.

One has the impression that neutrino physics theorists are stuck, for their conceptial understanding of QM, with Heisenberg in 1930. It is as though Dirac’s  and Feynman’s  discoveries of the true role of space-time in QM were never made. At this very early period in the history of QM, statements like the following can often be found :

There exists a body of exact mathematical laws but these cannot
be interpreted as expressing simple relationships between objects
in space and time.

Or, :

EITHER:

* Causal relationship expressed by mathematical laws
* BUT: physical description of space-time impossible

OR

* Phenomena described in terms of space and time
* BUT: Uncertainty Principle

The Feynman path amplitude formulation actually gives a mathematical law expressing a causal relationship (i.e. ordered in time) between elementary processes in space-time! The above quotations a very reminiscent of the woolly arguments invoked in blanket fashion by Giunti to avoid thinking clearly about the space time aspects of the QM of neutrino oscillations.

For example, misuse of the Heisenberg Uncertanity relations :

The wave packet treatment of neutrino oscillations is also necessary for a correct description of the momentum and energy uncertainties necessary for the coherent production and detection of different massive neutrinos  whose interference generates the oscillations

It is not the wave functions of the neutrinos that are ‘coherent in the oscillation phenomenon, but rather the path amplitudes for the whole experiment. The neutrinos are only unobserved intermediate states in these amplitudes.

or :

One is free to define ‘so-called space velocities $`\overline{v}=x/t`$ of massive
neutrinos which are trivially identical. However I think that such a
definition is useless. Neutrinos are not classical objects for which
$`v=x/t`$. It is pretty clear that the uncertainty principle forbids any
relation of this type

In fact, according to quantum field theory, on-shell particles do propagate classically in space-time; only highly virtual ones show different behaviour. Notice the blanket invocation of ‘Heisenberg uncertainty’ to block all further discussion. Actually the relation $`v=x/t`$ is absolutely essential for the correcct quantum mechanical description of neutrino oscillations.

and

Of course energy and momentum are not exactly defined in the production process, in order to allow the coherent production of different massive
neutrinos, but it is pretty unlikely that the uncertainty in energy and
momentum could generate the equal energy constraint.

In relativistic quantum field theory energy and momentum are both exactly conserved at all vertices. ‘Heisenberg uncertainty’ (see below) only appears as a smearing of the masses of virtual particles about pole values. Another, desperate, last resort invocation of the ‘Uncertainty Principle’!

As discussed in detail in my paper Ref. , the misuse of spatial wavepackets by Giunti and many other authors is related to a confusion of QM with classical wave theory, for which the concepts of phase velocity, group velocity and wave packets are useful and meaningful concepts. In the Dirac-Feynman formulation of QM they become almost irrelevant. Still, the Heisenberg Uncertainty Principle is an important part, when correctly used, of QM. I describe below its application to neutrino production in pion decay. In fact the neutrinos in $`\pi \mu \nu `$ (not in $`\pi e\nu `$) are described by a momentum wave packet. It has however nothing to do with the Fourier transform of the arbitary gaussian spatial function that is parachuted into the theory by Giunti and many other authors.

It is a consequence of quantum field theory that on-shell particles, or virtual particles propagating over macroscopic spatial separations, do so in a classical manner . The dominant trajectories in Feynman path integrals are the classical ones. Consider the space-time description of pion decay, $`\pi \mu \nu `$, in QM. The Heisenberg Uncertainty Principle (so dear to the hearts of neutrino physics theorists!) should be correctly invoked in three different ways in the description of this process: (i) the decay width $`\mathrm{\Gamma }`$ and the decay lifetime $`\tau `$ certainly respect the energy-time Uncertainty Relation $`\mathrm{\Gamma }\tau =1`$; (ii) as a result of the finite pion lifetime the physical mass of the pion differs from the pole mass. In accordance with the energy-time Uncertainty Relation the mass smearing, given by a Breit-Wigner amplitude, is inversely proportional to the decay lifetime; (iii) due to the finite muon decay lifetime, the physical mass of the muon is also smeared around its pole value. Since the muon is unobserved, overall energy-momentum conservation leads to a neutrino momentum wave-packet in the case of $`\pi \mu \nu `$ but not for $`\pi e\nu `$, since the electron has an infinite lifetime. The effects of this wave packet are calculated in my paper Ref. . They are tiny. The upshot is, that once the neutrinos are produced (respecting, of course, the law $`\mathrm{\Gamma }\tau =1`$) they propagate in space-time as essentially classical particles. This wave packet is the only one that is relevant to pion decay. In particular there is no ad hoc gaussian spatial wave packet describing the source. The uncertainty in the position of the source particle, a property of the initial state, and therefore the same in all path amplitudes, gives only a small, incoherent, correction to the oscillation phase. Again, this effect is calculated in Ref. and found to be very small.

I would like to remark that I do agree with Giunti that the ‘equal energy’ hypothesis proposed in Ref.  is not relevant to resolving the evident contradiction between assumptions A2 and A3 of TNO as compared to A4. As shown in Ref.  sufficient assumptions to obtain the standard oscillation phase are A2 or A3 (since A2 implies A3). The result at, O($`m_\nu ^2`$), is independent of the kinematical assumptions such as A4 (equal momenta, equal energies or exact energy-momentum conservation). In fact, in the papers by Stodolsky  and Lipkin  cited ini Ref. , not onlt is the equal energy assumption made, but the temporal part of the neutrino propagator phase (necessarily non-zero in the laboratory system by Lorentz invariance) is, arbitarily and incorrectly, set to zero.

A misnomer which occurs widely in Refs.  is ‘plane wave’ for ‘space-time propagator’. This function has meaning only as a factor in a path amplitude. Interpreted as a wavefunction it is non square-integrable and therefore devoid of any physical significance.

The atomic physics experiment that provides a close analogy with a two flavour neutrino oscillation problem is called the ‘Photodetachment Microscope .A coherent source of electrons of fixed energy is provided by a negative ion beam irradiated by a laser. The detached electron moves in a constant electric field before detection. Just two classical trajectories link the point of emission to any point on a plane detector oriented perpendicularly to the electric field direction. Quantum interference effects are observed between the path amplitudes corresponding to the two trajectories. A good pedagogical description can be found in Ref.  where the appropriate path integral formulabbb A similar formula has recently been proposed for the neutrino oscillation problem by Paz̆ma and Vanko . The corresponding oscillation phase was not, however, derived.r :

$$\psi (\stackrel{}{r},t_f)=_{\mathrm{}}^{t_f}\mathrm{exp}[i\frac{ϵt_i}{\mathrm{}}]\mathrm{exp}[i\frac{S_{cl}(\stackrel{}{r},t_i,t_f)}{\mathrm{}}]𝑑t_i$$

is given.

In this formula $`ϵ`$ is the energy of the detached electron and $`S_{cl}`$ the classical action corresponding to an electron trajectory. Note particularly the time integral on the RHS of the equation. The first exponential function is the propagator of the coherent source (analagous to that of a coherent neutrino source) the second represents the propagator of the electron in the electric field. In practice it is well approximated by the contributions of the two classical trajectories mentioned above, corresponding to values of $`t_i`$ with a fixed separation. These are the analogues of the propagators of different neutrino mass eigenstates. A typical value of the difference in $`t_i`$ between the two trajectories, quoted in Ref.  is 160 psec for a time-of-flight of 117 nsec. The laws of physics are the same in this and in any neutrino oscillation experiment. Contrary to the claims in TNO and Ref. \], no ad hoc ‘wave packets’ are required for a correct quantum mechanical description of such systems. The standard formula for the oscillation phase which, as discussed in Ref.  may well be correct for heavy quark oscillations, is then in most cases, not so  for neutrino oscillations.

Finally, it is remarked that as described in Ref.  a clear experimental discrimination between the standard and path amplitude formulae of the oscillation phase is provided by long-linebase ‘$`\nu _\mu `$ disappearence’ experiments with pion and kaon source beams.


format: text

text:
# Pacman (I): A new Algorithm to calculate Faraday Rotation Maps

## 1 Introduction

Magnetic fields are a common phenomenon on different astrophysical scales. One possibility to observe astrophysical magnetic fields and thus, gain insight into their structure, ordering scale and field strength is given by the Faraday rotation effect. It arises whenever linearly polarised radio emission passes through a magnetised medium. The polarisation angle $`\phi `$ of the radiation will be rotated due to the interaction between the magnetic field component parallel to the propagation direction of radiation and the radio emission itself. The resulting change in polarisation angle is proportional to the wavelength squared $`\lambda ^2`$ which can be expressed by $`\phi =\mathrm{RM}\lambda ^2+\phi ^0`$, where the proportionality constant RM in this relation is called Faraday rotation measure and $`\phi ^0`$ is the intrinsic polarisation angle of the emission at the radio source.

The latter two quantities, RM and $`\phi ^0`$, can be measured through multi-frequency polarisation observations. The RM can be described as the line of sight integral over the electron density and the magnetic field component parallel to the line of sight. Therefore, the RM contains valuable information on the magnetic field in the foreground of the polarised radio source. The intrinsic polarisation angle $`\phi ^0`$ and thus, the intrinsic magnetic field direction at the radio source gives insight into the magnetic field configuration at the source.

The statistical analysis of RM measurements in terms of correlation functions and equivalently power spectra as developed by Vogt & Enßlin (2003) and Enßlin & Vogt (2003) requires that the RMs are unambiguously determined. Thus, any ambiguous RM can lead to misinterpretation of the data investigated. For the calculation of RM and $`\phi ^0`$ using the relationship $`\phi =\mathrm{RM}\lambda ^2+\phi ^0`$ mentioned above, a least squares fit is normally applied to the polarisation angle data. Since the measured polarisation angle $`\phi `$ is constrained only to values between 0 and $`\pi `$ leaving the freedom of additions of $`\pm n\pi `$, where $`n`$ is an integer, the determination of RM and $`\phi ^0`$ is ambiguous, causing the so called $`n\pi `$-ambiguity. Therefore, a least squares fit has to be applied to all possible $`n\pi `$-combinations of the polarisation angle data at each data point of the polarised radio source while searching for the $`n\pi `$-combination for which $`\chi ^2`$ is minimal.

In principle, $`\chi ^2`$ can be decreased to infinitely small numbers by increasing RM substantially. Vallée & Kronberg (1975) and Haves (1975) suggested to avoid this problem by introducing an artificial upper limit for $`|\mathrm{RM}|\mathrm{RM}_{\mathrm{max}}`$. Since this is a biased approach, Ruzmaikin & Sokoloff (1979) proposed to assume that no $`n\pi `$-ambiguity exists between the measurements of two closely spaced wavelengths taken from a whole wavelength data set. The standard error of the polarisation measurements is then used to constrain the possible $`n\pi `$-combinations for the least squares fit for the subsequent observed frequencies. We refer to these methods as the “standard fit” algorithms, as they are currently the most widely used methods. However, these methods might still give multiple acceptable solutions for data with low signal-to-noise, requiring the solution to be flagged and all the information carried by these data is lost. Furthermore, it still can happen that the algorithm chooses a wrong RM and imprints spurious artefacts on the RM and the $`\phi ^0`$ maps.

Recently, a completely different approach was proposed by Sarala & Jain (2001) which takes the circular nature of the polarisation angle into account. The authors apply a maximum likelihood method to spectral polarisation data. Although this approach is not biased towards any RM value, it is rather designed for a large number of observed wavelengths. Similarly, de Bruyn (1996) and Brentjens & de Bruyn (2004) propose a RM-synthesis via wide-band low-frequency polarimetry. However, typically the observations are only performed at three or four wavelengths especially for extended (diffuse) radio sources.

Here, we propose a new approach for the unambiguous determination of RM and $`\phi ^0`$. We assume that if regions exhibit small polarisation angle gradients between neighbouring pixels in all observed frequencies simultaneously, then these pixels can be considered as connected. Note, that we calculate the gradient modulo $`\pi `$, which implies that polarisation angles of 0 and $`\pi `$ are regarded as having the same orientation and thus, the cyclic nature of polarisation angles is reflected. Information about one pixel can be used for neighbouring ones and especially the solution to the $`n\pi `$-ambiguity should be the same.

In cases of small gradients, assuming continuity in polarisation angles allows us to assign an absolute polarisation angle for each pixel with respect to a reference pixel within each observed frequency. This assignment process has to be done for each spatially independent patch of polarisation data separately, such as each side of a double-lobed radio source. The reference pixel is defined to have a unique absolute polarisation angle and the algorithm will start from this pixel to assign absolute polarisation angles with respect to the reference pixel while going from one pixel to its neighbours. Figuratively, the algorithm eats its way through the set of available data pixel. It might become clear now why we have chosen to call the algorithm *Pacman*111The computer code for Pacman is publicly available at http://dipastro.pd.astro.it/~cosmo/Pacman (Polarisation Angle Correcting rotation Measure ANalysis).

Pacman reduces the number of least squares fits in order to solve for the $`n\pi `$-ambiguity. Preferably, the reference point is chosen to have a high signal-to-noise ratio so that in many areas, it is sufficient to solve the $`n\pi `$-ambiguity only for a small number of neighbouring pixels simultaneously and to use this solution for all spatially connected pixels. Pixels with low signal-to-noise will profit from their neighbouring pixels allowing a reliable determination of the RM and $`\phi ^0`$.

In Sect. 2, we describe in detail the idea and the implementation of the Pacman algorithm. In Sect. 3, we test Pacman on artificially generated RM maps and demonstrated its ability to solve the $`n\pi `$-ambiguity properly. However, the application to observational polarisation data and statistical characterisation of the resulting maps will be presented in the second paper (Vogt et al.), to which we will refer as Paper II.

## 2 The new Pacman algorithm

### 2.1 The Idea

As described in the introduction, the Faraday rotation measure $`\mathrm{RM}_{ij}`$ at each point with map pixel coordinate $`(ij)`$ of the source, is usually calculated by applying a least squares fit to measured polarisation angles $`\phi _{ij}(k)`$ observed at frequency $`k1\mathrm{}f`$ such that

$$\phi _{ij}(k)=\mathrm{RM}_{ij}\lambda _k^2+\phi _{ij}^0,$$
(1)
where $`\phi _{ij}^0`$ is the intrinsic polarisation angle at the polarised source.

Since every measured polarisation angle is observationally constrained only to a value between 0 and $`\pi `$, one has to replace $`\phi _{ij}(k)`$ in the equation above by $`\stackrel{~}{\phi }_{ij}(k)=\phi _{ij}(k)\pm n_{ij}(k)\pi `$, where $`n_{ij}(k)`$ is an integer, leading to the so called $`n\pi `$-ambiguity. Taking this into account, a least squares fit to calculate RMij and $`\phi _{ij}^0`$ at each pixel has to be applied by allowing all possible combinations of $`n_{ij}(k)\pi `$ while determining the $`n_{ij}(k)`$ for which the $`\chi _{ij}^2`$ is minimal. The presence of observational noise might cause a standard least squares fit, as suggested by Vallée & Kronberg (1975), Haves (1975) or Ruzmaikin & Sokoloff (1979), to choose a spurious RM value especially for areas of low signal-to-noise.

The idea of Pacman is to reduce the number of pixels for which the $`n\pi `$-ambiguity has to be individually solved. This is done by splitting the solution of the $`n\pi `$-problem into two problems, a local and a global one,

$$n_{ij}(k)=\stackrel{~}{n}_{ij}(k)+n(k),$$
(2)
where $`\stackrel{~}{n}_{ij}(k)`$ is the local solution, linking polarisation angles of neighbouring pixels within a frequency map, and $`n(k)`$ is the global solution to the problem, linking polarisation angles of the different frequencies. The local part $`\stackrel{~}{n}_{ij}(k)`$ is determined by construction of absolute polarisation angle maps for each frequency with respect to a high signal-to-noise reference pixel being defined to possess a unique polarisation angle. The term absolute polarisation angle is to be understood as a value determined relative to the reference pixel by adding $`\pm n\pi `$ to the measured polarisation angles in order to remove jumps of the order of $`\pi `$ in the measured polarisation angle map of each observed frequency. The global $`n\pi `$ ambiguity is solved for a high signal-to-noise area surrounding the reference pixel resulting in $`n(k)`$. This is then also the solution of the global $`n\pi `$-problem for all spatially connected points which are assigned absolute polarisation angles with respect to this reference area.

The splitting of the problem in a local and a global one is possible if the real polarisation angle $`\overline{\phi }_{ij}(k)`$ ($`\overline{\phi }_{ij}(k)=\mathrm{RM}_{ij}\lambda _k^2+\phi _{ij}^0`$) is a smooth quantity which does not change more than $`\pm \pi /2`$ between neighbouring pixels. In Sect. 2.2, it is described how Pacman ensures that only pixels fulfilling this condition are used.

The source might consists of several spatially independent areas of polarisation, to which we will refer as patches in the following. The $`n\pi `$-ambiguity has to be solved for each of these patches separately, requiring separate reference pixels to be defined for the construction of absolute polarisation angles.

The advantage of Pacman is that the global $`n\pi `$-ambiguity is solved only for pixels having the highest signal-to-noise ratios. Therefore, noisier pixels which are situated at the margin of the source profit from an already defined global solution $`n(k)`$ to the $`n\pi `$-ambiguity, making a reliable determination of RMij and the intrinsic polarisation angle $`\phi _{ij}^0`$ for these pixels possible.

### 2.2 The Basic Algorithm

The algorithm starts with reading the maps of polarisation angles and of its errors, processes them, and when finished, it saves the different calculated maps. A flow chart of our Pacman algorithm is exhibited in Fig. 1, which shows schematically the procedure the algorithm follows in order to determine the solution to the $`n\pi `$-ambiguity and to calculate the various maps.

After loading of the various polarisation data for the different frequencies, the algorithm calculates a quality measure map. Different actions are involved in the processing of the data. One of the processes is building a patch which involves the construction of absolute polarisation angle maps. First, the reference pixel for a patch of polarisation is determined from the quality measure map. This step is followed by solving the local $`n\pi `$-ambiguity in each frequency map for pixels being spatially connected to the defined patch reference pixel and having small gradients in polarisation angle to the neighbouring pixels at all frequencies. This process is mainly governed and stopped by quality requirements.

When the patch building process is stopped then the algorithm solves the global $`n\pi `$-ambiguity across the frequency maps individually for a selected set of superior quality pixels (voters) within the patch. The majority of votes is used as the global patch $`n\pi `$-solution and the resulting RM and $`\phi ^0`$ values are stored. If there are still unprocessed high quality pixel left, the algorithm starts again to build a new patch. Otherwise Pacman finishes and saves the resulting complete maps.

The values for RM and $`\phi ^0`$ are calculated following a weighted least squares fit expressed by

$`\mathrm{RM}_{ij}`$ $`=`$ $`{\displaystyle \frac{\overline{S_{ij}}\overline{\lambda ^2\phi _{ij}}\overline{\lambda ^2}\overline{\phi _{ij}}}{\overline{S_{ij}}\overline{\lambda ^4}\overline{\lambda ^2}^2}}`$ (3)
$`\phi _{ij}^0`$ $`=`$ $`{\displaystyle \frac{\overline{\phi _{ij}}\overline{\lambda ^4}\overline{\lambda ^2}\overline{\lambda ^2\phi _{ij}}}{\overline{S_{ij}}\overline{\lambda ^4}\overline{\lambda ^2}^2}}`$ (4)
where $`\overline{S_{ij}}=_{k=1}^f1/\sigma _{k_{ij}}^2`$, $`\overline{\phi _{ij}}=_{k=1}^f\phi _{ij}(k)/\sigma _{k_{ij}}^2`$, $`\overline{\lambda ^2}=_{k=1}^f\lambda _k^2/\sigma _{k_{ij}}^2`$, $`\overline{\lambda ^2\phi _{ij}}=_{k=1}^f\lambda _k^2\phi _{ij}(k)/\sigma _{k_{ij}}^2`$ and $`\overline{\lambda ^4}=_{k=1}^f\lambda _k^4/\sigma _{k_{ij}}^2`$. In these relations, $`\sigma _{k_{ij}}`$ is the standard error of the polarisation angle $`\phi _{ij}(k)`$ at the pixel coordinate $`(ij)`$ at the $`k`$th wavelength $`\lambda _k`$.

Since the performance of the algorithm is mainly governed by quality requirements, we need some quality measure in order to rank the map pixels. One good candidate is the expected uncertainty $`\sigma _{ij}^{\mathrm{RM}}`$ of any RM value obtained in an error weighted least squares fit which is calculated by

$$\sigma _{ij}^{\mathrm{RM}}=\sqrt{\frac{\overline{S_{ij}}}{\overline{S_{ij}}\overline{\lambda ^4}\overline{\lambda ^2}^2}},$$
(5)
where the terms are defined as above.

The uncertainty $`\sigma _{ij}^{\mathrm{RM}}`$ of the RM value is used to assign each pixel an initial quality $`q_{ij}=\sigma _{ij}^{\mathrm{RM}}`$ since $`\sigma _{ij}^{\mathrm{RM}}`$ accounts for the statistic used to determine the RM maps. Hence, small values of $`q_{ij}`$ indicate high quality pixels. The quality is then used to determine the way Pacman goes through the data which will be preferably from high to low quality data pixels.

For the construction of absolute polarisation angle maps for each frequency, Pacman starts at the best quality pixel having the smallest value $`q_{ij}`$, which is defined to be the reference pixel. For this first point, the measured polarisation angle $`\phi _{ij}(k)`$ is defined to possess a unique polarisation angle value $`\overline{\phi }_{ij}(k)=\phi _{ij}(k)`$ for each observed frequency. It is important to note that the reference pixel is the same for all frequencies.

Then, Pacman compiles what we call the border list (BL). It contains pixels being direct neighbours to points which have been already assigned an absolute polarisation angle. In the following, the set of adjacent pixels ($`i\pm 1j\pm 1`$) to the pixel ($`ij`$) will be referred to as direct neighbours ($`i^{}j^{}`$). Beginning with the direct neighbours of the reference pixels the border list is continuously updated during the progression of the algorithm. A pixel can be rejected from the border list if the standard error $`\sigma _{k_{ij}}`$ of the polarisation angle for any frequency at this pixel exceeds a certain limit $`\sigma _k^{\mathrm{max}}`$ which is set at the beginning of the calculation. However, this requirement can be relaxed (see Sect. 2.7).

Having defined the reference pixel, Pacman assigns absolute polarisation angles to pixels within the border list always starting with the pixel having the best quality, i.e. the lowest $`q_{ij}`$ value. For this pixel, the algorithm solves the local $`n\pi `$-ambiguity with respect to $`\stackrel{~}{n}_{ij}(k)𝒁`$ by minimising the expression

$$\sigma _{ij}^\mathrm{\Delta }=\underset{\{i^{}j^{}\}D_{ij}}{}\left[\left(\phi _{ij}(k)\pm \stackrel{~}{n}_{ij}(k)\pi \right)\overline{\phi }_{i^{}j^{}}(k)\right]^2,$$
(6)
where $`D_{ij}`$ is the set of all direct neighbours $`(i^{}j^{})`$ to the pixel $`(ij)`$ which have already been assigned an absolute polarisation angle. The resulting $`\stackrel{~}{n}_{ij}(k)`$ value determines the absolute polarisation angle $`\overline{\phi }_{ij}(k)=\phi _{ij}(k)\pm \stackrel{~}{n}_{ij}(k)\pi `$ which has the smallest difference to the already defined absolute polarisation angles of adjacent pixels $`\overline{\phi }_{i^{}j^{}}(k)`$. For each pixel, this is done at each frequency $`k1\mathrm{}.f`$ independently but simultaneously. Thus, Pacman goes the same way through the data in each frequency. Note that one can introduce a value $`\sigma _{\mathrm{max}}^\mathrm{\Delta }`$ which causes Pacman to reject the pixel if $`\sigma _{ij}^\mathrm{\Delta }`$ exceeds $`\sigma _{\mathrm{max}}^\mathrm{\Delta }`$ (see Sect. 2.6). Additionally the border list is updated to include direct neighbours to the recently processed pixel $`(ij)`$ which have not yet been assigned an absolute polarisation angle.

Pacman repeats this process for the best remaining pixel in the border list which has the lowest $`q_{ij}`$ and so on, until the whole patch consists of spatially connected pixels with assigned absolute polarisation angles and no acceptable neighbouring pixels remain in the border list.

At this stage, Pacman solves the global $`n\pi `$-ambiguity $`n(k)`$ by applying a standard least squares fit to a set of best constrained pixels within the patch (i.e. the surrounding area of the reference pixel). Pacman solves for each of these best constrained pixels the global $`n\pi `$-ambiguity independently by either using the method of Haves (1975); Vallée & Kronberg (1975) or of Ruzmaikin & Sokoloff (1979) minimising the $`\chi ^2`$ for all possible combinations of $`n\pi `$. It searches then for the solution $`n(k)`$, which has been derived for the majority of pixels. This democratically determined solution $`n(k)`$ is used to perform a least squares fit of the other pixels defined within this patch. Note that in an extreme case, only the reference pixel could be used to solve the global $`n\pi `$-ambiguity.

If there are still good quality pixels left in the map which were not connected to the first patch of good quality data, the algorithm begins a new patch of absolute polarisation angles, starting again with the best remaining pixel with smallest $`q_{ij}`$. A new patch of absolute polarisation angles is constructed. As a result, the source will be divided in a set of spatially independent polarisation patches, and the global $`n\pi `$-ambiguity is solved once for each patch separately. A $`q_{ij}^{\mathrm{thresh}}`$ can be introduced to prevent Pacman from starting new patches if the pixels remaining have quality values above that threshold value whose value is chosen in the beginning.

### 2.3 Improving Quality

The solution of the local $`n\pi `$-ambiguity for pixel $`(ij)`$ becomes more reliable as more direct neighbours of these pixels have their absolute polarisation angle $`\overline{\phi }_{i^{}j^{}}`$ defined. This is taken into account by modifying the quality $`q_{ij}`$ of the pixels added to the border list to

$$1/q_{ij}=1/\sigma _{ij}^{RM}+\frac{\alpha }{n_{D_{ij}}^\beta }\left(\underset{i^{}j^{}D_{ij}}{}1/\sigma _{i^{}j^{}}^{RM}\right),$$
(7)
where $`n_{D_{ij}}`$ is the number of already defined neighbours $`(i^{}j^{})`$, and $`\alpha `$ and $`\beta `$ are free parameters. This ensures that pixels having more already defined neighbours are considered before others although the original data point might have a smaller signal quality, i.e a higher $`q_{ij}`$. Values between 0 and 1 for the free parameter $`\beta `$ are reasonable and yield good results. We used $`\alpha =\beta =1`$ and observed that Pacman goes through the acceptable data points of the polarisation patches in a uniform manner.

### 2.4 Restricting Gradients

The algorithm might be faced with a situation where a high quality pixel can be influenced by a very poor quality pixel, for example when the polarised radio source consists of two radio lobes each having a good signal-to-noise ratio which are connected by a bridge containing only low quality pixels. The Pacman algorithm would start by defining absolute polarisation angles from one of the two lobes, eventually reaching the bridge of low signal-to-noise and entering the second lobe from there. In such cases, it might happen that within the area of low quality data pixels a distinct determination of absolute polarisation angles by solving the local $`n\pi `$-ambiguity using Eq. (5) is no longer possible. The entire second lobe would then suffer, and wrong solutions would be introduced.

In order to avoid such situations, the algorithm is restricted to accept in the border list only neighbours which have a lower quality than the one under consideration. This forces the algorithm to go always from high to low quality pixels, leading to an artificial splitting of connected regions in the map into different patches.

However, such a strict rule would lead to heavy fragmentation and is not preferable. Therefore, we introduced a parameter $`g`$ to relax this strict limitation such that a new pixel is only accepted in the border list when the relation

$$q_{i^{}j^{}}>gq_{ij}$$
(8)
is fulfilled. This relation is always accounted for when adding a new pixel to the border list. We found that $`g`$ between 1.1 and 1.5 is a good choice.

### 2.5 Topological Defects

The algorithm can also be faced with situations which we call topological defects. These defects can be understood by supposing a ring-like polarised structure for which all polarisation angles point towards the centre of the ring. Starting at any point of the ring to define absolute polarisation angles and following the ring structure, a jump will appear on the border of the first to the last defined absolute polarisation angles after having performed a full circle. Scheuer et al. (1977) describe this problem which they encountered by their analysis of polarisation data.

In observational data, topological defects are often more complex structured. The artificial jump introduced will cause difficulties in the solution of the local $`n\pi `$-ambiguity because this procedure relies on already defined absolute polarisation angles of neighbouring pixels. When confronted with this situation, our algorithm divides the list of direct neighbours into sublists possessing similar polarisation angles. The sublist containing the most pixels is then used to assign the absolute polarisation angle by solving the local $`n\pi `$-ambiguity for the pixel under consideration.

The locations of the polarisation angle steps of topological defects are somewhat artificial since they depend on the actual path of the algorithm through the data. Since the algorithm processes all frequencies simultaneously, the steps of topological defects are at the same positions within all frequency maps, and therefore do not cause any further problems. However, if Pacman encounters a jump in the polarisation angle to all possible neighbour sublists in any of the frequencies under consideration, Pacman will reject this pixel and this pixel will not be considered for this patch but queued back for consideration for the next patches to be constructed. In our experience, these topological defects are rare events, but it is necessary to take them into account for the solution of the local $`n\pi `$-ambiguity.

### 2.6 Spurious Points

The intrinsic polarisation angle distribution might show strong gradients extending over a few pixels. This could lead to situations, where our algorithm can not solve the local $`n\pi `$-ambiguity at all frequencies simultaneously. Therefore, Pacman refuses to assign an absolute value $`\overline{\phi }_{ij}`$ to pixels when Eq. (4) yields a value $`\sigma _{ij}^\mathrm{\Delta }`$ above a certain threshold $`\sigma _{\mathrm{max}}^\mathrm{\Delta }`$, which can be set at the beginning of the calculation.

The pixels in the regions where this might occur most often have a low signal-to-noise ratio. In our experience, such situations always occur at strongly depolarised areas, leading to blanked regions in the RM distribution.

### 2.7 Multi-Frequency Fits

The aim of any RM derivation algorithm should be to calculate RMs for an area as large as possible using as much information as is available. On the other hand, for radio observations the total radio intensity decreases with increasing frequency. This can lead to the problem that the area of acceptable polarisation data at a high frequency is much smaller than at a lower frequency. This is especially true for (diffuse) extended radio sources. Furthermore, the limit of allowed standard errors $`\sigma _k^{\mathrm{max}}`$ of the polarisation angle might be exceeded for only one frequency leaving the values for the other frequencies still in the acceptable range.

Pacman can take this into account and performs the RM fit by omitting the polarisation angles at frequencies which do not meet the quality requirements. In order to do that, Pacman uses the standard errors $`\sigma _{k_{ij}}`$ of polarisation angles to define independently for each frequency, areas in which the RM map is permitted to be produced. An additional parameter $`k_{\mathrm{min}}f`$ is introduced which describes the minimum number of frequencies allowing the freedom to use any combination of the minimum of frequencies.

In cases of $`k_{\mathrm{min}}<f`$, Pacman will start to determine the solution to global and local $`n\pi `$-ambiguity only for the pixels fulfilling the quality requirements at all $`f`$ frequencies. After finishing that, the algorithm proceeds to include pixels satisfying the quality criteria in less than $`f`$ frequencies. For these pixels, the same patch building procedure as described in Sect. 2.2 applies with some modification in order to prevent the final maps from heavy fragmentation. The best quality pixel among the remaining pixel is picked but before starting a new patch, Pacman tests if the pixel under consideration adjoins a patch which has already been processed before. If the pixel adjoins such a patch, Pacman tries to solve the local $`n\pi `$-ambiguity following Eq.(6) and applies the patch solution of the global $`n\pi `$-ambiguity to the pixel. If the pixel is neither adjoined to an already processed patch nor the local $`n\pi `$-ambiguity solvable (i.e. $`\sigma _{ij}^\mathrm{\Delta }`$ exceeds $`\sigma _{\mathrm{max}}^\mathrm{\Delta }`$) then a new patch is initialised by this pixel.

Another possibility is to force Pacman to use certain frequencies in all circumstances. Thus, if the quality requirements are not fulfilled for these particular frequencies, there will be no RM value determined for the pixel under consideration. This has the advantage that one can use relatively close frequencies as a basis and then include other frequencies at points when a reliable polarisation signal is detected. The advantage of this technique is discussed in Paper II.

### 2.8 Additional Information

Apart from the resulting RM and $`\phi ^0`$ maps, Pacman provides sets of additional information about the data in order to estimate the reliability of the results. A patch map which contains all patches used is one example. Such a map is very useful, especially if one requires a minimal number of pixels in a patch in order to accept any calculated RM values from a particular patch. A map which includes rejected and thus flagged pixels, can also be obtained. Probably more important are the final $`\chi ^2`$\- and $`\sigma ^{\mathrm{RM}}`$\- maps which are also provided by *Pacman*. This information allows one to understand the reliability of the RM maps obtained and can be used for further evaluation and analysis of the RM maps.

## 3 testing the Pacman algorithm

In order to demonstrate the ability of Pacman to solve the $`n\pi `$-ambiguity, we use artificially generated maps. To generate them, we start from a polarisation data set from Abell 2255 (Govoni et al., 2002), which was kindly provided by Federica Govoni. We calculated the $`\mathrm{RM}`$ and $`\phi _0`$ maps from this data set and assume that these two maps consist of exact values. We generated then polarisation maps at four different frequencies which would exactly result from the $`\mathrm{RM}`$ and $`\phi _0`$ maps. As error maps, we generated Gaussian deviates for each frequency which were multiplied by the original error maps provided by Federica Govoni. The so generated error maps were added to the generated frequency maps in order to provide realistic mock observations.

Pacman and the standard fit algorithm were used to calculate the corresponding ’observed’ $`\mathrm{RM}`$ and $`\phi _0`$ maps which then were compared pixel by pixel to the initial exact maps. The result of this comparison is shown in Fig. 2. On the right panel, the comparison between the $`\mathrm{RM}_{\mathrm{pacman}}`$ values of the Pacman map and the values $`\mathrm{RM}_{\mathrm{sim}}`$ of the initial map is shown. The scatter in the data is due to the noise which was added to the frequency maps. On the left panel of Fig. 2, the pixel by pixel comparison between the values of the standard fit $`\mathrm{RM}_{\mathrm{stdfit}}`$ map and the initial $`\mathrm{RM}_{\mathrm{sim}}`$ map is shown. Again the scatter is due to the added noise. However, one can clearly see the points at $`\pm `$ 1000 rad m-2 which deviate from the initial data and are due to the wrongly solved $`n\pi `$-ambiguity in the case of the standard fit. Thus, this test demonstrates that Pacman yields reliable results to the $`n\pi `$-ambiguity.

## 4 Conclusions

We have presented a new algorithm for the calculation of Faraday rotation maps from multi-frequency polarisation data sets. Unlike other methods, our algorithm uses global information and connects information about individual neighbouring pixels with one another. It assumes that if regions exhibit small polarisation angle gradients between neighbouring pixels in all observed frequencies simultaneously, then these pixels can be considered as connected, and information about one pixel can be used for neighbouring ones, and especially, the solution to the $`n\pi `$-ambiguity should be the same. We like to stress, that this is a very weak assumption, and – like all other criteria used within the Pacman algorithm – only depends on the polarisation data at hand and the signal-to-noise ratio of the observations.

Our Pacman algorithm is especially useful for the calculation of RM and $`\phi ^0`$ maps of extended radio sources. Global algorithms as implemented in Pacman are preferable for the calculation of RM maps and needed if reliable RM values are desired from low signal-to-noise regions of the source.

Our Pacman algorithm reduces $`n\pi `$-artefacts in noisy regions and makes the unambiguous determination of RM and $`\phi ^0`$ in these regions possible. Pacman allows to use all available information obtained by the observation. Any statistical analysis of the RM maps will profit from these improvements.

A more detailed description of the application of Pacman to observations and various statistical tests performed on resulting RM maps in comparison to the standard fit algorithms are given in Paper II.

## acknowledgements

We want to thank Federica Govoni for allowing us to use their data on Abell 2255 to test the new algorithm. We like to thank Tracy Clarke, James Anderson, Phil Kronberg and an anonymous referee for useful comments on the manuscript. K. D. acknowledges support by a Marie Curie Fellowship of the European Community program ”Human Potential” under contract number MCFI-2001-01227.


text:
# A Comparison of different measures for dynamical event-mean transverse momentum fluctuation

thanks: This work is supported in part by the NSFC under project 90103019.

## Abstract

Various measures for the dynamical event-mean transverse momentum fluctuation are compared with the real dynamical fluctuation using a Monte Carlo model. The variance calculated from the $`G`$-moments can reproduce the dynamical variance well, while those obtained through subtraction procedures are approximate measures for not very low multiplicity. The $`\mathrm{\Phi }_{p_\mathrm{t}}`$ can also serve as an approximate measure after divided by the square root of mean multiplicity.

Multiparticle production; Dynamical fluctuation; Event-by-event analysis

The successful operation of the relativistic heavy ion collider rhic at bnl has stimulated an intensive study on the theoretically predicted quark gluon plasma (qgp) phase transition. It is expected that the appearance of phase transition will result in an anomalous behavior in the event-by-event fluctuations of transverse momentum Mrowski Shuryak Pisarski .

Since the multiplicity in a single event is finite, there are inevitably statistical fluctuations (sf) in the event-by-event analysis of transverse momentum. Various methods PHENIX Voloshin CERES  have been proposed to get rid of the statistical fluctuations and measure the dynamical ones. Most of them are based on a subtraction procedure, i.e. first to estimate the variance of statistical fluctuations and then subtract it from that of the experimental data. The basic problem in this kind of methods is how to reliably estimate the statistical variance.

Another approach proposed recently FuLiu  is to eliminate the sf directly from the experimental data. This method can reproduce the dynamical variance exactly provided the sf are Poissonian, i.e. the particles are emitted independently.

The aim of this paper is to compare the different measures proposed in the market and, in particular, to examine the reliability of the estimation of statistical variance. The so called $`\mathrm{\Phi }_{p_\mathrm{t}}`$ measure marek  will also be discussed and a new measure $`\mathrm{\Phi }_{p_\mathrm{t}}^{(\mathrm{r})}`$ will be proposed, which can largely reduce the strong multiplicity dependence of $`\mathrm{\Phi }_{p_\mathrm{t}}`$ itself.

Let $`p(p_\mathrm{t})`$ be the dynamical transverse-momentum distribution in single events and

$$p_m=_{\delta _m}p(p_\mathrm{t})dp_\mathrm{t},(m=1,2,\mathrm{},M)$$
(1)
be the corresponding “coarse-grained” distribution BP . The realization of $`p(p_\mathrm{t})`$ in experiment is the distribution of the total number $`n`$ of particles in the $`p_\mathrm{t}`$ region $`\mathrm{\Delta }`$, and the expression

$$q_m=n_m/n(m=1,2,\mathrm{},M)$$
(2)
is an evaluation of $`p_m`$, where $`n_m`$ is the number of particles falling into the $`m`$th bin. Thus the dynamical and experimentally-measured event-mean transverse momentum are, respectively,

$$\overline{p_\mathrm{t}}_{\mathrm{dyn}}=_\mathrm{\Delta }p_\mathrm{t}p(p_\mathrm{t})dp_\mathrm{t}=\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_mp_m,$$
(3)
$$\overline{p_\mathrm{t}}_{\mathrm{exp}}=\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_mq_m=\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_m\frac{n_m}{n},$$
(4)
where $`(p_\mathrm{t})_m`$ is the $`p_\mathrm{t}`$ value in the $`m`$th bin. The event-space moments of $`\overline{p_\mathrm{t}}_{\mathrm{dyn}}`$ and $`\overline{p_\mathrm{t}}_{\mathrm{exp}}`$ are

$$C_p^{\mathrm{dyn}}(\overline{p_\mathrm{t}})=\overline{p_\mathrm{t}}_{\mathrm{dyn}}^p=\left(\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_mp_m\right)^p,$$
(5)
$$C_p^{\mathrm{exp}}(\overline{p_\mathrm{t}})=(\overline{p_\mathrm{t}}_{\mathrm{exp}})^p=\left(\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_m\frac{n_m}{n}\right)^p,$$
(6)
respectively. These two are evidently unequal due to the existence of sf. In particular, the experimental moments $`C_p^{\mathrm{exp}}(\overline{p_\mathrm{t}})`$ depend crucially on the particle multiplicity NA49 , which is mainly due to sf. Therefore, subtraction procedures are commonly used to extract the dynamical fluctuation of $`\overline{p_\mathrm{t}}`$.

The basic idea of the subtraction procedure is that the variances $`\sigma ^2(\overline{p_\mathrm{t}})=C_2(\overline{p_\mathrm{t}})C_1(\overline{p_\mathrm{t}})^2`$ of dynamical and statistical fluctuations are additive

$$\sigma _{\overline{p_\mathrm{t}}\mathrm{data}}^2=\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^2+\sigma _{\overline{p_\mathrm{t}}\mathrm{stat}}^2,$$
(7)
where $`\sigma _{\overline{p_\mathrm{t}}\mathrm{data}}^2`$, $`\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^2`$ and $`\sigma _{\overline{p_\mathrm{t}}\mathrm{stat}}^2`$ are the experimental, dynamical and statistical variances, respectively. This additivity is true if the df and sf are independent, which is a reasonable assumption. Therefore, the main problem in the subtraction procedure for extracting the dynamical variance is to estimate the statistical ones. Various methods have been proposed for this purpose.

For example, in Ref. PHENIX  the results from mixed events are considered as the baseline for the random distribution and the difference in the fluctuation from a random distribution defined as

$$d=\omega _{\mathrm{data}}\omega _{\mathrm{baseline}}$$
(8)
is taken as a measure of the dynamical fluctuation. In Eq.(8)

$$\omega =\frac{\sqrt{\overline{p_\mathrm{t}}^2\overline{p_\mathrm{t}}^2}}{\overline{p_\mathrm{t}}}=\frac{\sqrt{\sigma _{\overline{p_\mathrm{t}}}^2}}{\overline{p_\mathrm{t}}},$$
(9)
$`\overline{p_\mathrm{t}}`$ is the mean transverse momentum in a single event and $`\mathrm{}`$ denotes the average over event sample.

Alternatively, in Ref. Voloshin  the statistical variance of event mean $`p_\mathrm{t}`$, under the assumption of independent particle production, is estimated as

$$\sigma _{\overline{p_\mathrm{t}}\mathrm{stat}}^2=\frac{\sigma _{p_\mathrm{t}\mathrm{incl}}^2}{n},$$
(10)
and, therefore, the dynamical variance is equal to

$$\sigma _{\overline{p_\mathrm{t}}_{\mathrm{VKR}}}^2=\sigma _{\overline{p_\mathrm{t}}\mathrm{data}}^2\sigma _{\overline{p_\mathrm{t}}\mathrm{stat}}^2=\sigma _{\overline{p_\mathrm{t}}\mathrm{data}}^2\frac{\sigma _{p_\mathrm{t}\mathrm{incl}}^2}{n}.$$
(11)
A somewhat different expression is used in Ref. CERES , where the experimental variance is weighted by event multiplicity $`N_j`$:

$$\sigma _{\overline{p_\mathrm{t}}_{\mathrm{CERES}}}^2=\frac{_{j=1}^𝒩N_j(\overline{p_\mathrm{t}}^j\overline{p_\mathrm{t}})^2}{_{j=1}^𝒩N_j}\frac{\sigma _{p_\mathrm{t}\mathrm{incl}}^2}{n},$$
(12)
where $`𝒩`$ is the total number of events.

A widely used measure for the non-statistical mean $`p_\mathrm{t}`$ fluctuation is the $`\mathrm{\Phi }_{p_\mathrm{t}}`$ proposed in Ref. marek

$$\mathrm{\Phi }_{p_\mathrm{t}}\sqrt{Z^2/n}\sqrt{z^2},$$
(13)
where $`z`$ and $`Z`$ are defined as $`zp_\mathrm{t}p_\mathrm{t}`$ for each particle and $`Z_{i=1}^nz_i=n(\overline{p_\mathrm{t}}p_\mathrm{t})`$ for each event, respectively note2 . The second term in the r.h.s. of Eq. (13) is the square root of the inclusive variance $`\sigma _{p_\mathrm{t}\mathrm{incl}}^2=(p_\mathrm{t}p_\mathrm{t})^2`$. Assuming that the multiplicity fluctuation is uncorrelated with the $`p_\mathrm{t}`$ fluctuation, we get from Eq.(13)

$$\mathrm{\Phi }_{p_\mathrm{t}}=\sqrt{n^2/n}\sigma _{\overline{p_\mathrm{t}}\mathrm{data}}\sigma _{p_\mathrm{t}\mathrm{incl}}.$$
(14)
This equation is evidently similar to Eq.(11).

All of the above measures have the same structure, being based on a subtraction procedure, i.e. to subtract the variance of $`\overline{p_\mathrm{t}}`$ or a quantity related to it, that will be expected from a pure statistical system, from the same quantity obtained in experiment. These measures will, of course, vanish for a pure statistical system, and a non-vanishing value of them will indicate the existence of dynamical effect. Therefore, the measures based on the subtraction procedure, as those listed above, will at least qualitatively measures the effect of df.

In order to have an idea on how the various subtraction methods work quantitatively, let us carry on a Monte Carlo simulation using a toy model. In this model the distribution of $`p_\mathrm{t}`$ is taken as

$$p(p_\mathrm{t})=\frac{4}{a^2}p_\mathrm{t}\mathrm{e}^{2p_\mathrm{t}/a},$$
(15)
with a Gaussian distributed parameter $`a`$ ($`\sigma ^2(a)=0.24`$). In total $`20\times 500,000`$ events are generated for 20 average multiplicities. The resulting dynamical and experimental variances $`\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^2`$ and $`\sigma _{\overline{p_\mathrm{t}}\mathrm{exp}}^2`$ are plotted in Fig.1 as dashed line and open triangles, respectively. For comparison, the $`\sigma _{\overline{p_\mathrm{t}}\mathrm{VKR}}^2`$ and $`\sigma _{\overline{p_\mathrm{t}}\mathrm{CERES}}^2`$ are also plotted in the same figure as upward and downward solid triangles, respectively.

The differences between the experimental and dynamical variances $`\sigma _{\overline{p_\mathrm{t}}\mathrm{exp}}^2\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^2`$ are due to sf. It can be seen from the figure that the $`\sigma _{\overline{p_\mathrm{t}}\mathrm{VKR}}^2`$ indeed reduces the sf, especially when the average multiplicity is not very low; while the $`\sigma _{\overline{p_\mathrm{t}}\mathrm{CERES}}^2`$ overshoots the influence of sf.

On the other hand, the open circles shown in the figure reproduce the dynamical variances very well. They are obtained through a totally different approach FuLiu , which is based on the elimination of the sf directly from the experimental data.

For the first order moment the elimination of sf is straightforward,

$$C_1^{(\mathrm{dyn})}(\overline{p_\mathrm{t}})=G_1(\overline{p_\mathrm{t}})=\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_m\frac{n_m}{n},$$
(16)
provided the sf is Poissonian. The key point in eliminating sf in the second order moment is to expand the 2nd power in the definition Eq.(5) of $`C_2^{(\mathrm{dyn})}(\overline{p_\mathrm{t}})`$,

$$C_2^{(\mathrm{dyn})}(\overline{p_\mathrm{t}})=\left(\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_mp_m\right)^2=\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_m^2p_m^2+\underset{mm^{}}{\overset{M}{}}(p_\mathrm{t})_m(p_\mathrm{t})_m^{}p_mp_m^{}.$$
(17)
Using the formulae

$$\underset{m=1}{\overset{M}{}}f_mp_m^p=\underset{m=1}{\overset{M}{}}f_m\frac{n_m(n_m1)\mathrm{}(n_mp+1)}{n^p},$$
(18)
which holds for Poissonian sf, we then get $`C_2^{(\mathrm{dyn})}(\overline{p_\mathrm{t}})=G_2(\overline{p_\mathrm{t}})`$, where

$$G_2(\overline{p_\mathrm{t}})=\underset{m=1}{\overset{M}{}}(p_\mathrm{t})_m^2\frac{n_m(n_m1)}{n^2}+\underset{mm^{}}{\overset{M}{}}(p_\mathrm{t})_m(p_\mathrm{t})_m^{}\frac{n_mn_m^{}}{n^2}.$$
(19)
This guarantees the equality of $`G`$-variance $`\sigma _{\overline{p_\mathrm{t}}G}^2=G_2(\overline{p_\mathrm{t}})G_1(\overline{p_\mathrm{t}})^2`$ to dynamical ones $`\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^2=C_2^{(\mathrm{dyn})}(\overline{p_\mathrm{t}})C_1^{(\mathrm{dyn})}(\overline{p_\mathrm{t}})^2.`$

On the contrary, the sf cannot be exactly gotten rid of using the subtraction procedure. The main problem lies in the estimation of statistical variance by Eq.(10). This equality holds only for a pure statistical system without any dynamical fluctuation. When dynamical fluctuations are present the $`\sigma _{\overline{p_\mathrm{t}}\mathrm{stat}}^2`$ and $`\sigma _{p_\mathrm{t}\mathrm{incl}}^2/n`$ are, in general, unequal. For an exhibition, the $`\sigma _{\overline{p_\mathrm{t}}\mathrm{stat}}^2=\sigma _{\overline{p_\mathrm{t}}\mathrm{exp}}^2\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^2`$, cf. Fig.1, and $`\sigma _{p_\mathrm{t}\mathrm{incl}}^2/n`$ are plotted in Fig.2 as solid circles and solid triangles, respectively. They are approximately equal to each other only when the average multiplicity is not very low, e.g. when $`n>12`$.

In order to further clarify this point, we have made a pure statistical model with $`p_m=1/M`$, cf. Eq.1. In this case $`\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^20`$. The corresponding $`\sigma _{\overline{p_\mathrm{t}}\mathrm{stat}}^2=\sigma _{\overline{p_\mathrm{t}}\mathrm{exp}}^2`$ and $`\sigma _{p_\mathrm{t}\mathrm{incl}}^2/n`$ are plotted in Fig.2 as open circles and open triangles, respectively. The differences between them are small up to $`n=5`$. This shows that the discrepancy between $`\sigma _{\overline{p_\mathrm{t}}\mathrm{VKR}}^2`$ ($`\sigma _{\overline{p_\mathrm{t}}\mathrm{CERES}}^2`$) and $`\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^2`$ is mainly due to the inappropriate estimation at low multiplicity of $`\sigma _{\overline{p_\mathrm{t}}\mathrm{stat}}^2`$ by $`\sigma _{p_\mathrm{t}\mathrm{incl}}^2/n`$ for a system containing df.

Let us now turn to the discussion on the $`\mathrm{\Phi }_{p_\mathrm{t}}`$ measure marek NA49 . The Monte Carlo results of $`\mathrm{\Phi }_{p_\mathrm{t}}`$ are plotted in Fig.3 as solid circles. It is strongly multiplicity dependent, and is, therefore, not a good measure for df. On accounting of Eq.(14) we divide $`\mathrm{\Phi }_{p_\mathrm{t}}`$ by the square root of average multiplicity and define

$$\mathrm{\Phi }_{p_\mathrm{t}}^{(\mathrm{r})}\frac{\mathrm{\Phi }_{p_\mathrm{t}}}{\sqrt{n}}\sigma _{\overline{p_\mathrm{t}}\mathrm{data}}\frac{\sigma _{p_\mathrm{t}\mathrm{incl}}}{\sqrt{n}},$$
(20)
cf. Eq.(11).

The Monte Carlo results of $`\mathrm{\Phi }_{p_\mathrm{t}}^{(\mathrm{r})}`$ are shown in Fig.3 as upward solid triangles together with $`\sigma _{\overline{p_\mathrm{t}}_{\mathrm{CERES}}}=\sqrt{\sigma _{\overline{p_\mathrm{t}}_{\mathrm{CERES}}}^2}`$ (downward triangles) and $`\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}`$ (dashed line). It can be seen from the figure that $`\mathrm{\Phi }_{p_\mathrm{t}}^{(\mathrm{r})}`$ has a much weaker multiplicity dependence and is closer to $`\sigma _{\overline{p_\mathrm{t}}_{\mathrm{dyn}}}`$ than $`\mathrm{\Phi }_{p_\mathrm{t}}`$. Therefore, $`\mathrm{\Phi }_{p_\mathrm{t}}^{(\mathrm{r})}`$ is a better measure of df than $`\mathrm{\Phi }_{p_\mathrm{t}}`$ itself.

In this paper, various measures for the dynamical transverse momentum fluctuation in event-by-event analysis are compared with the dynamical fluctuation using a Monte Carlo model.

It turns out that the $`\sigma _{\overline{p_\mathrm{t}}G}^2`$ calculated from the G-moments coincides with the dynamical variance $`\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}^2`$, showing that the G-moment method is effective in eliminating the statistical fluctuations coming from the finite number of particle in a single event and thus provides a good measure for the dynamical fluctuations.

The $`\sigma _{\overline{p_\mathrm{t}}\mathrm{VKR}}^2`$ proposed by S. A. Voloshin, V. Koch and H. G. Ritter Voloshin  and its revised version $`\sigma _{\overline{p_\mathrm{t}}\mathrm{CERES}}^2`$ CERES  are found to be good approximate measures for the dynamical fluctuation when the multiplicity is not very low. Therefore, their application to relativistic heavy ion experiments is justified.

The $`\mathrm{\Phi }_{p_\mathrm{t}}`$ divided by the square root of average multiplicity ($`\mathrm{\Phi }_{p_\mathrm{t}}^\mathrm{r}=\mathrm{\Phi }_{p_\mathrm{t}}/\sqrt{n}`$) is another approximate measure for the dynamical fluctuation, which depends on multiplicity weaker than $`\mathrm{\Phi }_{p_\mathrm{t}}`$ and is closer to the square root of dynamical variance $`\sigma _{\overline{p_\mathrm{t}}\mathrm{dyn}}`$.

Acknowledgement The authors thank Wu Yuanfang and Liu Feng for helpful discussions and comments.


text:
# Shear band formation in granular media as a variational problem

(; version 1.0)

## Abstract

Strain in sheared dense granular material is often localized in a narrow region called shear band. Recent experiments in a modified Couette cell Fenistein2003a ; Fenistein2003b  provided localized shear flow in the bulk away from the confining walls. The non-trivial shape of the shear band was measured as the function of the cell geometry. First we present a geometric argument for narrow shear bands which connects the function of their surface position with the shape in the bulk. Assuming a simple dissipation mechanism we show that the principle of minimum dissipation of energy provides a good description of the shape function. Furthermore, we discuss the possibility and behavior of shear bands which are detached from the free surface and are entirely covered in the bulk.

granular flow, shear band, variational principle, least dissipation, optimization

Granular media constitute an interesting field of research from the point of view of both basic science and application. The intrinsically nonlinear and dissipative nature of the interaction between the particles leads to a great deal of interesting phenomena including force networks, different kinds of instabilities, clustering, complex flow properties DryGranu98 ; Jaeger96 . One of the most apparent instabilities occurring in granular media is the formation of shear bands: At slow shear rate strain is not distributed throughout the sample but appears in a localized fashion along a rather narrow interface between two essentially unstrained parts. Shear-banding was the subject of various experimental and theoretical studies in the last few years DryGranu98 ; Mueth00 ; Hartley03 ; Thompson91 ; Torok00 ; Schwedes03 ; Veje99  and still presents significant difficulties for theoretical descriptions.

Recently universal geometrical properties of shear bands were discovered in a Couette geometry modified such that shear localization near boundaries was avoided Fenistein2003a ; Fenistein2003b . The experimental setup was a cylindrical container filled with grains up to height $`H`$. The bottom was split into an outer ring rotating with the container wall, and a stationary disk of radius $`R_s`$ in the center (Fig. 1)111The experiments Fenistein2003a ; Fenistein2003b  were performed both with “disk” geometry and with a split-bottomed Couette geometry, where an additional stationary inner cylinder is present. In this paper only the “disk” geometry is discussed.. Thus the outer and the inner part of the material were rotated relative to each other which created a shear band with cylindrical symmetry: It started at the perimeter of the stationary bottom disk and extended through the bulk up to the free surface. On the surface the angular velocity of the granular material as the function of the radius was measured. It follows to high accuracy an error function characterized by two parameters: the width $`W`$ and the center position $`R_c`$ of the shear zone. The width grows with increasing height $`H`$ while the radius $`R_c`$ gets smaller. Interestingly, the surface position $`R_c`$ proved to be very robust: It depends only on two length parameters $`H`$ and $`R_s`$, but not on the particle properties nor on the shear rate. This is contrary to the width of the band, which is affected by the size and shape of the grains but is insensitive to changes of the slip radius $`R_s`$. The fact that $`R_c`$ and $`W`$ depend on different control parameters suggests that the two quantities can be studied separately.

In this Letter we address the problem of the position of the shear band. By suitable choice e.g. of the particle size the width can be made arbitrarily small compared to $`H`$ and $`R_s`$. This justifies to model the shear band as an infinitely thin layer that represents the boundary between two blocks of material within which no flow occurs. We derive the shape of the band from a *variational principle*. Note that optimization principle was applied already to shear band formation in a different context Torok00 .

Surface - bulk relation. In the experiments the position $`R_c(R_s,H)`$ on the free surface was found Fenistein2003b  to scale like

$$R_c(R_s,H)=R_s(1(H/R_s)^\alpha )$$
(1)
with $`\alpha 2.5`$ for the experimentally accessible values of $`H/R_s`$. It is much more difficult to measure the position of the shear band in the bulk, $`r(h)`$, for fixed $`R_s`$ and $`H`$. Nonetheless, the experimental data clearly show that the bulk profile follows another form than $`R_c(R_s,H)`$, and the bulk radius at height $`h`$ depends also on the filling height $`H`$ (Fig. 2).

We show that *the bulk profile is determined once the surface positions $`R_c(R_s,H)`$ are given*. Let us take a system with total height $`H`$ and find the position of the shear band $`r`$ at height $`h<H`$. The subsystem above $`h`$ can be regarded as a smaller system with height $`(Hh)`$ and with slip radius $`r`$ at the bottom. Pressure and boundary conditions are the same, and the difference in the width is neglected in the narrow band approximation. We conclude that

$$R_c(R_s,H)=R_c(r,Hh).$$
(2)
Knowing the function $`R_c(R_s,H)`$ thus allows to calculate the shape $`r(h)`$ of the shear band throughout the whole system. Note, if the simple size scaling holds for the surface radius that $`R_c/R_s`$ depends only on $`H/R_s`$ (as it was found in Fenistein2003b ) then it follows immediately also for the bulk function due to Eq. (2):

$$\frac{r}{R_s}=f(\frac{H}{R_s},\frac{h}{R_s}).$$
(3)
Based on Eq. (2) an explicit functional form of the bulk profile can be obtained using the experimental fit function of $`R_c`$:

$$h=Hr\left[1\frac{R_s}{r}(1(H/R_s)^\alpha )\right]^{1/\alpha }.$$
(4)
The resulting curves for some filling heights are plotted in Fig. 2 using $`\alpha =2.5`$ and the comparison with the experimental data shows very good agreement.

Variational principle. In order to describe the form of the shear band our idea is to apply the *principle of least dissipation* Onsager1931  (which is a common treatment of time-independent irreversible phenomena statphys ). Therefore we require a steady state flow that matches the outer constraints but provides the minimum rate of energy dissipation.

Applying this to the cylindrical geometry within the narrow band approximation the question of the shape is traced back to a variational problem among the functions $`r(h)`$ (where $`H`$ and $`R_s`$ are kept fixed) with the condition $`r(0)=R_s`$ while the other boundary at $`H`$ is free. The dissipation rate is given by the sliding velocity between the two sides $`r(h)\omega `$ times the shear stress $`\sigma _{tn}`$ integrated over the whole shear band. Up to a constant factor, the expression to be minimized is

$$_0^Hr^2\sqrt{1+\left(dr/dh\right)^2}\sigma _{tn}𝑑h=\text{min.}$$
(5)
This quantity represents not only the *dissipation rate* but also the *mechanical torque* that the rotating and stationary part of the system exert on each other. Therefore the *least dissipation* for this specific geometry is equivalent to the *minimal torque* which gives further justification of this approach: it is plausible that the yielding surface is established where the resistance against the outer constraint is the smallest, i.e. where the material is the weakest Torok00 .

Sliding model. For the shear stress in Eq. 5 we will use a very simple sliding model. The shear stress in the yielding surface is taken similar to the Coulomb friction between two solid bodies: It acts against the sliding direction, its magnitude is proportional to the normal pressure pressing the two sides against each other, but it is independent of the sliding velocity. We assume hydrostatic pressure i.e. proportional to the depth 222Here we ignore anisotropy effects.. Thus Janssen effect is neglected, which is naturally justified if H is smaller than the container width. In our dynamical situation, however, we expect that the applicability of the hydrostatic pressure can be extended even for larger filling heights: The shear band (due to many collisions and slip events) acts as a source of small vibrations in the whole system and can cause slight creep at the particle-wall contacts inhibiting the particles to keep their original anchoring position. Finally they transmit their load to the next particle below rather than to the side wall and therefore the whole weight will be carried by the bottom.

This sliding model leads to the variational problem

$$_0^Hr^2\sqrt{1+\left(dr/dh\right)^2}(Hh)𝑑h=\text{min.}.$$
(6)
The solutions which minimize the integral have automatically the scaling property given by Eq. (3). The reason of this data collapse is that taking a $`\lambda `$ times larger system (i.e. taking $`\lambda R_s`$, $`\lambda H`$, $`\lambda r(h/\lambda )`$ instead of $`R_s`$, $`H`$, $`r(h)`$) changes the value of the integral only by a constant factor ($`\lambda ^4`$) thus it represents the same variational problem.

Numerical results. The function $`r(h)`$ is discretized and the minimization is performed numerically based on genetic optimization: $`r(h)`$ is varied randomly but only the changes that lower the value of the left hand side in Eq. (6) ,i.e., reduce the dissipation, are admitted. During the optimization the noise level is continuously decreased and the final state $`r(h)`$ is regarded as one local minimum of the variational problem. The landscape where the minimum has to be found is simple (see later).

Results shown in Fig. (3) reproduce nicely the qualitative behavior found in the experiment: the concave shear bands appear in the bulk and build up a convex confining shape of the surface positions as the filling height is varied. The shear radii at a fixed bulk height $`h`$ and also at the top get smaller with increasing filling height.

The realistic bulk profiles provided by the principle of minimum dissipation can be interpreted easily. The benefit gained by having such a curved cylindrical shape (slimmer at the top) in place of a regular cylinder is twofold: First, the surface of the shear band can be reduced by letting the regular cylinder deform and by pulling the shape towards the center a bit at a fixed bottom radius, similarly to a soap membrane spanned between a ring and a plate where the plane of the ring is parallel to the plate. Second, smaller radius results in smaller sliding velocity (or thinking of minimizing the torque it represents smaller lever). Therefore one can roughly think of the bulk profiles shown in Fig. (3) as equilibrium situations where the reduction of the sliding velocity is counterbalanced by the increase in the surface due to going beyond the minimum surface.

The quantitative agreement with the experimental fit function (Eq. 1) is also surprisingly good given the crude assumptions we made and the fact that our model contains no free fit parameters. The difference of the theoretical and the experimental values of $`R_c`$ is less than $`20`$ % of $`R_sR_c`$.

We can more easily analyze the limit $`H/R_s0`$ than in the experiment, where $`R_s`$-values are limited by the container size and a power law fit to the data at low filling heights is not very precise. From our variational principle we find an exponent $`\alpha =2`$ for small $`H/R_s`$. In the region $`0.4<H/R_s<0.7`$ there is no clean power law any more and the effective exponent increases. In our model this seems to be due to a phase transition at about $`H/R_s0.7`$ to a different shape of the shear band, which we discuss next.

What is the predicted behavior of the shear band for filling heights larger then those reported in Fenistein2003a ; Fenistein2003b  ? For large enough $`H/R_s`$ the class of the “open” solutions discussed so far is replaced by a new type of solutions: The shear band, instead of running up to the free surface and having a circle on the top as its upper edge, closes forming a cupola-like shape (still with bottom radius $`R_s`$). In that case the material covered by the “closed” shear band is at rest while the material around and above (including the whole free surface) is rotating. Several “open” and “closed” profiles can be seen in Fig. (4b) obtained for various values of $`H`$. Fig. (4a) shows the upper radius of the shear bands which is characteristic of the “open” solutions but becomes zero for the “closed” ones. For these cupola shapes 333The “closed” profiles appear in the integral of Eq. (6) in the way that the function $`r(h)`$ becomes zero for $`h_{\text{top}}hH`$. Therefore this region has no contribution to the dissipation and practically the integral is meant over the interval $`[0,h_{\text{top}}]`$. a more relevant parameter is their heights $`h_{\text{top}}H`$ in the center, plotted in Fig. (4c). For “open” profiles $`h_{\text{top}}`$ equals simply the system height.

Interesting is the behavior of the parameter $`h_{\text{top}}`$. It is a monotonically decreasing function of the filling height once it is detached from $`H`$, i.e. large filling heights press the shear band to the bottom. Solving the Euler-Lagrange equation for the variational problem neglecting terms of higher than first order in $`dh/dr`$ and taking the limit of $`Hh_{\text{top}}`$ one obtains

$$h(r)=h_{\text{top}}\frac{1}{6H}r^2,h_{\text{top}}=\frac{1}{6H}R_s^2.$$
(7)
Fig. (4d) shows the numerical solution of $`h_{\text{top}}`$ : it is in excellent agreement with the approximate analytical solution. In the limit $`H/R_s\mathrm{}`$ the material forms one solid block and the sliding occurs on the surface of the bottom disk.

The $`H`$ dependence of $`h_{\text{top}}`$ has the following reason: Larger $`h`$ corresponds to lower pressure and thus to smaller shear-resistance; therefore it is worth to raise the shear band into a cupola-shape even if its surface becomes larger than the disk at the bottom. Note, however, that stronger gravity would not affect the shape: It gives a larger pressure gradient but this constant factor has no impact on the variational problem. By contrast, increasing $`H`$ (or, alternatively, applying additional pressure at the surface) makes the relative pressure change smaller near the bottom which results in weaker uplifting “force”. This is why the cupola height $`h_{\text{top}}`$ approaches zero for large $`H`$ (respectively large pressure).

At intermediate $`H`$ (around $`0.7R_s`$) there exists a small region of filling heights where the variational problem (Eq. 6) exhibits two local minima. The behavior of the local minima corresponds to a first order phase transition, where the two phases are the two types of shear band profile (Fig 5). Following the “closed” solution as $`H`$ is decreased its height grows till the top of the cupola reaches the surface at system height $`H_1`$. At this point the “closed” solution becomes unstable, $`R_c`$ at the top runs suddenly to a non-zero value (Fig. 5a). Regarding the other direction (increasing $`H`$) the upper edge of the “open” shape is pulled towards the center then at height $`H_2`$ it shrinks to one point which is followed by a jump into a “closed” shape (discontinuity in $`h_{\text{top}}`$ see in Fig. 5b). The interval $`[H_1,H_2]`$ defines the region where two local minima exist, outside this region there is only one phase possible. Local minima are physically meaningful in the presence of a kinetic barrier. The latter can be provided by the difference between the static and dynamic friction coefficients, usually present in granular systems.

In this Letter we have presented a theoretical analysis of recent experiments  on shear band formation in granular media. We used the approximation of narrow shear bands. First we showed a geometric argument which related the position of the surface endpoints of the band to the bulk shape. We calculated the shape of shear bands from a variational principle and the results are in good agreement with the experiments. The theory provides with a number of predictions. i) There is a transition in the shape of the shear band as a function of the filling height $`H`$: For low values of $`H`$ the shear band is an open curved cylinder which ends on the surface while for large $`H`$ a cupola is formed. ii) This transition is of the first order, accompanied by a hysteresis. iii) The height of the cupola is proportional to $`R_s/H`$. These predictions should be experimentally accessible.

An interesting question we did not address in this Letter is the nucleation kinetics at the first order phase transition from open to closed shapes. It is conceivable that the finite width of the shear bands plays a role here. Further progress in this direction requires a continuum theory going beyond the narrow band approximation we used here.

We would like to thank to D. Fenistein and M. van Hecke for useful discussions. Partial support by grant OTKA T035028 and by the German-Hungarian Cooperation Fund is acknowledged.