text:
# Rashba precession in quantum wires with interaction

## Abstract

Rashba precession of spins moving along a one-dimensional quantum channel is calculated, accounting for Coulomb interactions. The Tomonaga–Luttinger model is formulated in the presence of spin-orbit scattering and solved by Bosonization. Increasing interaction strength at decreasing carrier density is found to enhance spin precession and the nominal Rashba parameter due to the decreasing spin velocity compared with the Fermi velocity. This result can elucidate the observed pronounced changes of the spin splitting on applied gate voltages which are estimated to influence the interface electric field in heterostructures only little.

PACS : 71.10.Pm, 71.70.Ej, 73.20.Dx

According to recent speculations spin could replace the electric charge to carry information in future electronic devices . The ‘spin transistor’ proposed by Datta and Das  might switch faster than traditional transistors since during operation it avoids redistributing charges. The idea is based on the Rashba effect  which causes spins to precess as they move along a heterostructure  so that the conductance depends on the final spin orientation relative to the magnetization of the ferromagnetic drain contact . The strength of the Rashba effect is proportional to the electric field acting perpendicular to the electron plane which, when varied by a gate, changes the final spin orientation and thus the transport properties of the device. Advantage is taken of long spin coherence times and lengths , found in semiconductors. In the attempt of experimental realization considerable progress has been achieved meanwhile to inject finite spin densities . Polarizations of 90% have been reported in GaAs .

Beating patterns varying with gate voltages have indeed been observed in Shubnikov–de Haas (SdH) measurements . Estimates show, however, that the the gate voltage adds an electric field contribution to a much stronger intrinsic field at the interface which does not suffice to explain the observed variations in spin splittings by a factor of two. One should notice that the gate voltage not only changes the strength of the electric field but at the same time it also alters the carrier density in the heterostructure as is directly monitored by the SdH oscillations. Without interactions this would not change the Rashba precession within the effective mass approximation  but, since particularly in semiconductors the strength of Coulomb interactions change with density through the $`r_\mathrm{s}`$–parameter, the question arises whether interactions influence the Rashba precession. This is demonstrated in the present Communications. In principle, the two effects of the gate, changing field strength and density, can be separated experimentally by varying the voltage at a front gate and a back gate independently .

Usually, Rashba spin precession is described as a band structure effect, resulting from spin splitting , as calculated in the original work by Rashba  for the homogeneous two-dimensional case. Such a single particle approach, however, cannot account for electron-electron interactions. Including many body effects requires going beyond band structure theory. Particularly convenient  and powerful  is the Tomonaga–Luttinger (TL) model to incorporate interactions exactly for the low energy and long wave length properties in systems of one spatial dimension. From the practical point of view, quantum wires seem most efficient to realize the spin transistor since guiding the electron waves avoids angular dispersion of momenta  and diminishes spin scattering events . Spin injection into carbon nanotubes, which are known to show pronounced TL behavior  has been demonstrated . Here we investigate a quantum wire with spin-orbit interaction in the one-subband limit and apply the TL–model. We expect our main result that interactions enhance Rashba precession to carry over to more channel situations. Indeed, a recent careful analysis within improved versions of the RPA approximation has revealed enhancement of the Rashba precession in the two dimensional electron system . One might comprehend these findings as a consequence of ‘repulsion’ between the two spin split bands by the interaction.

Spin splitting originates from spin–orbit coupling which in narrow gap heterostructures, such as those based on InAs, is dominated  by the Rashba term

$$H^{\mathrm{so}}=\alpha (\sigma _xp_z\sigma _zp_x)$$
(1)
where $`\sigma _{x,z}`$ are Pauli matrices and $`p`$ is the electron momentum in the $`x`$$`z`$–plane of the heterostructure. Here, we use coordinates indicated in the inset of Fig. 1. The Rashba parameter $`\alpha `$ is mostly determined by the intrinsic electric field perpendicular to the interface. Typical values of 1—3$`\times 10^{11}`$eVm  are reported for InAs .

Taking the confining potential in the lateral $`z`$–direction as parabolic, the eigenenergies of

$$H=\frac{(p_z+m\alpha \sigma _x)^2+p_x^2}{2m}+\frac{m}{2}\omega _0^2z^2\alpha \sigma _zp_x$$
(2)
determine the spin split single particle dispersion . Here, $`\omega _0`$ is the subband energy and effective mass ($`m`$) approximation is assumed. The eigenfunctions are plane waves of momentum $`k`$ in $`x`$–direction along the wire and, to good approximation, given as slightly modified oscillator functions $`\varphi _n`$ for the subband index $`n`$

$$\psi _{kns}(x,z)=\mathrm{e}^{\mathrm{i}kx}\varphi _n(z)(\mathrm{cos}(m\alpha z)|s+\mathrm{i}\mathrm{sin}(m\alpha z)|s).$$
(3)
On the wire axis the $`\psi _{kns}(x,z)`$ are spin polarized along the $`z`$–direction ($`s=\pm 1`$) while the spins acquire a non-zero out of plane $`y`$–component $`\mathrm{sin}(2m\alpha z)=\mathrm{sin}(2\sqrt{\eta }z/\mathrm{})`$ away from the axis. This describes a texture in the lateral spin density distribution. Here, $`\mathrm{}=(m\omega _0)^{1/2}`$ is the oscillator length and the parameter $`\eta =\frac{m}{\omega _0}\alpha ^2=\frac{2\epsilon _\mathrm{F}}{\omega _0}(\alpha /v_\mathrm{F})^2`$ is introduced for later use ($`\epsilon _\mathrm{F}`$ is the Fermi energy above the subband edge and $`v_\mathrm{F}`$ the Fermi velocity). In the mostly used compounds $`\alpha /v_\mathrm{F}<0.2`$  (in GaAs it is even much smaller), so that always $`\eta <0.1`$ in the lowest subband.

At $`k=0`$ different spin orientations are degenerate but when $`k0`$ Rashba splitting occurs, yielding

$$E_{k0s}=m^{}\alpha ^2+\frac{\omega _0}{2}+\frac{1}{2m^{}}(k+sm^{}\alpha ^{})^2,s=\pm 1$$
(4)
for the energies of the lowest subband. Expression (4) is obtained using the basis (3) after expanding Laguerre polynomials arising from band mixing by the spin-orbit term in powers of $`\eta `$. Up to the order $`𝒪(\eta ^{5/2})`$ one can incorporate the effect of subband mixing, described by $`\eta `$, into renormalized values for effective mass $`m^{}=m(1+8\eta ^2)`$ and Rashba parameter $`\alpha ^{}=\alpha (1\eta )`$. Note that a similar reasoning would be true also for other than parabolic dispersion relations where effects of inter-subband mixing could still be incorporated into renormalized values for the kinetic energy parameter and for $`\alpha `$. The dispersion relation (4) resembles the one obtained in two dimensions  but restricted to one dimensional $`k`$–space, cf. Fig. 1. For higher subbands the leading energy correction $`2n\eta `$, even in the strictly parabolic dispersion case, depends on the subband index $`n`$ which somewhat weakens the optimistic conjecture expressed by Datta and Das  that the proposed spin transistor would not loose sensitivity in multi-mode operation or at elevated temperatures or voltages. The Rashba phases in higher subbands differ slightly, $`\eta `$, from the phase acquired in the ground subband.

Without interactions the difference $`\mathrm{\hspace{0.25em}2}m^{}\alpha ^{}`$ in momenta between $`s=+1`$ and $`s=1`$ carriers, being independent of Fermi energy, is the origin for the Rashba precession of spins initially polarized e.g. along the $`+x`$–direction. After traversing the Rashba distance $`\lambda _\mathrm{R}=\pi /m^{}\alpha ^{}`$ the spin is reversed. Contrary to the observations , Rashba splitting, and therefore $`\lambda _\mathrm{R}`$, is independent of carrier density $`\mathrm{\hspace{0.25em}2}k_\mathrm{F}/\pi `$ in (4) at given interfacial electric field ; $`k_\mathrm{F}`$ is the Fermi momentum. Also, we note that the Fermi velocities of both spin components are equal.

The term $`sm^{}\alpha ^{}`$ inside the brackets in Eq. (4), shifting particle momentum, formally resembles a vector potential associated with the magnetic flux through an Aharanov–Bohm ring, leading there  to a persistent charge current in the ground state. Similarly, in (4) this term causes a non-zero persistent spin current $`J_\sigma `$ when periodic boundary conditions are imposed .

This similarity carries over to the TL–model for the quantum wire at low energies which in Bose variables reads

$$H=\underset{\nu =\rho ,\sigma }{}\frac{\pi }{4L}\left(v_{\nu \text{N}}N_\nu ^2+v_{\nu \text{J}}J_\nu ^2\right)+\left(\underset{q0}{}H_q\right)m^{}\alpha ^{}v_\mathrm{F}J_\sigma .$$
(5)
Here, $`N_\nu =N_{\nu \text{R}}+N_{\nu \text{L}}`$ and $`J_\nu =N_{\nu \text{R}}N_{\nu \text{L}}`$ denote particle numbers and currents, respectively, $`N_{\mathrm{L}/\mathrm{R}}`$ are the number of left/right going particles; $`L`$ is the wire length. Topological excitations of $`N_\nu `$ or $`J_\nu `$ as well as density excitations at momentum $`q`$, described by $`H_q`$, preserve charge ($`\nu =\rho `$) spin ($`\nu =\sigma `$) separation (as already mentioned, here we have to high accuracy equal velocities of both spins, cf. ; this is not the case, for example, in the presence of a Zeeman field , when spin up and spin down velocities differ at the Fermi energy due to the quadratic energy dispersion). The interaction is exactly included in the TL–model by renormalized values for its parameters such as the velocities $`v_{\nu \text{N}}`$ and $`v_{\nu \text{J}}`$ which now differ from $`v_\mathrm{F}`$. In (5) we have already omitted the backscattering in spin sector originating from the exchange of two electrons near two opposite Fermi points of opposite spins. Such a term would couple the topological sector with the charge and spin density excitations. For repulsive interactions back scattering is known to be irrelevant at low energies , also when persistent currents are present . Therefore, the $`H_q`$ will decouple, leaving the the correlation exponents $`K_\nu `$ ($`K_\nu 1`$ without interactions) and the power law decays of various correlation functions, that can be calculated within the TL–model, unaffected  by the Rashba coupling $`\alpha `$, cf. Eq. (7) below, contrary to the result obtained in Ref. . Spin-orbit coupling, and therefore the Rashba term (the last term in (5)), does not depend explicitly on the interaction.

Electrons injected at $`x=0`$ with spins polarized along the wire axis will perform Rashba precession. Most directly this can be monitored using the correlation function $`f(x)=\frac{1}{2}(\mathrm{\Psi }_{}(x)+\mathrm{\Psi }_{}(x))(\mathrm{\Psi }_{}^+(0)+\mathrm{\Psi }_{}^+(0))`$ where the Fermi operator $`\mathrm{\Psi }_s(x)\mathrm{e}^{\mathrm{i}\sqrt{\pi /2}(\varphi _\rho (x)+s\varphi _\sigma (x))}_{r=\pm }\mathrm{e}^{\mathrm{i}r\sqrt{\pi /2}(\theta _\rho (x)+s\theta _\sigma (x))}`$ is expressed through the momentum

$`\varphi _\nu (x)=\sqrt{{\displaystyle \frac{\pi }{2}}}J_\nu {\displaystyle \frac{x}{L}}+{\displaystyle \underset{q0}{}}\varphi _{\nu q}(x)`$

and the density like Bose fields

$`\theta _\nu (x)=\sqrt{{\displaystyle \frac{\pi }{2}}}N_\nu {\displaystyle \frac{x}{L}}+{\displaystyle \underset{q0}{}}\theta _{\nu q}(x)`$

in the usual way . The $`q0`$ components yield power law decay while the $`q=0`$ components are relevant for charge or spin stiffness and the persistent current. The result for the desired correlation function is

$`f(x)`$ $`=`$ $`{\displaystyle \frac{k_\mathrm{F}}{\pi }}|k_\mathrm{F}x|^{(K_\rho +1/K_\rho +K_\sigma +1/K_\sigma )/4}`$ (7)
$`\times \mathrm{sin}(|k_\mathrm{F}x|)\mathrm{cos}({\displaystyle \frac{\pi }{2}}j_\sigma x).`$
The square $`|f(L)|^2`$ is proportional to the probability for spin polarization in $`+x`$–direction at a distance $`L`$ from the source and, by a similar reasoning as in , to the current accepted by an ideal ferromagnetic drain which is polarized parallel to the source in $`+x`$–direction. A more sophisticated full transport calculation has been carried out recently  in the absence of spin-orbit coupling but would go beyond the scope of the present work. The period in the slowly oscillating last term in (7), $`j_\sigma J_\sigma /L=\frac{2}{\pi }m^{}\alpha ^{}v_\mathrm{F}/v_{\sigma \text{J}}`$, where the expectation value $`.`$ refers to the ground state of (5), determines the Rashba length

$$\lambda _\mathrm{R}=\frac{\pi }{m^{}\alpha ^{}}\frac{v_{\sigma \text{J}}}{v_\mathrm{F}}.$$
(8)
This quantity depends on the velocity of the spin current $`v_{\sigma \text{J}}`$ and is now altered compared to its value in the absence of interactions. We note, that in charge sector the corresponding velocity $`v_{\rho \text{J}}=v_\mathrm{F}`$ remains constant as a consequence of Galilei invariance of the quantum wire  at carrier densities much smaller than the inverse lattice constant of the underlying semiconductor lattice. Galilei invariance does not hold in spin sector where particles of opposite spins moving in opposite directions will experience some drag by the interaction. Therefore, $`v_{\sigma \text{J}}`$ differs from $`v_\mathrm{F}`$.

Rigorous relations hold  among the velocities $`v_{\nu \text{J}}=v_\nu K_\nu =v_{\nu \text{N}}K_\nu ^2`$ in TL–liquids, determined by the correlation exponents. Here, $`v_\nu `$ is the velocity of density excitations at small $`q`$. In many cases spin rotation invariance holds and fixes $`K_\sigma =1`$. Rashba coupling breaks this SU(2) invariance. Fortunately, the corresponding energy $`\mathrm{\hspace{0.25em}2}m^{}\alpha ^{}v_\mathrm{F}`$ is small compared to the Fermi energy (or to the typical Coulomb energy $`v_{\rho \text{N}}/v_\mathrm{F}`$) so that $`K_\sigma `$ will deviate from unity at most slightly and $`v_{\sigma \text{J}}v_\sigma v_{\sigma \text{N}}`$ to good accuracy. Much more important than the differences between those three velocities are their deviation from $`v_\mathrm{F}`$ occurring already in the SU(2) symmetric case as seen for example in the Hubbard chain . How long range interactions alter $`v_\sigma /v_\mathrm{F}`$ has been investigated perturbatively  and, recently, by extensive quantum Monte Carlo studies . Starting from values close to unity at high carrier densities, $`r_\mathrm{s}<0.5`$, $`v_\sigma /v_\mathrm{F}`$ can drop below 0.5 when $`r_\mathrm{s}1.4`$. Those densities are easily reached in present day quantum wires . In the limit of very small particle densities $`2k_\mathrm{F}/\pi `$ the spin velocity is expected to vanish like $`v_\sigma k_\mathrm{F}^2`$.

Thus, according to (8), Rashba precession is enhanced by repulsive interactions. Qualitatively, one might understand this result on the mean field level as a consequence of the Fock contribution precipitating repulsions between opposite rather than same spins. This enhances spin splitting and the nominal value of the Rashba parameter. In principle this argument applies also to situations with more than one subband occupied or even two dimensional heterostructures. Also in those cases we would expect that Coulomb interactions amplify Rashba precession as it has been confirmed in Ref. .

The amount of Rashba enhancement expected at given carrier density can be measured independently through the Zeeman spin susceptibility $`\chi =(^2E_0/B^2)/L=(L^2E_0/N_\sigma ^2)^1=2/\pi v_{\sigma \text{N}}`$ by monitoring the ‘exchange enhancement’ of the effective $`g`$–factor; $`E_0`$ is the ground state energy of (5). In InAs the magnetic energy $`B`$ translates quite accurately as 1 Ry per Tesla when taking the $`g`$–factor  as $`|g|=13`$. Assuming again $`K_\sigma 1`$ yields

$$\lambda _\mathrm{R}=2/m^{}\alpha ^{}v_\mathrm{F}\chi .$$
(9)
In conclusion, we have established a theory beyond describing Rashba precession as a single particle band structure effect. We have considered a quantum wire and the TL–model to incorporate interactions exactly. Increasing repulsion between carriers along with decreasing particle densities is found to reduce the Rashba length $`\lambda _\mathrm{R}`$ over which spins complete cycles as they move along the wire. Accordingly, the nominal value of the Rashba parameter increases, as determined by Shubnikov–de Haas measurements. This is demonstrated to be a consequence of decreasing spin velocities. The latter could be measured independently through the magnetic susceptibility w.r.t. a Zeeman field. Contrary to the relatively small influence of gates on the strength of the interfacial electric field this interaction induced contribution can explain variations of $`\lambda _\mathrm{R}`$ by a factor of 2. It would be valuable to experimentally separate the influence of the field strength from the carrier density by applying voltages independently to a front and a back gate . Moreover, the front gate could screen the long range part of the Coulomb interaction and thereby serve to vary the microscopic interaction strength.

I thank Charles Creffield, Dirk Grundler, Allan MacDonald, Ulrich Merkt, and, especially, Ulrich Zülicke for valuable discussions, and M. E. Raikh for drawing my attention to the work of Ref. .

FIGURE 1

Fig. 1: Energy dispersion in the lowest spin split subband of a quantum wire with parabolic confining potential. On the wire axis the spins $`s=,`$ are polarized in the plane of the heterostructure.


format: text

text:
# Resumming the Light Hemisphere Mass and Narrow Jet Broadening distributions in e+e- annihilation**footnote *Work supported in part by the UK Particle Physics and Astronomy Research Council and by the EU Fourth Framework Programme ‘Training and Mobility of Researchers’, Network ‘Quantum Chromodynamics and the Deep Structure of Elementary Particles’, contract FMRX-CT98-0194 (DG 12 - MIHT).

## 1 Introduction

Over the last decade, event shape distributions in e+e- annihilation have provided some of the most precise tests of QCD. Most of the attention has focussed on three-jet event shape variables such as thrust, wide jet broadening, heavy hemisphere mass, C parameter etc.. These variables may be generated by single gluon emission from the underlying $`q\overline{q}`$ pair in the e+e- annihilation process which first contributes at $`𝒪(\alpha _s)`$. The analysis of such events is rather sophisticated. In the first instance, fixed order QCD perturbation theory  is employed (at next-to-leading order (NLO)). However, when the observable $`O`$ is small and $`L=\mathrm{ln}(1/O)`$ is large the distribution suffers large logarithmic contributions of the form $`𝒪(\alpha _s^nL^{2n1})`$ that render the fixed order perturbation theory insufficient. In many cases the dominant infrared logarithms can be resummed using the coherent branching formalism . Furthermore, significant power suppressed effects are present and have been phenomenologically studied . To accurately describe the data, all of these contributions – fixed order, infrared resummation and power correction – are found to be necessary. Perhaps the most important measurement extracted from three-jet event shape data is the determination of the strong coupling constant.

Four-jet event shape observables also contain useful information about QCD. They are more sensitive to the triple gluon vertex and therefore the true gauge structure of QCD  and also to the presence of other light coloured particles such as the gluino  that can be pair produced by gluon splitting. However, these variables have received much less attention partly because they are suppressed by an additional power of $`\alpha _s`$ requiring a second gluon to be radiated but also because the theoretical description is much less developed. Recently however, four separate general purpose Monte Carlo programs have been developed to estimate the next-to-leading order, $`𝒪(\alpha _s^3)`$, corrections to four-jet event shapes: MENLO PARC , DEBRECEN  and MERCUTIO  employing on the one-loop helicity amplitudes for $`e^+e^{}4`$ partons  and EERAD2  based on the interference of the one-loop matrix element with tree level . For most four-jet event shape variables, the situation is the same as for the three-jet event shape variables; the NLO corrections are very large and, for renormalisation scales of the order of the centre-of-mass energy $`Q`$, still undershoot the data significantly . This points at the presence of large infrared effects as well as significant power corrections. With the exception of the four-jet rate in certain schemes  and near-to-planar three-jet events , the issue of resumming infrared logarithms for specific four-jet event shape variables has not yet been addressed. Similarly, power suppressed effects have only been studied for the out-of-event plane momentum distribution . It is the goal of this Letter to resum the leading $`𝒪(\alpha _s^nL^{2n1})`$ and next-to-leading $`𝒪(\alpha _s^nL^{2n2})`$ infrared logarithms to all orders in the coupling constant for the light hemisphere mass $`\rho _L`$ and narrow jet broadening $`B_N`$ distributions.

To be more precise let us first define these variables. We first separate the event at centre-of-mass energy $`Q=\sqrt{s}`$ into two hemispheres $`H_1`$, $`H_2`$ divided by the plane normal to the thrust axis $`𝐧_T`$. Particles that satisfy $`𝐩_i.𝐧_T>0`$ are assigned to hemisphere $`H_1`$, while all other particles are in $`H_2`$. Jet broadening measures the summed scalar momentum transverse to the thrust axis in one of the hemispheres while the hemisphere mass is the invariant mass of the hemisphere,

$`B_N`$ $`=`$ $`\underset{i=1,2}{\mathrm{min}}{\displaystyle \frac{_{𝐩_kH_i}|𝐩_k\times 𝐧_T|}{2_k|𝐩_k|}}`$ (1)
$`\rho _L`$ $`=`$ $`{\displaystyle \frac{1}{s}}\underset{i=1,2}{\mathrm{min}}\left({\displaystyle \underset{𝐩_kH_i}{}}p_k\right)^2.`$ (2)
Note that this definition of the light hemisphere mass is the common modification of the original variables suggested by Clavelli . These four-jet event shape variables are intimately connected to their three-jet event shape counterparts, the wide jet broadening and heavy hemisphere mass,

$`B_W`$ $`=`$ $`\underset{i=1,2}{\mathrm{max}}{\displaystyle \frac{_{𝐩_kH_i}|𝐩_k\times 𝐧_T|}{2_k|𝐩_k|}}`$ (3)
$`\rho _H`$ $`=`$ $`{\displaystyle \frac{1}{s}}\underset{i=1,2}{\mathrm{max}}\left({\displaystyle \underset{𝐩_kH_i}{}}p_k\right)^2`$ (4)
that have the property of exponentiation. That is to say that the fraction of events where the observable $`O=B_W`$ or $`O=\rho _H`$ has a value less than $`O`$ obeys the exponentiated form,

$`R(O,\alpha _s(Q^2))`$ $`=`$ $`{\displaystyle _0^O}{\displaystyle \frac{1}{\sigma }}{\displaystyle \frac{d\sigma }{dO}}𝑑O`$ (5)
$`=`$ $`C(\alpha _s(Q^2))\mathrm{\Sigma }(O,\alpha _s(Q^2))+D(O,\alpha _s(Q^2))`$ (6)
where

$`C(\alpha _s(Q^2))`$ $`=`$ $`1+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}C_n\alpha _s^n`$ (7)
$`\mathrm{ln}\left(\mathrm{\Sigma }(O,\alpha _s(Q^2))\right)`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{m=1}{\overset{n+1}{}}}G_{nm}\alpha _s^nL^m`$ (8)
$`=`$ $`Lg_1(\alpha _sL)+g_2(\alpha _sL)+\alpha _sg_3(\alpha _sL)+\mathrm{}`$
$`D(O,\alpha _s(Q^2))`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}D_n\alpha _s^n.`$ (9)
Here $`C_n`$ and $`G_{nm}`$ are constants and the perturbatively calculable coefficients $`D_n0`$ as $`O0`$. Knowledge of $`G_{nn+1}`$ (or equivalently $`g_1`$) together with $`G_{11}`$ allows resummation of terms in $`\mathrm{\Sigma }`$ down to $`𝒪(\alpha _s^nL^{2n1})`$. Partial inclusion of other subleading logarithms may be accomplished by retaining all knowledge of $`g_2`$. The accuracy of $`R`$ also depends on the number of terms known in $`C(\alpha _s(Q^2)`$, with each known term giving information about two towers of logarithms. For example, $`C=1`$ resums the $`𝒪(\alpha _s^nL^{2n})`$ and $`𝒪(\alpha _s^nL^{2n1})`$ terms in $`R`$, $`C_1`$ resums the $`𝒪(\alpha _s^nL^{2n2})`$ and $`𝒪(\alpha _s^nL^{2n3})`$ terms and so on.

## 2 Coherent branching

The coherent branching formalism allows the resummation of soft and collinear logarithms due to the emission of gluons from a hard parton. As a specific example, let us consider the jet mass distribution $`J^a(Q,k^2)`$ as the probability of producing a final state jet with invariant mass $`k^2`$ from a parent parton $`a`$ produced in a hard process at the scale $`Q^2`$. For an initial quark this is

$`J^q(Q^2,k^2)=\delta (k^2)+{\displaystyle \underset{0}{\overset{Q^2}{}}}{\displaystyle \frac{\mathrm{d}\stackrel{~}{q}^2}{\stackrel{~}{q}^2}}{\displaystyle \underset{0}{\overset{1}{}}}\mathrm{d}z\mathrm{\Theta }(z^2(1z)^2\stackrel{~}{q}^2Q_0^2)P^{qq}[\alpha _s(z^2(1z)^2)\stackrel{~}{q}^2),z]`$ (10)
$`\times `$ $`[{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}\mathrm{d}q^2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}\mathrm{d}k_{}^{}{}_{}{}^{2}\delta (k^2z(1z)\stackrel{~}{q}^2{\displaystyle \frac{k_{}^{}{}_{}{}^{2}}{z}}{\displaystyle \frac{q^2}{1z}})J^q(z^2\stackrel{~}{q}^2,k_{}^{}{}_{}{}^{2})J^g((1z)^2\stackrel{~}{q}^2,q^2)`$
$`J^q(\stackrel{~}{q}^2,k^2)]`$
normalised such that

$$_0^{\mathrm{}}𝑑k^2J^a(Q^2,k^2)=1$$
(11)
and where the next-to-leading order $`qqg`$ splitting kernel in the $`\overline{MS}`$ scheme with $`N_f`$ flavours is

$$P_{qq}[\alpha _s,z]=\frac{\alpha _s}{2\pi }C_F\frac{1+z^2}{1z}\left(1+\frac{\alpha _s}{2\pi }K\right)+\mathrm{},$$
(12)
where

$$K=C_A(\frac{67}{18}\frac{\pi ^2}{6})\frac{5}{9}N{}_{f}{}^{}.$$
(13)
Eq. (10) has a simple physical interpretation. The first term is the possibility that the originating parton does not emit any radiation. The transverse momentum of the parton is therefore unchanged. Alternatively, the quark may branch into a quark and a gluon subject to the phase space constraints of two body decay and which subsequently undergo further emissions. This is described by the term proportional to $`J^qJ^g`$. The last term is due to virtual corrections and ensures that soft and collinear singularities are regularised. A similar equation holds for $`J^g`$ and involves the $`ggg`$ and $`gq\overline{q}`$ splitting kernels. Solving the integral equation for $`J^a`$ is equivalent to resumming the infrared logarithms. In particular, the probability that an isolated parton $`a`$ forms a jet with mass less than $`\rho Q^2`$ is given by

$$\mathrm{\Sigma }^a(\rho ,\alpha _s(Q^2))=_0^{\rho Q^2}𝑑k^2J^a(Q^2,k^2).$$
(14)
In Ref. , $`\mathrm{ln}(\mathrm{\Sigma }_H^a)`$ is solved to next-to-leading logarithmic accuracy.

Similarly we can define the function $`T^a(Q,𝐤_t;p_t)`$  which describes the distribution of the summed scalar transverse momentum $`p_t`$ in a jet of parton type $`a`$ produced with vector transverse momentum $`𝐤_t`$, at scale Q. The structure is identical to Eq. (10).

## 3 The probabilistic interpretation

We can apply the coherent branching formalism to event shapes in the following schematic way illustrated in Fig. 1. The underlying configuration in e+e- annihilation is the production of a quark-antiquark pair aligned with the thrust axis (Fig. 1(a)). Each parton then undergoes soft and collinear gluon emission (denoted by the open cone at the head of the parton). This contribution describes small angle and soft emission accurately in two-jet-like events when $`B_W`$ and $`\rho _H`$ are small and gives rise to the exponentiated first term in Eq. (6). However, it does not describe the possibility of wide angle gluon emission shown in Fig. 1(b) where the gluon is the hardest parton. This three-jet-like matrix element correction is not logarithmically enhanced at $`𝒪(\alpha _s)`$ and first contributes at next-to-next-to-leading logarithmic accuracy and can therefore be neglected.

Let us first focus on the hemisphere masses. To next-to-leading logarithmic accuracy, the fraction of events with heavy hemisphere mass less than $`\rho _H`$ i.e. $`k_1^2<\rho _HQ^2`$ and $`k_2^2<\rho _HQ^2`$ is given by the two jet contribution

$`R_H(\rho _H,\alpha _s(Q^2))\underset{\rho _H1}{=}`$ (15)
$`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_1^2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_2^2J^q(Q^2,k_1^2)J^q(Q^2,k_2^2)\mathrm{\Theta }(\rho _HQ^2k_1^2)\mathrm{\Theta }(k_1^2k_2^2)`$
$`+{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_1^2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_2^2J^q(Q^2,k_1^2)J^q(Q^2,k_2^2)\mathrm{\Theta }(\rho _HQ^2k_2^2)\mathrm{\Theta }(k_2^2k_1^2),`$
and where the constraint that $`\rho _H1`$ has suppressed the three and more jet contributions. Following the steps of Ref. , we can rewrite this formula using the phase space restrictions as

$$R_H(\rho _H,\alpha _s(Q^2))\underset{\rho _H1}{=}\left[\mathrm{\Sigma }_H^q(\rho _H,\alpha _s(Q^2))\right]^2.$$
(16)
where

$$\mathrm{\Sigma }_H^q(\rho _H,\alpha _s(Q^2))=\underset{0}{\overset{\rho _HQ^2}{}}dk^2J^q(Q^2,k^2).$$
(17)
The fraction of events with light hemisphere mass less than $`\rho _L`$ also receives contributions from the two-jet configuration when $`k_1^2<\rho _LQ^2`$ and $`k_2^2>k_1^2`$ and vice versa. Altogether we have,

$`R_L(\rho _L,\alpha _s(Q^2))\underset{\rho _L1}{=}`$ (18)
$`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_1^2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_2^2J^q(Q^2,k_1^2)J^q(Q^2,k_2^2)\mathrm{\Theta }(\rho _LQ^2k_1^2)\mathrm{\Theta }(k_2^2k_1^2)`$
$`+{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_1^2{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dk_2^2J^q(Q^2,k_1^2)J^q(Q^2,k_2^2)\mathrm{\Theta }(\rho _LQ^2k_2^2)\mathrm{\Theta }(k_1^2k_2^2).`$
Simplifying the phase space constraints, we find

$$R_L(\rho _L,\alpha _s(Q^2))\underset{\rho _L1}{=}2\mathrm{\Sigma }_H^q(\rho _L,\alpha _s(Q^2))\left[\mathrm{\Sigma }_H^q(\rho _L,\alpha _s(Q^2))\right]^2.$$
(19)
The functions that resum the logarithms for the light hemisphere mass are the same as those that resum the logarithms for the heavy hemisphere mass. Now however, exponentiation in its purest form is spoiled because the final result is a sum of terms.

The analysis for the narrow jet broadening proceeds in the same way. We find that $`R_N`$, the probability of finding an event with a light hemisphere mass less than $`B_N`$, is given by

$$R_N(B_N,\alpha _s(Q^2))\underset{B_N1}{=}2\mathrm{\Sigma }_W^q(B_N,\alpha _s(Q^2))\left[\mathrm{\Sigma }_W^q(B_N,\alpha _s(Q^2))\right]^2$$
(20)
where the probability of obtaining a jet with summed scalar transverse momentum $`p_t`$ with respect to the jet axis less than $`2BQ`$ starting from a parton of type $`a`$, $`\mathrm{\Sigma }_W^a`$, is given by,

$$\mathrm{\Sigma }_W^a(B,\alpha _s(Q^2)=_0^{2BQ}T^a(Q,\mathrm{𝟎},p_t)dp_t.$$
(21)
## 4 All-orders resummation of large logarithms

In this section we discuss the all-orders resummation of leading logarithms $`𝒪(\alpha _s^nL^{2n})`$ and next-to-leading $`𝒪(\alpha _s^nL^{2n1})`$ logarithms to all orders in the coupling constant. From eqs. (19) and (20), we see that to determine $`R_L`$ and $`R_N`$ requires knowledge of $`\mathrm{\Sigma }_H^a`$ and $`\mathrm{\Sigma }_W^a`$ respectively. Both of these functions have the exponentiated form of Eq. (6) and the corresponding functions $`g_1(\alpha _sL)`$ and $`g_2(\alpha _sL)`$ determined.

Explicit expressions for $`\mathrm{\Sigma }_H`$ valid to this order are given in  and, introducing the renormalisation scale dependence in the standard manner and dropping the parton index, we reproduce them here for illustrative purposes,

$`\mathrm{\Sigma }_H(\rho ,\alpha _s(\mu ^2),\frac{Q^2}{\mu ^2})`$ $`=`$ $`{\displaystyle \frac{\mathrm{exp}[(\alpha _s(\mu ^2),L)]}{\mathrm{\Gamma }[1𝒮(\alpha _s(\mu ^2),L)]}}`$ (22)
$`=`$ $`{\displaystyle \frac{\mathrm{exp}[Lf_1(x)+f_2(x)+x^2f_1^{}(x)\mathrm{ln}(\mu ^2/Q^2)]}{\mathrm{\Gamma }[1f_1(x)xf_1^{}(x)]}}+𝒪(\alpha _s^nL^{n1})`$
where

$$L=\mathrm{ln}(1/\rho ),x=\beta _0\alpha _s(\mu ^2)L.$$
(23)
The functions $`f_1`$, $`f_2`$ and $`f_1^{}`$ are

$`f_1(x)`$ $`=`$ $`{\displaystyle \frac{A^{(1)}}{2\pi \beta _0x}}\left[(12x)\mathrm{ln}(12x)2(1x)\mathrm{ln}(1x)\right],`$ (24)
$`f_2(x)`$ $`=`$ $`{\displaystyle \frac{A^{(2)}}{2\pi ^2\beta _0^2}}\left[2\mathrm{ln}(1x)\mathrm{ln}(12x)\right]`$
$`+{\displaystyle \frac{B^{(1)}}{2\pi \beta _0}}\mathrm{ln}(1x){\displaystyle \frac{A^{(1)}\gamma _E}{\pi \beta _0}}\left[\mathrm{ln}(1x)\mathrm{ln}(12x)\right]`$
$`{\displaystyle \frac{A^{(1)}\beta _1}{2\pi \beta _0^3}}\left[\mathrm{ln}(12x)2\mathrm{ln}(1x)+{\displaystyle \frac{1}{2}}\mathrm{ln}^2(12x)\mathrm{ln}^2(1x)\right],`$
$`f_1^{}(x)`$ $`=`$ $`{\displaystyle \frac{A^{(1)}}{2\pi \beta _0x^2}}\left[\mathrm{ln}(12x)2\mathrm{ln}(1x)\right],`$ (26)
with

$$\beta _0=\frac{11C_A2N_f}{12\pi },\beta _1=\frac{17C_{A}^{}{}_{}{}^{2}5C_AN{}_{f}{}^{}3C_FN_f}{24\pi ^2},$$
(27)
For quarks,

$$A^{(1)}=C_F,A^{(2)}=\frac{1}{2}C_FK,B^{(1)}=\frac{3}{2}C_F,$$
(28)
with $`K`$ given by Eq. (13).

Altogether Eqs. (22) to (28) are sufficient to determine $`R_N`$ of Eq. (20) to next-to-leading logarithmic accuracy. The expression for $`R_N`$ therefore correctly sums to all orders in the strong coupling, $`\alpha _s`$, only the leading two towers of large logarithms from $`𝒪(\alpha _s^nL^{2n})`$ down to $`𝒪(\alpha _s^nL^{2n1})`$.

The resummed formula (19) does not include the $`𝒪(\alpha _s^2L^2)`$ or $`𝒪(\alpha _s^2L)`$ terms (just as the analogous formulae for resumming three jet variables do not include the $`𝒪(\alpha _s^2L^2)`$ or $`𝒪(\alpha _s^2L)`$ terms) present in the lowest order perturbative coefficient. Similarly, it does not produce the $`𝒪(\alpha _s^3L^4)`$ to $`𝒪(\alpha _s^3L)`$ terms that occur in the next-to-leading order perturbative coefficient. The perturbative calculation provides the $`\alpha _s^2`$ and $`\alpha _s^3`$ contributions exactly and therefore the most significant omitted term is $`𝒪(\alpha _s^4L^6)`$, see Fig. 2.

Precisely the same discussion applies to the narrow jet broadening. Using the coherent branching formalism the functions $`g_1(\alpha _sL)`$ and $`g_2(\alpha _sL)`$ have been determined for $`\mathrm{\Sigma }_W`$ and Catani et al (CTW) have provided analogous expressions for $`\mathrm{\Sigma }_W`$ that are given in . However, in doing so certain simplifying approximations concerning the recoil transverse momentum have been made. Dokshitzer and collaborators  have found that treating the quark recoil more carefully causes the CTW result for $`\mathrm{\Sigma }_W`$ to be adjusted by a multiplicative factor. This modifies $`g_2(\alpha _sL)`$ but leaves $`G_{11}`$ unchanged. The same correction is in principle needed for the narrow jet broadening but is relevant beyond the next-to-leading logarithmic accuracy. Nevertheless, we chose to apply the correction to the CTW result. We do not give the final form for $`\mathrm{\Sigma }_W`$ here, but instead refer the interested reader to Ref. . Inserting the recoil-corrected form for $`\mathrm{\Sigma }_W`$ in the resummed expression for $`R_N`$ (20) again allows resummation of the leading and next-to-leading logarithms as shown in Fig. 2.

## 5 Numerical results

As usual, the resummed result contains part of the fixed order perturbative contribution and the overlap must be removed by matching. This is done by expanding the resummed result as a series in the strong coupling constant and explicitly removing the terms corresponding to the fixed order calculations. This can be achieved in several ways of which $`R`$ matching and $`\mathrm{ln}(R)`$ matching are the most common. In the $`R`$ matching scheme, the coefficients of each of the unsummed logarithms present in the fixed order perturbative coefficients must be numerically extracted. For the four-jet event shape observables discussed here, this corresponds to determining the coefficient of $`\alpha _s^2L`$ from the lowest order perturbative contribution and the coefficients of $`\alpha _s^3L^2`$ and $`\alpha _s^3L`$ from the next-to-leading order contribution. This is impractical. However, in the more commonly used $`\mathrm{ln}(R)`$ matching scheme, it is assumed that the fixed order result exponentiates and therefore it is not necessary to make this extraction because any logarithmic terms remaining after subtracting the overlap from the fixed order contribution are exponentially suppressed. Of course, for the four-jet variables considered here, $`R`$ does not exponentiate. Nevertheless, in the $`\mathrm{ln}(R)`$ matching procedure, any remaining logarithmic terms are still exponentially suppressed. We therefore employ the $`\mathrm{ln}(R)`$ matching procedure.

To be more precise, the resummed result has the form

$$R_{resummed}=𝒞_12\mathrm{\Sigma }𝒞_2\mathrm{\Sigma }^2+𝒪(\alpha _s),$$
(29)
where $`𝒞_i`$ are power series in $`\alpha _s`$ and $`𝒪(\alpha _s)`$ represents the contribution from three-jet configurations. The $`𝒪(\alpha _s)`$ term of $`𝒞_1`$ and $`𝒞_2`$ are related by the constraint that there can be no logarithmic contributions at $`𝒪(\alpha _s)`$ \- the expansion of $`R`$ containing logarithms starts at $`𝒪(\alpha _s^2)`$. However, to the next-to-leading logarithmic accuracy relevant for this paper, we can set $`𝒞_1=𝒞_2=1`$. The $`\mathrm{ln}R`$ matching formula then dictates that

$$\mathrm{ln}(R)=\mathrm{ln}(R_{resummed})+\mathrm{ln}(R_{fixed})\mathrm{ln}(R_{expanded}),$$
(30)
where $`\mathrm{ln}R_{fixed}`$ and $`\mathrm{ln}R_{expanded}`$ are obtained by expanding the fixed order result and the $`\mathrm{ln}R_{resummed}`$ through to $`𝒪(\alpha _s^3)`$ respectively. Using the expansion of $`\mathrm{ln}(\mathrm{\Sigma })`$ as a series,

$$\mathrm{ln}(\mathrm{\Sigma })=\underset{n=1}{\overset{\mathrm{}}{}}\underset{m=1}{\overset{n+1}{}}G_{nm}\alpha _s^nL^m,$$

and with $`𝒞_1=𝒞_2=1`$, we see that,

$`\mathrm{ln}(R_{expanded})`$ $`=`$ $`(G_{11}L+G_{12}L^2)^2\alpha _s^2`$ (31)
$`2(G_{11}L+G_{12}L^2)(G_{21}L+G_{22}L^2+G_{23}L^3+{\displaystyle \frac{1}{2}}(G_{11}L+G_{12}L^2)^2)\alpha _s^3`$
$`+𝒪(\alpha _s^4).`$
At next-to-leading logarithmic accuracy, logarithmic terms from the fixed order contribution remain in $`\mathrm{ln}(R)`$ given by Eq. (30), however, they are exponentially suppressed when translated back to $`R`$.

We expect that at large values of the observable $`O`$ the resummed result is dominated by the fixed order calculation. However, the resummed expressions (19) and (20) valid in the small $`\rho _L`$ and $`B_N`$ limits do not contain information about the kinematic endpoints of the distributions. To ensure that the resummed result vanishes at the endpoint we make the substitution

$$\frac{1}{O}\frac{1}{O}\frac{1}{O^{\mathrm{max}}}+1,$$
(32)
where $`O^{\mathrm{max}}`$ corresponds to the endpoint of the distribution at the accuracy of the fixed order calculation. We use $`\rho _L^{\mathrm{max}}=0.167`$ and $`B_N^{\mathrm{max}}=0.204`$.

Numerical results for the light hemisphere mass and for the narrow jet broadening are shown in figures 3 and 4 respectively. Throughout we set $`\mu =Q=M_Z`$ and use $`\alpha _s(M_Z)=0.118`$ corresponding to the current world average. The next-to-leading order result which diverges at small values of the event shapes is taken from  and is evaluated at $`\mu =Q`$.

Figures 3 and 4 show that the resummations are extremely important for $`\rho _L<0.01`$ and $`B_N<0.02`$. Rather than the divergent fixed order prediction, we have the more physical resummed result that the probability of finding events with no radiation (very small values of $`\rho _L`$ and $`B_N`$) are vanishingly small. For the reference value of $`\mu =Q`$, the peak position of the $`\rho _L`$ distribution occurs at $`\rho _L=0.01`$ with a height of 0.33, while for $`B_N`$ the peak occurs at $`B_N=0.02`$ with a value of 0.48. At larger values, the resummation changes the NLO prediction by a more moderate amount indicating that uncalculated higher order corrections are under control. At very large values of $`O`$, the resummed and NLO predictions coincide because of the matching procedure of Eq. (32).

To illustrate the residual renormalisation scale dependence, we also show the effect of varying $`\mu `$ by a factor of 2 either side of the reference value $`\mu =Q`$. We see that around the peak region, different scale choices alter the prediction by $`\pm 10\%`$.

The infrared resummation significantly improves the perturbative prediction for the event shape observable. However, when comparing with experimental data we should be aware that important non-perturbative hadronisation corrections are present. The effect of hadronisation on the distribution is to shift the value of the observable away from the two-jet region,

$$OO+O_{\mathrm{NP}},$$
(33)
where the non-perturbative correction depends on the typical hadron scale $`𝒪(1\mathrm{GeV})`$ and is suppressed by a power of $`Q`$. In principle these power corrections can be estimated using the dispersive approach of Ref.  where a non-perturbative parameter $`\mu _I`$ is introduced to describe the running of $`\alpha _s`$ in the infrared region. For the associated three jet variables the non-perturbative corrections are typically estimated to be $`𝒪(1\mathrm{GeV}/Q)`$ for $`\rho _H`$  and $`𝒪(0.3\mathrm{GeV}\mathrm{ln}(1/B_W)/Q)`$ for $`B_W`$  and arise through the hadronisation of one of the two jets in the event. Because the four-jet event shapes are largely related to what happens in the second jet, we might expect that the hadronisation corrections are similar. To illustrate the potential effects of hadronisation in the four-jet event shapes, we just transfer these corrections directly so that,

$`\rho _L`$ $``$ $`\rho _L+{\displaystyle \frac{1\mathrm{GeV}}{Q}},`$ (34)
$`B_N`$ $``$ $`B_N+{\displaystyle \frac{0.3\mathrm{GeV}\mathrm{ln}(1/B_N)}{Q}}.`$ (35)
The simplified hadronisation correction applied to the resummed distributions for the light hemisphere mass and narrow jet broadening with $`\mu =Q`$ is shown in figures 5 and 6 respectively. To emphasize the dramatic effect the power correction has, we also show the uncorrected predictions. There are two effects. First the distribution is shifted to the right by an amount $`O_{\mathrm{NP}}`$ and second, the distribution is rescaled by a factor $`(O+O_{\mathrm{NP}})/O`$. In the region of the turnover where $`O`$ is of the same order as $`O_{\mathrm{NP}}`$ there is an enhancement of $`𝒪`$(100%). The hadronisation correction is smaller at larger values of $`O`$.

For comparison, we also show the charged hadron data collected by the DELPHI Collaboration  at the $`Z`$ resonance. We see remarkable agreement (strikingly so in view of the simplified hadronisation correction applied here). The only discrepancy is at very small values of $`O<O_{\mathrm{NP}}`$ where individual hadrons in the light/narrow hemisphere will significantly affect the value of $`O`$. We do not expect to successfully describe such events.

## 6 Conclusions

In this paper, we presented predictions for the light hemisphere mass and the narrow jet broadening distributions where the infrared logarithms have been resummed to next-to-leading logarithmic order using the coherent branching formalism. The resummed expressions do not exponentiate, but involve the difference of exponential factors so that the perturbative series starts at $`𝒪(\alpha _s^2)`$. The expressions presented in sections 3 and 4 resum the leading $`𝒪(\alpha _s^nL^{2n1})`$ and next-to-leading $`𝒪(\alpha _s^nL^{2n2})`$ infrared logarithms to all orders in the coupling constant. Taken in conjunction with the fixed order $`𝒪(\alpha _s^3)`$ non-logarithmic terms, the first neglected term is of $`𝒪(\alpha _s^4L^6)`$. The numerical results shown in Figs. 3 and 4 indicate that the resummation effects are sizeable at small values of the event shape parameter and that the resummation procedure significantly improves the perturbative prediction.

However, because important non-perturbative hadronisation corrections are present, resummation alone is insufficient to describe the experimental data. Applying simple power corrections similar to those obtained for the heavy hemisphere mass and the wide jet broadening gives good qualitative agreement with the available data from LEP. We anticipate that the improved theoretical description of the four-jet event shape distributions presented here can be combined with a more sophisticated hadronisation correction based on the dispersive approach of Ref. . The data from LEP can then be used to further test the structure of QCD in four-jet like events and extract values of the strong coupling constant.

## Acknowledgements

We thank Mike Seymour and David Summers for useful and stimulating discussions. SJB thanks the UK Particle Physics and Astronomy Research Council for research studentships. This work was supported in part by the EU Fourth Framework Programme ‘Training and Mobility of Researchers’, Network ‘Quantum Chromodynamics and the Deep Structure of Elementary Particles’, contract FMRX-CT98-0194 (DG-12-MIHT).


text:
# Methods of Determination of the Energy and Mass of Primary Cosmic Ray Particles at Extensive Air Shower Energies

## 1 Introduction

The cosmic ray (cr) energy spectrum extends from about 1 GeV to above $`10^{20}`$ eV. Over this wide range of energies the intensity drops by more than 30 orders of magnitude. Despite the enormous dynamic range covered, the spectrum appears rather structureless and can be well approximated by broken power-laws $`dN/dEE^\gamma `$. Up to energies of some $`10^{14}`$ eV the flux of particles is sufficiently high allowing measurements of their elemental distributions by high flying balloon or satellite experiments. Such studies have provided important implications for the origin and transport properties of cr’s in the interstellar medium. Two prominent examples are ratios of secondary to primary elements, such as the B/C-ratio, which are used to extract the average amount of matter cr-particles have traversed from their sources to the solar system, or are radioactive isotopes, e.g. 10Be or 26Al, which carry information about the average ‘age’ of cosmic rays. With many new complex experiments taking data or starting up in the near future and with a possibly new generation of balloons, this is a vital field of research and has been subject to separate talks presented at this Symposium.

Above a few times $`10^{15}`$ eV the flux drops to only one particle per m2 per year. This excludes any type of ‘direct observation’ even in the near future, at least if high statistics is required. Ironically, one of the most prominent features of the cr energy spectrum is the steepening of the slope from $`\gamma 2.7`$ to $`\gamma 3.1`$ at an energy just above some $`10^{15}`$ eV. This is known as the ‘knee’ of the spectrum. It was first deduced from observations of the shower size spectrum made by Kulikov and Khristianson et alin 1956  but it still remains unclear as to what is the cause of this spectral steepening. At an energy above $`10^{18}`$ eV the spectrum flattens again at what is called the ‘ankle’. Data currently exist, though with very poor statistics, up to $`310^{20}`$ eV and there seems to be no end to the energy spectrum .

The origin and acceleration mechanism of these, so called, ultra- and extremely high energy cosmic rays have been subject to debate for several decades. Mainly for reasons of the required power the dominant acceleration sites are generally believed to be shocks associated with supernova remnants. Naturally, this leads to a power law spectrum as is observed experimentally. Detailed examination suggests that this process is limited to $`E/Z\begin{array}{c}<\\ \end{array}10^{15}`$ eV. Curiously, this coincides well with the knee at $`E_{\mathrm{knee}}410^{15}`$ eV, indicating that the feature may be related to the upper limit of acceleration. The underlying picture of particle acceleration in magnetic field irregularities in the vicinity of strong shocks suggests the maximum energies of different elements to scale with their rigidity $`R=pc/Ze`$. This naturally would lead to an overabundance of heavy elements above the knee, a prediction to be proven by experiments. A change in the cr propagation with decreasing galactic containment at higher energies has also been considered. This rising leakage results in a steepening of the cr energy spectrum and again would lead to a similar scaling with the rigidity of particles but would in addition predict anisotropies in the arrival direction of cr’s with respect to the galactic plane. Besides such kind of ‘conventional’ source and propagation models  several other hypotheses have been discussed in the recent literature. These include the astrophysically motivated single source model of Erlykin and Wolfendale  trying to explain possible structures around the knee, as well as several particle physics motivated scenarios which try to explain the knee due to different kinds of cr-interactions. For example, photodisintegration at the source , interactions with gravitationally bound massive neutrinos , or sudden changes in the character of high-energy hadronic interactions during the development of extensive air shower (eas)  have been considered.

To constrain the SN acceleration model from the other proposed mechanisms precise measurements of the primary energy spectrum and particularly of the mass composition as a function of energy are needed. Significant progress has been made on these problems in recent years, but the situation is far from being clear.

The other target of great interest is the energy range around the Greisen-Zatsepin-Kuzmin (gzk) cut-off at $`E510^{19}`$ eV. Explanation of these particles requires the existence of extreme powerful sources within a distance of approximately 100 Mpc. Hot spots of radio galaxy lobes – if close enough – or topological defects from early epochs of the universe would be potential candidates. This topic has been addressed at this Symposium in some detail by Ostrowski, Zavrtanik, and others. Therefore, I will be very brief on experimental aspects relevant to measurements in this energy range.

## 2 Extensive Air Showers and Experimental Observables

An air shower is a cascade of particles generated by the interaction of a single high energy primary cosmic ray particle with the atmosphere. The secondary particles produced in each collision - in case of a primary hadron mostly charged and neutral pions - may either decay or interact with another nucleus, thereby multiplying the particles within an eas. After reaching a maximum ($`X_{max}`$) in the number of secondary photons, electrons, muons, and hadrons, the shower attenuates as more and more particles fall below the threshold for further particle production. A disk of relativistic particles extended over an area with a diameter of some tens of meters at $`10^{14}`$ eV to several kilometre at $`10^{20}`$ eV can then be observed at ground. This magnifying effect of the earth atmosphere allows to instrument only a very small portion of the eas area and to still reconstruct the major properties of the primary particles. It is a lucky coincidence that at the energy where direct detection of crs rays becomes impractical, the resulting air showers become big enough to be easily detectable at ground level.

Due to the nature of the involved hadronic and electromagnetic interactions and the different decay properties of particles, an eas has three components, electromagnetic, muonic, and hadronic. On average, a 1 PeV primary proton will produce about $`10^6`$ particles at sea-level, 80 % of which are photons, 18 % electrons, 1.7 % muons, and about 0.3 % hadrons. Neutrinos will also be produced by weak decays of particles, but they remain unseen by ordinary eas experiments. During their propagation through the atmosphere, relativistic charged particles will also produce Cherenkov light and, finally, excitations of the $`2P`$ and $`1N`$-band of N2 and N$`{}_{}{}^{+}{}_{2}{}^{}`$, respectively, will give rise to emission of fluorescence light.

Extracting the primary energy and mass from such measurements is not straightforward and a model must be adopted to relate the observed eas parameters to the properties of the primary particle. As shall be discussed below, the analysis is complicated by large fluctuations of eas observables and by the fact, that virtually all of the eas observables are sensitive to changes in the mass and energy of the primary particle so that a careful reconstruction of the energy spectrum requires knowledge of the chemical composition.

Ignoring these complications for a moment, some basic characteristics of eas, as observed in Fig. 1, can be deduced from very simple considerations. On average, the first interaction of the primary particle with a nucleus in the atmosphere - mostly nitrogen - will occur within its characteristic interaction length, $`\lambda _I1/\sigma `$. Due to the larger hadronic cross section of a Fe+N system as compared to p+N, the iron nucleus will - on average - interact higher in the atmosphere. Because of the strong absorption of the electromagnetic and hadronic shower components in the atmosphere, one expects

to detect fewer electrons and hadrons at ground in case of primary heavy nuclei. In a simplified picture, we furtheron may consider the primary iron nucleus as 56 nucleons each having 1/56th of the primary energy. This has two important consequences: i) the energy of produced particles within a single collision will be lower by about the same number, so that more pions and kaons have a chance to decay into muons before reinteracting with another air nucleus, ii) the multiplicity of produced particles within a single high energy collision roughly scales with the logarithm of the energy, so that 56 nucleons each with 1/56th of the energy, will produce a higher number of secondary particles, which again may decay into muons. As a result of the two effects one expects to observe more muons at ground in case of heavy primaries. Also, fluctuations in the number of particles will be smaller in case of heavy primary nuclei. Major experiments performing particle measurements at ground in the energy range of the knee include CASA-MIA , EAS-TOP , HEGRA , KASCADE , and Tibet-AS$`\gamma `$ , only the latter two of which currently take data. Most of the data in the gzk energy range are based on the Akeno Giant-Air-Shower-Experiment (AGASA) . Their primary observables are lateral particle density distributions at ground, $`\rho _{e,\mu ,h}(r)`$, and their integrals $`N_{e,\mu ,h}=_{r_1}^{r_2}2\pi r\rho _{e,\mu ,h}(r)𝑑r`$, yielding total (for extrapolations $`r_10`$ and $`r_2\mathrm{}`$) or truncated particle shower sizes. Truncated muon-sizes, $`N_\mu ^{\mathrm{tr}}`$, have been introduced by the KASCADE collaboration in order to avoid systematic uncertainties resulting from extrapolations of $`\rho _\mu (r)`$ to large distances not covered by the experiment.

Besides measuring particles at ground, experiments also aim at observing traces of the longitudinal shower development. A basic parameter here is the position of the shower maximum, $`X_{max}`$, which penetrates deeper into the atmosphere with increasing primary energy and decreasing primary mass. Experimentally, $`X_{max}`$ is often inferred from the lateral density distribution of Cherenkov photons, $`\rho _{Ch}`$. Experiments following this approach include CASA-BLANCA  and the HEGRA AIROBICC detectors . Also, air Cherenkov telescopes have been employed for sampling the shower maximum (HEGRA , DICE ), and more recently, $`X_{max}`$ has been reconstructed from muon tracks at ground by means of triangulation (HEGRA  and KASCADE (talk presented by C. Büttner at this Symposium)). A very elegant approach is the detection of air fluorescence light at large distances from the shower axis giving true information about the longitudinal development of the eas. The major drawback is the low light yield limiting its application to energies $`E\begin{array}{c}>\\ \end{array}10^{17}`$ eV. This technique pioneered by Fly’s Eye , is being employed by HIRes  and the Pierre Auger Observatory . A summary of the different experimental observables is given in table 2. Deep underground muon detectors are omitted here because of limited space and because of their only indirect relation to eas observations.

## 3 Techniques of Data Analysis and Selected Results

### 3.1 Single parameter methods

Until recently, most experiments have performed single parameter analyses using only one of the eas observables discussed above. The composition has then been inferred by comparing the mean value of the experimental observable to either different trial composition models or to mean values of simulated proton and iron primaries and interpolations to a hypothetical ‘mean mass’. Obviously, a significant piece of information available from the full distribution of shower observables is ignored in such analyses. The problems here are manifold, most importantly: i) the extracted mean mass is ambiguous, i.e. the mean values of the experimental distribution can be described by a single ‘mean mass’ component, by a mix of pure protons and iron, or by any other combination. Clearly, a ‘composition’ in a strict sense is not obtained. ii) There is no clue of how well the simulations describe the experimental data. Only if the mean value of the experimental observable falls outside the window spanned by

proton and iron simulations, doubts can be raised about the quality of the simulations (or apparatus). Furthermore, publications discussing single parameter results often do another simplification by comparing measured and simulated observables not at the level of the detector, but instead comparing quantities which are inferred from experimental data by means of eas simulations. A very prominent example of such kind of analyses are plots of $`X_{max}`$ as a function of energy as shown in Fig. 2. In this and many other cases, $`X_{max}`$ is based on measurements of the slope of the Cherenkov lateral distribution within about 120 m from the shower core. eas simulations using different primary particles are then performed to convert this experimental slope parameter to $`X_{max}`$. Finally, this biased and model dependent $`X_{max}`$-value is presented in such graphs and compared to $`X_{max}`$-values obtained directly from CORSIKA simulations. Clearly, comparing experimental and simulated slope parameters directly instead of introducing secondary quantities appears much more appropriate and avoids sources of systematic biases. Determination of $`E_{prim}`$ in experiments using open Cherenkov counters is mostly based on the Cherenkov light intensity observed at about 120 m from the shower core where influences of the unknown primary mass are minimal. However, a careful analysis of this procedure exhibits a residual mass dependence on the order of $`\pm 10`$ % in the knee energy range . Also, uncertainties to the energy scale by at least the same amount are caused by hadronic interaction models  but, again, are not considered in plots of the type of Fig. 2. Finally, the figure also nicely demonstrates the aforementioned inabilities to judge the quality of the different hadronic interactions models or to infer a true composition. Despite this criticism, the qualitative feature of the data – a somewhat lighter ‘mean mass’ towards the knee and a heavier one above – appears common to all interaction models.

With increasing computer power, eas simulations have advanced in quality and quantity, i.e. both the complexity to which details of hadronic interactions and propagation effects in the atmosphere are taken into account has been improved as well as the number of simulated events available for comparison. Therefore, analyses of eas data could progress from comparing only mean values to full distributions of observables, the benefits of which have been discussed already. Two examples are presented in Fig. 3. The left hand side of the top panel presents the event-by-event distribution of the inner Cherenkov-slope (previously used as a measure for $`X_{max}`$, see above) in an energy range $`810^{14}E210^{15}`$ eV  and the right hand side the distribution of $`\mathrm{log}N_\mu ^{\mathrm{tr}}/\mathrm{log}N_e`$ obtained from KASCADE for zenith angles $`18^{}`$ \- $`25^{}`$ and similar energies . Both data sets are compared to CORSIKA simulations of different primary particles using the QGSJET model. Note, that the left hand tails of the distributions are well accounted for by proton simulations while the right hand tails are accounted for by iron simulations. The KASCADE distributions cover more than three orders of magnitude and thus provide some level of confidence in the results of the simulations. Both data sets also require at least 3 different primary masses to account for the full shape of the experimental distributions. However, it is also observed that the Cherenkov data alone do not provide sufficient discrimination power to distinguish the nitrogen from the iron distribution. The bottom panel shows the extracted mass fractions as a function of energy. Clearly, a heavier composition is observed at energies above the knee. Also, the relative abundance of iron is fairly similar in both data sets but the major difference is a reversed proton and Helium abundance resulting in a somewhat lighter composition from KASCADE data. To resolve this differences, several cross checks of the two data sets appear expedient, e.g. tests of the zenith angle dependence, still higher statistics in simulations to verify the tails of the distributions, etc. Such work is in progress now.

### 3.2 Multi-component methods and non-parametric approaches

Multi-component detector installations can measure several eas parameters simultaneously on an event-by-event basis. Proper combinations can then be identified to provide an estimate of the primary energy and mass (see Table 1,2). The analysis of multivariate parameter distributions also needs to account for influences of fluctuations which may be different in different observables. Non-parametric Bayesian methods and neural network approaches are well suited for these purposes and they also specify the uncertainty of the results in a quantitative way . Simpler approaches, like the knn-method have also been employed , but they suffer from various deficiencies like dependencies of the results on the density and thereby on the statistics of simulated data points.

First comprehensive non-parametric approaches with promising results have been presented at this Symposium by M. Roth et al. In these techniques, each event is represented by an observation vector $`x=(N_\mu ^{\mathrm{tr}}`$, $`N_e`$, $`N_h\mathrm{})^T`$ of $`n`$ suitable eas observables. This vector serves as input to a procedure based on Bayesian or neural network decision rules by which the observed event is assigned to a given class of elemental groups, say p, O, and Fe. The first step in such an analysis is a calculation of likelihood (probability density) distributions $`p(x|\omega _i)`$ by means of Monte Carlo simulations using different primary particles. Here, $`p(x|\omega _i)`$ describes the probability to observe the multi-dimensional observation vector $`x`$ in an event belonging to a given class $`\omega _i\{`$p, O, Fe$`\}`$. The term ‘non-parametric’ indicates that the representations of the distributions (like probability density functions of Bayes classifiers or weights of neural networks) are no more specified by a-priori chosen functional forms, but are constructed through the analysis process by the given (simulated) data distributions themselves. Using the Bayes theorem one can then translate the likelihood $`p(x|\omega _i)`$ for finding an event $`x`$ in a given class $`\omega _i`$ to the actually required distribution $`p(\omega _i|x)`$ representing the probability of class $`\omega _i`$ being associated to a measured event $`x`$. In this step, an assumption is made about an a priori knowledge of the relative abundance of each class, a subject being handled by Bayesian inference procedures. If there is no further knowledge, the prior probabilities should be assumed to be equal (Bayes’ Postulate).

Analysing an experimentally measured event, a statistical decision on the primary particle type/energy is then to be made. The applied Bayesian approach of the statistical inference provides the Bayes optimal way of combining prior and experimental knowledge and the Bayes theorem specifies how such modification should be made . In simplest case, the Bayes optimal decision rule is to classify $`x`$ into class $`\omega _i`$, if $`p(\omega _i|x)>p(\omega _j|x)`$ for all classes $`\omega _j\omega _i`$. In this way one is also able to specify the actual quality of the decision by classification and misclassification matrices of the results. Unfortunately, because of limited computer time, analyses making use of all these pieces of important information have not yet been performed for high statistics data. However, good agreement to faster neural network approaches was reported in Ref. .

It is important to realize that the composition extracted that way yields a somewhat heavier composition as compared to the results from Fig. 3. This effect may possibly be caused by different treatments of correlations among the fluctuating eas observables. Furthermore, detailed investigations also show that within a given interaction model, different sets of eas observables lead to different mass compositions . Generally, inclusion of hadronic observables tends to shift the composition to a heavier one, a result already suggested by earlier studies of the KASCADE collaboration . These inconsistencies appear to be caused by deficiencies of the employed hadronic interaction models. Both kinds of effects are subject to further studies and, at present, any quantitative result on chemical composition may be subject to further changes.

Results on chemical composition are often presented in terms of the mean logarithmic mass. However, this quantity has some severe drawbacks as any information about individual abundances is lost. Furthermore, in absence of an unique prescription on how to calculate the average from abundances of different mass groups, uncertainties in $`\mathrm{ln}A`$ may be on the order of 0.2 and more. This has been source for some confusion in the past and can best be avoided by providing tables of reconstructed abundances of

mass groups (including correlated uncertainties).

A compilation of the all-particle energy spectrum reconstructed from various experiments is presented in figure 4. The agreement appears reasonable and deviations are mostly explained by uncertainties in the energy scale by up to 25 %, e.g. CASA MIA data  were shifted upwards in energy by 20 % to yield a better agreement to the other data sets. This is likely to be explained by the interaction model SIBYLL 1.6  employed by the authors of Ref.  but which has been proven to provide only a poor description of the experimental data . The lines correspond to a simultaneous fit of the electron and muon size spectra of KASCADE, assuming the all-particle spectrum to be described by a sum of proton and iron primaries . Interestingly, a knee is only reconstructed for the light component and no indication of a break is seen in the heavy component up to about $`10^{17}`$ eV. This important finding giving direct support to the picture of acceleration in magnetic fields (see above) will be target of future studies with improved experimental capabilities. For example, KASCADE and EAS-TOP have just started a common effort to install the EAS-TOP scintillators at the site of Forschungszentrum Karlsruhe providing a 12 times larger acceptance as compared to the original KASCADE experiment and still taking advantage of the multi-detector capabilities.

## 4 Summary and Outlook

Much progress has been made in reconstructing the parameters of primary cr-particles from eas observables. The advent of multiparameter measurements and realistic eas simulations now allow for much more detailed investigations taking account for eas fluctuations and detector effects. Different experiments agree fairly well on the all-particle energy spectrum but there is still some uncertainty in extracting the chemical composition. This may partly be explained by uncertainties of hadronic interaction models employed in eas simulations and there appears a demand for more tests of hadronic interaction models by means of accelerator and/or eas data. Other sources of systematic uncertainties appear to be caused by incomplete data analysis techniques and biases affecting mostly single parameter measurements. As discussions and presentations at this Symposium have shown, there is well founded hope to solve these imperfections thereby answering the question about the sources of crs and the origin of the knee already in the very near future.

#### Acknowledgement

I would like to express my gratitude to the organizing committee for their hospitality and inviting me to give this talk. During the preparation I greatly benefited from fruitful discussions with my colleagues from the KASCADE collaboration, with C. Pryke, J. Matthews, S. Swordy, and many others. I would also like to thank M. Giller for valuable comments to this manuscript.

## References