text:
# OFF-FORWARD PARTON DISTRIBUTIONS AND IMPACT PARAMETER DEPENDENT PARTON STRUCTURE

## 1 Introduction

Deeply virtual Compton scattering experiments in the Bjorken limit allow measuring generalized or off-forward parton distributions (OFPDs)

$`\overline{p}^+{\displaystyle \frac{dx^{}}{2\pi }p^{}|\overline{\psi }(\frac{x^{}}{2})\gamma ^+\psi (\frac{x^{}}{2})|pe^{ix\overline{p}^+x^{}}}`$ $`=`$ $`H(x,\xi ,t)\overline{u}(p^{})\gamma ^+u(p)`$
$`+E(x,\xi ,t)\overline{u}(p^{}){\displaystyle \frac{i\sigma ^{+\nu }\mathrm{\Delta }_\nu }{2M}}u(p),`$
$`\overline{p}^+{\displaystyle \frac{dx^{}}{2\pi }p^{}|\overline{\psi }(\frac{x^{}}{2})\gamma ^+\gamma _5\psi (\frac{x^{}}{2})|pe^{ix\overline{p}^+x^{}}}`$ $`=`$ $`\stackrel{~}{H}(x,\xi ,t)\overline{u}(p^{})\gamma ^+\gamma _5u(p)`$
$`+\stackrel{~}{E}(x,\xi ,t)\overline{u}(p^{}){\displaystyle \frac{\gamma _5\mathrm{\Delta }^+}{2M}}u(p),`$
where $`x^\pm =x^0\pm x^3`$ and $`p^+=p^0+p^3`$ refer to the usual light-cone components, $`\overline{p}=\frac{1}{2}\left(p+p^{}\right)`$, $`\mathrm{\Delta }=pp^{}`$, and $`t\mathrm{\Delta }^2`$. The “off-forwardness” (or skewedness) in Eqs. (1,1) is defined as $`\xi \frac{\mathrm{\Delta }^+}{p^+}`$. From the point of view of parton physics in the infinite momentum frame, these OFPDs have the physical meaning of the amplitude for the process that a quark is taken out of the nucleon with longitudinal momentum fraction $`x`$ and then inserted back into the nucleon with a four momentum transfer $`\mathrm{\Delta }^\mu `$ . OFPDs play dual roles and in a certain sense they interpolate between form factors and conventional parton distribution functions (PDFs) : for $`\xi =t=0`$ one recovers conventional PDFs, i.e. longitudinal momentum distributions in the infinite momentum frame (IMF), while when one integrates $`H(x,\xi ,t)`$ over $`x`$, one obtains a form factor, i.e. the Fourier transform of a position space density (in the Breit frame!). One of the new physics insights that one can learn from these OFPDs is the angular momentum distribution . Others include meson distribution amplitudes. 111For a discussion of this connection in the context of $`QCD_{1+1}`$ see Ref. . Beyond that it is not yet entirely clear what kind of new physics insights one gains from studying these generalized parton distributions. The main reason, why the physical interpretation of OFPDs is still somewhat obscure is due to the fact that the initial and final state in Eq. (1) are not the same and therefore, in general, these OFPDs cannot be interpreted as a ‘density’ but rather their physical significance is that of a probability amplitude.

In this note, we will discuss the limit $`\xi 0`$, but $`t0`$, where

$$H(x,t)H(x,\xi =0,t)\text{and}\stackrel{~}{H}(x,t)\stackrel{~}{H}(x,\xi =0,t)$$
(3)
do have a simple interpretation in terms of a density, namely as the Fourier transform of impact parameter dependent parton distributions w.r.t. the impact parameter, i.e.

$`H(x,\stackrel{}{\mathrm{\Delta }}_{}^2)`$ $`=`$ $`{\displaystyle d^2b_{}q(x,\stackrel{}{b}_{})e^{i\stackrel{}{b}_{}\stackrel{}{\mathrm{\Delta }}_{}}}`$
$`\stackrel{~}{H}(x,\stackrel{}{\mathrm{\Delta }}_{}^2)`$ $`=`$ $`{\displaystyle d^2b_{}\mathrm{\Delta }q(x,\stackrel{}{b}_{})e^{i\stackrel{}{b}_{}\stackrel{}{\mathrm{\Delta }}_{}}}`$ (4)
## 2 Impact parameter dependent PDFs and OFPDs

PDFs are usually defined as matrix elements between plane wave states, which extend throughout space. Therefore, before we can introduce the notion of impact parameter dependent PDFs, we need to define localized wave packets. For our purposes, it is most suitable to consider a wave packet $`|\mathrm{\Psi }`$ which is chosen such that it has a sharp longitudinal momentum $`p_z`$, but whose position is a localized wave packet in the transverse direction

$$|\mathrm{\Psi }=\frac{d^2p_{}}{\sqrt{2E_\stackrel{}{p}(2\pi )^2}}\mathrm{\Psi }(\stackrel{}{p}_{})|\stackrel{}{p}.$$
(5)
Although this definition is completely general, what we have in mind for the states $`|p`$ are for example nucleon states, which are of course extended objects, i.e. Eq. (5) describes wave packets of particles that are themselves already extended objects. We will come back to this point further below. Clearly,

$$F_\mathrm{\Psi }(x,\stackrel{}{b}_{})\frac{dx^{}}{2\pi }\mathrm{\Psi }|\overline{\psi }(\frac{x^{}}{2},\stackrel{}{b}_{})\gamma ^+\psi (\frac{x^{}}{2},\stackrel{}{b}_{})|\mathrm{\Psi }e^{ix\overline{p}^+x^{}},$$
(6)
describes the probability to find partons with momentum fraction $`x`$ at transverse (position) coordinate $`\stackrel{}{b}_{}`$ in this wave packet. Note that in Eq. (6) we have implicitly assumed that we work in light-cone gauge $`A^+=0`$. In any other gauge, one needs to insert a (straight line) ‘gauge string’ connecting the points $`(\frac{x^{}}{2},\stackrel{}{b}_{})`$ and $`(\frac{x^{}}{2},\stackrel{}{b}_{})`$ in order to render Eq. (6) manifestly gauge invariant.

What we will show in the following is that, for a suitably localized wave packet, $`F_\mathrm{\Psi }(x,\stackrel{}{b}_{})`$ can be related to OFPDs with $`\xi =0`$. Using Eq. (5), one finds

$`f_\mathrm{\Psi }(x,\stackrel{}{\mathrm{\Delta }}_{})`$ $``$ $`{\displaystyle d^2\mathrm{\Delta }_{}e^{i\stackrel{}{\mathrm{\Delta }}_{}\stackrel{}{b}_{}}F_\mathrm{\Psi }(x,\stackrel{}{b}_{})}`$ (7)
$`=`$ $`{\displaystyle \frac{d^2p_{}\mathrm{\Psi }^{}(\stackrel{}{p}_{}^{})\mathrm{\Psi }(\stackrel{}{p}_{})}{\sqrt{2E_\stackrel{}{p}2E_\stackrel{}{p}^{}}}𝑑x^{}e^{ix\overline{p}^+x^{}}p^{}|\overline{\psi }(\frac{x^{}}{2},\stackrel{}{0}_{})\psi (\frac{x^{}}{2},\stackrel{}{0}_{})|p}`$
$`=`$ $`{\displaystyle }{\displaystyle \frac{d^2p_{}\mathrm{\Psi }^{}(\stackrel{}{p}_{}^{})\mathrm{\Psi }(\stackrel{}{p}_{})}{\sqrt{2E_\stackrel{}{p}2E_\stackrel{}{p}^{}}}}H(x,\xi =0,t).`$
where $`\stackrel{}{p}_{}^{}=\stackrel{}{p}_{}+\stackrel{}{\mathrm{\Delta }}_{}`$ and $`p_z^{}=p_z`$, i.e. $`\xi =0`$.

The physics of Eq. (7) is the following: the non-trivial intrinsic structure of the target particle is expressed in the OFPD $`H(x,\xi ,t)`$. However, as the convolution of $`H(x,\xi ,t)`$ with the wave function reflects, the impact parameter dependent PD in the state $`\mathrm{\Psi }`$ gets spread out in position space due to the fact that the wave function $`\mathrm{\Psi }`$ does in general have a nonzero width in position space.

Intuitively, one would like to chose a wave packet that is point-like in position space (i.e. constant in momentum space) so that the $`\stackrel{}{b}_{}`$-dependence in Eq. (6) is only due to the intrinsic structure of the target particle but not due to the wave packet used to ‘nail it down’. However, one need to be careful in this step (and being able to properly address these issues was the sole reason for working with wave packets) because as soon as one localizes a particle to a region of space smaller than its Compton wavelength, its motion within the wave packet becomes relativistic and therefore the structure of the particle gets affected by Lorentz contraction as well as other relativistic effects.

### 2.1 Nonrelativistic limit

It is very instructive to consider the nonrelativistic (NR) limit first, where none of these complications occur. Formally, the simplification arises since $`E_\stackrel{}{p}=E_\stackrel{}{p}^{}=M`$ and therefore $`\mathrm{\Delta }^2=\stackrel{}{\mathrm{\Delta }}^2=\stackrel{}{\mathrm{\Delta }}_{}^2`$ First of all, this means that one can pull $`H(x,0,\stackrel{}{\mathrm{\Delta }}_{}^2)`$ out of the integral in Eq. (7), yielding

$$f_\mathrm{\Psi }(x,\stackrel{}{\mathrm{\Delta }}_{})=H(x,0,\stackrel{}{\mathrm{\Delta }}_{}^2)\frac{d^2p_{}\mathrm{\Psi }^{}(\stackrel{}{p}_{}^{})\mathrm{\Psi }(\stackrel{}{p}_{})}{2M}.$$
(8)
In order to proceed further, we choose a wave packet that is very localized in transverse position space. Specifically, we choose a packet whose width in transverse momentum space is much larger than a typical QCD scale. That way, the $`\stackrel{}{\mathrm{\Delta }}_{}`$-dependence on the r.h.s. of Eq. (8) is mostly due to the matrix element and not due to the wave packet $`\mathrm{\Psi }`$. Therefore, by making the wave packet very localized in position space one obtains

$$f_\mathrm{\Psi }(x,\stackrel{}{\mathrm{\Delta }}_{})=H(x,\stackrel{}{\mathrm{\Delta }}_{}^2).$$
(9)
The dependence on the detailed shape of the wave packet has disappeared once it is chosen localized enough and it is thus legitimate to identify the Fourier transformed (w.r.t. $`\stackrel{}{\mathrm{\Delta }}_{}`$) $`\xi =0`$ OFPD with the impact parameter dependence of the parton distribution in the target particle itself.

### 2.2 Relativistic corrections

In order to have a unique definition of impact parameter dependent PDFs, i.e. a definition which does not depend on the exact shape of the wave packet, one would like to make the wave packet as small as possible in position space and certainly much smaller than the spatial extension of the hadron in the $``$ direction — otherwise Eq. (6) is dominated by the shape of the wave packet and not by the intrinsic $`\stackrel{}{b}_{}`$-dependence. However, once the wave packet is smaller than about a Compton wavelength of the target then relativistic corrections can no longer be ignored and may even become more important as the corrections due to the spatial extension of the wave packet.

This problem is well known from form factors, where the identification of the form factor with the Fourier transform of a charge distribution in position space is uniquely possible only for momenta that are much smaller than the target mass. One can also rephrase this statement in position space by saying that the identification of the Fourier transform of the form factor with a charge distribution works only if one looks at scales larger than about one Compton wavelength of the target. Details below this scale may depend on the Lorentz frame and are therefore not physical.

For OFPDs the problem is very similar. If one works for example in the rest frame then one faces the same kind of corrections that also affect the form factors, in the sense that the $``$ resolution is limited to about one Compton wavelength. More details can be found in Ref.  and are omitted here since the natural frame to interpret parton distribution functions is the infinite momentum frame (IMF), and as we will discuss in the next section, these relativistic corrections play no important role in the IMF.

### 2.3 Infinite momentum frame

In the following we chose a wave packet as in Eq. (5) but with $`p_zM`$, where $`M`$ is the target mass. Then one can for example expand 222Note that, by choice of the wave packet, $`\stackrel{}{\mathrm{\Delta }}\stackrel{}{p}^{}\stackrel{}{p}`$ is always transverse i.e. $`\mathrm{\Delta }_z=0`$.

$`E_{\stackrel{}{p}+\stackrel{}{\mathrm{\Delta }}}|p_z|+{\displaystyle \frac{\stackrel{}{p}_{}^2+2\stackrel{}{p}_{}\stackrel{}{\mathrm{\Delta }}_{}+\stackrel{}{\mathrm{\Delta }}_{}^2}{2|p_z|}}`$ (10)
and therefore the energy transfer

$$E_{\stackrel{}{p}+\stackrel{}{\mathrm{\Delta }}}E_\stackrel{}{p}\frac{2\stackrel{}{p}_{}\stackrel{}{\mathrm{\Delta }}_{}+\stackrel{}{\mathrm{\Delta }}_{}^2}{2|p_z|}$$
(11)
vanishess as $`p_z\mathrm{}`$. A more detailed analysis works as follows: Let us denote the typical momentum scale in the wave packet by $`\mathrm{\Lambda }_\stackrel{}{p}_{}`$ and the $``$ resolution that we are interested in by $`\mathrm{\Lambda }_\stackrel{}{\mathrm{\Delta }}_{}`$ (normally, $`\mathrm{\Lambda }_\stackrel{}{\mathrm{\Delta }}_{}`$ will be a typical hadronic momentum scale, i.e. on the order of $`1GeV`$). Then of course what we need to satisfy is $`\mathrm{\Lambda }_\stackrel{}{p}_{}\mathrm{\Lambda }_\stackrel{}{\mathrm{\Delta }}_{}`$. At the same time we need $`|p_z|\mathrm{\Lambda }_\stackrel{}{p}_{},\mathrm{\Lambda }_\stackrel{}{\mathrm{\Delta }}_{}`$, but this is no problem if we go to the IMF.

As a direct consequence of Eq. (11), the time-like component of the momentum transfer can be ignored, i.e. one can approximate $`\mathrm{\Delta }^2\stackrel{}{\mathrm{\Delta }}^2=\stackrel{}{\mathrm{\Delta }}_{}^2`$ and therefore one can again pull $`H(x,0,\stackrel{}{\mathrm{\Delta }}_{}^2)`$ out of the integral in Eq. (7), and one finds

$$f_\mathrm{\Psi }(x,\stackrel{}{\mathrm{\Delta }}_{})=H(x,0,\stackrel{}{\mathrm{\Delta }}_{}^2)\frac{d^2p_{}\mathrm{\Psi }^{}(\stackrel{}{p}_{}^{})\mathrm{\Psi }(\stackrel{}{p}_{})}{\sqrt{2E_\stackrel{}{p}2E_\stackrel{}{p}^{}}}.$$
(12)
Finally making use of the fact that we chose a wave packet that is very localized, i.e. $`\mathrm{\Lambda }_\stackrel{}{p}_{}\mathrm{\Lambda }_\stackrel{}{\mathrm{\Delta }}_{}`$, we can approximate $`\mathrm{\Psi }(\stackrel{}{p}_{}^{})\mathrm{\Psi }(\stackrel{}{p}_{})`$ in Eq. (12) yielding

$$f_\mathrm{\Psi }(x,\stackrel{}{\mathrm{\Delta }}_{})=H(x,0,\stackrel{}{\mathrm{\Delta }}_{}^2).$$
(13)
As in the NR case, the dependence on the wave packet has disappeared and the identification has become unique.

We have thus accomplished to show that, OFPDs at $`\xi =0`$ have in the IMF the simple physical interpretation as Fourier transforms of impact parameter dependent parton distributions with respect to the impact parameter. In other words, OFPDs, in the limit of $`\xi 0`$, allow to simultaneously determine the longitudinal momentum fraction and transverse impact parameter of partons in the target hadron in the IMF.

## 3 The transverse center of momentum

In the previous chapter, when we introduced the notion of PDFs as a function of the impact parameter, we did not specify with respect to which point (or line) the impact parameter is defined. In NR scattering experiments, the impact parameter usually refers to the center of mass, but it is perhaps not a priori obvious how to generalize the concept of a $``$ center of mass to the IMF.

In the infinite momentum or light-front (LF) frame there exists a residual Galilei invariance under the purely kinematic $``$ boosts

$`x_i`$ $``$ $`x_i^{}=x_i`$
$`\stackrel{}{k}_i`$ $``$ $`\stackrel{}{k}_i^{}\stackrel{}{k}_i+x_i\mathrm{\Delta }\stackrel{}{P}_{},`$ (14)
where we denote the longitudinal momentum fraction and $``$ momentum of the $`i`$-th parton in a given Fock component by $`x_i`$ and $`\stackrel{}{k}_i`$ respectively, the LF-Hamiltonian transforms covariantly, i.e. $`P^{}\frac{\stackrel{}{P}_{}^2}{2P^+}`$ remains constant, which resembles very much NR boosts

$$\stackrel{}{k}_i\stackrel{}{k}_i^{}=\stackrel{}{k}_i+m_i\mathrm{\Delta }\stackrel{}{v}=\stackrel{}{k}_i+\frac{m_i}{M}\mathrm{\Delta }\stackrel{}{P},$$
(15)
with $`E\frac{\stackrel{}{P}^2}{2M}`$ remaining constant. Because of this similarity between NR boosts and $``$ boosts in the IMF, many familiar results about NR boosts can be immediately transferred to the IMF.

For this work, the most important result concerns the physical interpretation of the form factor as the Fourier transform of the charge distribution in a frame where the center of mass

$$\stackrel{}{R}_{CM}\underset{i}{}\frac{m_i}{M}\stackrel{}{r}_i$$
(16)
is at the origin. One can also rephrase this result by stating that the Fourier transform of the form factor is the distribution of the charge as a function of the distance from $`\stackrel{}{R}_{CM}`$.

By comparing Eqs. (14) and (15) it becomes clear that in the IMF the momentum fractions $`x_i`$ play the role of the mass fraction $`\frac{m_i}{M}`$. Hence it is not very surprising that the IMF analog to the NR center of mass $`\stackrel{}{R}_{CM}`$ is given by the transverse center of momentum

$$\stackrel{}{R}_{}\underset{i}{}x_i\stackrel{}{r}_i,$$
(17)
i.e. a weighted average of $``$ positions, but where the weight factors are given by the (light-cone) momentum fractions and not the mass fractions 333Note that for NR systems the momentum fractions are given by the mass fractions..

The Fourier transform of the $`\xi =0`$ (i.e. ‘unskewed’) OFPD

$$F(x,\stackrel{}{b}_{})d^2\mathrm{\Delta }_{}e^{i\stackrel{}{b}_{}\stackrel{}{\mathrm{\Delta }}_{}}H(x,0,\stackrel{}{\mathrm{\Delta }}_{}^2)$$
(18)
can thus be interpreted as the (light-cone) momentum distribution of partons as a function of the $``$ separation from the $``$ center of energy-momentum $`\stackrel{}{R}_{}=_ix_i\stackrel{}{r}_i`$ of the target. 444The above results can be verified by starting with $`\xi =0`$ OFPDs calculated from a LF Fock space expansion for the hadron state .

The above observation about the $``$ center of momentum has one immediate consequence for the $`x1`$ behavior of $`F(x,\stackrel{}{b}_{})`$. Since the weight factors in the definition of $`\stackrel{}{R}_{}`$ are the momentum fractions, any parton $`i`$ that carries a large fraction $`x_i`$ of the target’s momentum will necessarily have a $``$ position $`\stackrel{}{r}_i`$ that is close to $`\stackrel{}{R}_{}`$. Therefore the transverse profile (i.e. the dependence on $`\stackrel{}{b}_{}`$) of $`F(x,\stackrel{}{b}_{})`$ will necessarily become more narrow as $`x1`$, i.e. we expect that partons become very localized in $``$ position as $`x1`$. By Fourier transform, this also implies that the slope of $`H(x,0,t)`$ w.r.t. $`t`$ at $`t=0`$, i.e.

$$\stackrel{}{b}_{}^24\frac{\frac{d}{dt}H(x,0,t)|_{t=0}}{H(x,0,0)}$$
(19)
should in fact vanish for $`x1`$!

## 4 Explicit example (model calculation)

In order to illustrate the results from the previous sections in a concrete example, we make a Gaussian ansatz for the $`\stackrel{}{k}_i`$ dependence of the light-cone wave function in each Fock component

$$\mathrm{\Psi }(x_i,\stackrel{}{k}_{})\mathrm{exp}\left[a^2\underset{i=1}{\overset{N}{}}\frac{\left(\stackrel{}{k}_ix_i\stackrel{}{P}_{}\right)^2}{x_i}\right]\mathrm{\Psi }(x_i),$$
(20)
where $`\mathrm{\Psi }(x_i)`$ is not further specified. The above model ansatz allows to carry out the $`\stackrel{}{k}_{}`$ momentum integrals explicitly, yielding

$$H(x,0,t)=\mathrm{exp}\left[\frac{a^2}{2}\frac{1x}{x}t\right]f(x),$$
(21)
i.e.

$$F(x,\stackrel{}{b}_{})\frac{x}{1x}\mathrm{exp}\left[\frac{1}{2a^2}\frac{x}{1x}\stackrel{}{b}_{}^2\right]f(x).$$
(22)
Of course, the above model (20) clearly oversimplifies the actual complexity of the nucleon’s Fock space wavefunction, but the qualitative behavior of the final result (22) seems reasonable: as predicted in the previous section, the width of $`F(x,\stackrel{}{b}_{})`$ in $``$ position space goes to zero as $`x1`$. For decreasing $`x`$, the distribution widens monotonically until it diverges for $`x0`$. Of course, the divergence of the $``$ width as $`x0`$ is probably an artifact of the model ansatz which seems is more suitable for intermediate to large values of $`x`$. Nevertheless, it seems reasonable that the $``$ width increases with decreasing $`x`$ since for example the pion cloud, which is expected to contribute only at smaller values of $`x`$, should me more spread out than the valence partons, which are expected to dominate at larger values of $`x`$.

Notice also that both the above model calculation as well as the pion cloud picture suggest that the $`\stackrel{}{b}_{}`$-dependence in $`F(x,\stackrel{}{b}_{})`$ does not factorize, which of course also translates back into a lack of factorization of the $`t`$ dependence in $`H(x,0,t)`$. A similar lack of factorization was also observed in Ref. .

## 5 Practical aspects

From the experimental point of view, $`\xi =0`$ is not directly accessible in DVCS since one needs some longitudinal momentum transfer in order to convert a virtual photon into a real photon. There are several ways around this difficulty. Of course, one way could be to access $`\xi =0`$ in real wide angle Compton scattering . However, it should also be possible to perform DVCS experiments at finite $`\xi `$ and to extrapolate to $`\xi =0`$.

Extrapolation to $`\xi =0`$ is greatly facilitated by working with moments since the $`\xi `$ dependence of the moments of OFPDs is in the form of polynomials . For example, for the even moments of $`H(x,\xi ,t)`$ one finds

$`H_n(\xi ,t)`$ $``$ $`{\displaystyle _1^1}𝑑xx^{n1}H(x,\xi ,t)={\displaystyle \underset{i=0}{\overset{\left[\frac{n1}{2}\right]}{}}}A_{n,2i}(t)\xi ^{2i}+C_n(t)`$
$`=`$ $`A_{n,0}(t)+A_{n,2}(t)\xi ^2+\mathrm{}+A_{n,n2}(t)\xi ^{n2}+C_n(t)\xi ^n,`$
i.e. for example

$$_1^1𝑑xxH(x,\xi ,t)=A_{2,0}(t)+C_2(t)\xi ^2.$$
(24)
Since the $`H_n`$ have such a simple functional dependence on $`\xi `$, one can thus use measurements of the moments of OFPDs at nonzero values of $`\xi `$ to determine (fit) the “form factors” $`A_{n,2i}(t)`$ and $`C(t)`$. Once one has determined these invariant form factors, one can go back and evaluate Eq. (5) for $`\xi =0`$, yielding

$$H_n(\xi =0,t)=A_{n,0}(t),$$
(25)
which then allows one to determine the impact parameter dependence of the $`nth`$ moment of the parton distribution in the target, using

$$q_n(\stackrel{}{b}_{})_1^1𝑑xx^{n1}q(x,\stackrel{}{b}_{})=d^2q_{}A_{n,0}(\stackrel{}{\mathrm{\Delta }}_{}^2)e^{i\stackrel{}{\mathrm{\Delta }}_{}\stackrel{}{b}_{}}.$$
(26)
A very similar procedure can be applied to the spin dependent OFPD $`\stackrel{~}{H}(x,\xi ,t)`$, whose even moments of $`\stackrel{~}{H}(x,\xi ,t)`$ satisfy

$`\stackrel{~}{H}_n(\xi ,t)`$ $``$ $`{\displaystyle _1^1}𝑑xx^{n1}\stackrel{~}{H}(x,\xi ,t)={\displaystyle \underset{i=0}{\overset{\left[\frac{n1}{2}\right]}{}}}\stackrel{~}{A}_{n,2i}(t)\xi ^{2i}`$
$`=`$ $`\stackrel{~}{A}_{n,0}(t)+\stackrel{~}{A}_{n,2}(t)\xi ^2+\mathrm{}+\stackrel{~}{A}_{n,n1}(t)\xi ^{n1}.`$
Similar to the unpolarized case, the $`\frac{n+1}{2}`$ form factors of the $`n^{th}`$ moment can be obtained from measurements of $`\stackrel{~}{H}`$ for $`\frac{n+1}{2}`$ different values of $`\xi `$ (and the same values of $`t`$) and the impact parameter dependence of the $`n^{th}`$ moment of the polarized parton distribution reads

$$\mathrm{\Delta }q_n(\stackrel{}{b}_{})_1^1𝑑xx^{n1}\mathrm{\Delta }q(x,\stackrel{}{b}_{})=d^2q_{}\stackrel{~}{A}_{n,0}(\stackrel{}{\mathrm{\Delta }}_{}^2)e^{i\stackrel{}{\mathrm{\Delta }}_{}\stackrel{}{b}_{}}.$$
(28)
Of course, this procedure becomes rather involved for high moments, but the steps outlined above seem to be a viable way of determining the impact parameter dependence of low moments of parton distributions from DVCS data.

## 6 Summary and outlook

Off-forward parton distributions for $`\xi 0`$, i.e. where the off-forwardness is only in the $``$ direction, can be identified with the Fourier transform of impact parameter dependent parton distributions w.r.t. the impact parameter $`b`$, i.e. the $``$ distance from the center of (longitudinal) momentum in the IMF. This identification is very much analogous to the identification of the charge form factor with the Fourier transform of a charge distribution in position space.

The $`\xi 0`$ limit of OFPDs is difficult to access in DVCS. However, as we illustrated in Sec.5, one can also measure $`x`$-moments for nonzero values of $`\xi `$ and use those to construct moments for $`\xi =0`$.

Knowing the impact parameter dependence allows one to gain information on the spatial distribution of partons inside hadrons and to obtain new insights about the nonperturbative intrinsic structure of hadrons. For example, the pion cloud of the nucleon is expected to contribute more for large values of $`b`$. Shadowing of small $`x`$ parton distributions, is probably stronger at small values of $`b`$ since partons in the geometric center of the nucleon are more effectively shielded by the surrounding partons that partons far away from the center. These and many other models and intuitive pictures for the parton structure of hadrons give rise to predictions for the impact parameter dependence of PDFs that reflect the underlying microscopic dynamics of these models. Through the experimental measurement of OFPDs using DVCS, it will for the first time become possible to obtain experimental information on the impact parameter dependence of parton structure which will thus provide much more comprehensive tests on our understanding of nonperturbative parton structure.

## Acknowledgments

This work was supported by a grant from the DOE (DE-FG03-95ER40965) and in part by TJNAF.

## References


format: text

text:
# A Measurement of the Average Longitudinal Development Profile of Cosmic Ray Air Showers Between 𝟏𝟎^𝟏𝟕eV and 𝟏𝟎^𝟏𝟖eV

## I Introduction

The atmospheric nitrogen fluorescence light technique plays an increasingly important role in extremely high energy cosmic ray observation. The success of the Fly’s Eye and HiReshires  experiments encourages the employment of this technique in future projects like the Auger experimentAuger , the Telescope Array proposalTA  and space based experiments like EUSO EUSO  and OWLOWL . In a recent report on the change of cosmic ray composition in the vicinity of 5$`\times `$1017 eV from HiRes/MIA joint dataprl , the authors have pointed out that the reconstruction of the shower longitudinal development depends upon the assumed functional form of the features of the shower development. How much do the characteristic parameters of the shower longitudinal development such as the shower maximum and its location, shower rise and shower decay constant depend on which function is employed? In the new technique of shower reconstruction developed for the HiRes monocular data HR1 , the shower development function plays an even more fundamental role in the shower geometry reconstruction. In effect, the shower geometry is varied until the correct shower width predicted by the average shower function is obtained. The use of an appropriate shower development function is crucial in the determination of shower energy because this is based on the integral of this function. In practice, empirical functions based on data at lower energies or based on theoretical electromagnetic cascade calculation are used at the highest energies, e.g. above 1017eV. However, none of these has been experimentally tested at these energies in the atmosphere. Values of parameters used in trial functions such as the $`\lambda `$ in the Gaisser-Hillas functiongh  and $`L_0`$ (for a hadronic shower) in the Greisen function Greisen  have never been measured in experiments at such a high energy. This is sufficient motivation to carefully test these well known functions and determine those parameters directly from experimental data. The HiRes/MIA joint experiment meets the necessary conditions for such a high resolution measurement of the shower longitudinal development.

The HiRes prototype detector provides a degree by degree measurement of the intensity of the fluorescence light produced by the shower electrons along the shower axis. This light intensity, after corrections, is proportional to the shower development profile. The angular bin signals, however, suffer fluctuations and are contaminated by the direct and scattered Cerenkov light component produced by the same electrons. The source of the fluctuation is sky noise and the statistical fluctuation due to finite sampling of the light. A method to minimize the bin by bin fluctuation is to take averages of the angular bin signals over all the events. This will establish an almost perfect “average shower longitudinal development profile” which can be used to test the various trial functions. In order to do so, however, all showers must be “aligned” and “normalized” so that the showers can be compared with each other in spite of variations of position of shower maximum and shower energy. To carry out the alignment and normalization properly requires that the shower geometry must be well determined. In addition, the Cerenkov components of the signals must be subtracted. Fluorescence light is proportional to the shower size and the proportionality coefficient is geometry related. In the HiRes/MIA joint experiment, the shower geometry is well determined with the help of muon front timing information from MIA, the wide elevation angle range (up to 70) of the HiRes prototype detector allows for a broad range of shower development to be detected, and the short shower-detector distance (about 3 km) minimizes uncertainties related to the atmospheric attenuation and scattering. We use the same data set selected for the study of the composition of primary cosmic raysprl . The criteria used to select this sample are well justified there based on a full detector simulation apj .

In the following sections, we will describe the experiment and data set briefly, explain how to minimize the contribution of Cerenkov light and how to subtract the Cerenkov light component in each degree bin, propose a way to normalize the individual showers, and present the average profile and compare it with well-known functions.

## II HiRes/MIA Joint Experiment and Data Set

The HiRes prototype is situated at $`112^{}\text{W}\text{longitude}`$ and $`40^{}\text{N}\text{latitude}`$ and at a vertical atmospheric depth of $`850`$ g/cm2. It ishires  composed of 14 optical reflecting telescopes. They image the extensive air shower(EAS) as it progresses through the field of view of the telescope. Nitrogen fluorescence light (300–400 nm) is emitted at an atmospheric depth $`X`$ in proportion to the number of charged particles in the EAS at that depth. The electrons radiate Cerenkov light in the same band simultaneously. Because of the very foreward angular distribution of the Cerenkov light, only those events pointing toward the detector will give rise to signals dominated by Cerenkov light. However, the Cerenkov light can be scattered by molecular and aerosol particles into almost all pixels along the shower track no matter what the shower geometry is. Therefore, the signal recorded by the HiRes detector consists of the sum of fluorescence light, Cerenkov light and various kinds of noise.

The MIA detector is situated about 3.3 km from HiRes and 150 m lower in elevation. It casamia  consists of 2500 m2 of active area distributed in 16 patches of 64 scintillation counters and measures the EAS muon arrival times with a precision of 4 ns and records all hits occurring within 4 $`\mu \text{s}`$ of the system trigger. Combining such well determined timing from MIA with the triggered HiRes pixel directions and HiRes timing information allows precise reconstruction of the shower geometry. This is essential in the determination of the atmospheric depth at which the shower passes the triggered HiRes pixels and the optical path length between the shower and the detector.

HiRes/MIA coincident data were collected on clear moonless nights between Aug. 23, 1993 and Aug. 24, 1996. The total coincident exposure time was 2878 hours corresponding to a duty cycle of 10.9%. 4034 coincident events were observed. For events passing a set of coincidence assurance cuts the shower trajectory, including arrival direction and core location for each event was obtained in an iterative procedure using the information from both HiRes and MIA Brian . The accuracy of the shower axis determination depends on the number of observed muons, the HiRes angular track length and the core distances from MIA and HiRes. 2491 events are reconstructed via this procedure. Monte Carlo(MC) studies apj  show that the median shower direction error is $`0.85^{}`$ with a median core location error of 45 m.

To insure data quality and maintain good resolution we require that for each event the depth of shower maximum $`X_m`$ is visible within a minimum observed slant depth interval of 250 g/cm2, that the track subtends at least 20$`^{}`$, that the accumulated gap between the fired pixels is less than 40% of the total gramage spanned by the shower, that the uncertainty in $`X_m`$, $`\mathrm{\Delta }X_m`$, is less than 50 g/cm2, that the reduced $`\chi ^2`$ for the profile fit not exceed 10 and that the MIA to core distance $`R_{p_{MIA}}`$ is less than 2000 m. Additionally, we require the minimum viewing angle of the fired tube, $`\theta _h`$, to be greater than 20 in order to delete events which possess large direct Cerenkov light components. These cuts, as summarized in TABLE 1, leave a sample of 488 events. The energy distribution peaks at 3$`\times `$1017 and the average Cerenkov light fraction is 15%. There is no event in this sample in which direct Cerenkov light contributes more than 1% of the total signal.

We note that the shower development profile has been fit to the Gaisser-Hillas(G-H) function in order to carry out some of the cuts. However, the use of this function as well as the fit parameters derived from it is for the purpose of event selection only. Hereafter, only the raw data for these selected events is used in the analysis.

## III Cerenkov Light Subtraction

HiRes pixel signals are rearranged into a one dimensional series of one degree bins along the shower track by considering the light acceptance of the mirrors according to the geometric position of the shower. As mentioned above, these HiRes bin signals include both fluorescence and Cerenkov light, although the total fraction of direct Cerenkov light is kept low by the imposition of the cuts described above. Bins which correspond to the deeper positions in the atmosphere tend to have more scattered Cerenkov light because the Cerenkov beam accumulates with increasing shower development and the molecular and aerosol scattering gets stronger as the shower penetrates deeper into the atmosphere. The raw signal versus depth, therefore, is not in itself a representation of the shower longitudinal development. We must subtract the scattered Cerenkov light from the observed signal. The resultant corrected signal composed only of the fluorescence light contribution is proportional to the shower size at the corresponding atmospheric depth except for the natural and electronic noise component.

The Cerenkov light produced in an angular bin is calculated by multiplying the average Cerenkov light yield of a single electron with the shower size at the corresponding depth $`X`$. The average yield is taken over the energy of shower electrons above the Cerenkov threshold, $`E_t`$, and over the emitting angle. The energy spectrum of shower electrons, the angular distribution of Cerenkov light and the pressure dependence of $`E_t`$ and yield are carefully considered. Details can be found in the Ref. FEnim  and its references.

The contribution of direct Cerenkov light is minimized by the data cuts. The contribution to the bin signal from scattered light from the accumulated Cerenkov beam is estimated as follows: we assume that the Cerenkov light beam in the angular bin is due to shower electrons in the previous adjacent angular bin. After taking into account the attenuation of the Cerenkov light beam traveling through the air between the bins, the contribution of scattered Cerenkov light from the beam into the detector is estimated based on the Rayleigh scattering theory and Mie scattering by aerosol.

The observed signal in this angular bin can then be modified by subtracting the scattered Cerenkov light. The corrected signal can be converted into a shower size for the sake of computing the contribution to the Cerenkov light beam in the next adjacent angular bin. This recurrent procedure requires knowledge of the Cerenkov light beam component in the first angular bin in the field of the view of the detector. If the first signal bin were at the beginning of the shower development, the Cerenkov light component would be zero. The observed signal in this bin would be pure fluorescence which can be converted directly into a shower size. We assume that the first detected bin corresponds to an early enough stage in the shower development so that the zero-Cerenkov-light-component assumption is still approximately true. This subtraction scheme systematically underestimate the Cerenkov fraction in each subsequent bin because the first detected signal from the shower is not at the very beginning of the shower cascade. However, we will find that this systematic error is negligible and the proposed Cerenkov light subtraction is sufficient for our purposes.

## IV The Normalization of Air Shower Development Profile

As mentioned above, all the showers must be aligned in atmospheric depth and be normalized in shower size before one can take the average in bins over events. Most of our events follow a transition curve, namely the shower size increases rapidly as it develops in the atmosphere until reaches its maximum, $`N_m`$, then it starts decreasing because shower electrons begin to lose more energy by ionization than by radiation of high energy gammas. The maximum size of a shower is approximately proportional to the energy of the primary particle. After subtracting the Cerenkov light components, the signals are converted into shower sizes by taking the geometry and atmosphere related light collecting efficiency and attenuation effect into account. The showers are normalized with respect to energy by their maximum size, i.e. $`N(X)/N_m`$, denoted as $`n(X)`$. All the shower sizes are normalized to one as they reach their maxima.

Air shower development fluctuates in atmospheric depth due to fluctuations in the hadronic multiparticle production (which depends on the nature of the primary nucleus) and electromagnetic processes. As a consequence, the position of the shower maximum, $`X_m`$, varies from event to event. One can express the shower longitudinal development as a function of the shower “age”, defined as $`s=\frac{3X}{X+2X_m}`$, instead of depth, $`X`$. By using $`s`$, the shower longitudinal development is universally described as a rising phase from 0 (initial position of the shower) to 1 (shower maximum) and the decay phase from 1 to 3 (infinite depth). The physical shower has an effective extent from 0 to 2, however.

After the Cerenkov light subtraction, the observed signal profiles of individual events span a range of atmospheric depth and show the shower maximum through their convex-upwards shape. The best way to find the position of shower maximum in the presence of statistical fluctuations is to carry out a local fit with a rather general function. We use a parabolic function to fit the data in order to determine both the location and the size of the shower maximum. The only constraint to this function is that it must be convex upwards. Here, “local fitting” means that the fitting procedure is only applied to those data points near the maximum to minimize the bias in the location of the maximum caused by using such a symmetric trial function. The method can be iterated to remove points far from the maximum. We find that the iteration affects the location of shower maximum by less than 4% however and we do not use it in what follows.

With shower maximum and size at maximum determined, the showers are normalized in size and aligned in depth. FIG 1 shows all 488 showers plotted with normalized shower size, $`n`$, versus shower age, $`s`$. In order to more easily see the change of the density of the dots in the scatter plot, we let each entry represents three measurements. 6069 entries in total out of 18230 bins signals are picked randomly and plotted. We note that there is a 20% fluctuation in $`n`$ near the shower maximum. The distribution of $`n`$ at shower maximum is shown in FIG 2. The residuals in $`n`$, defined as the ratio between the deviation from the expectation value and the error associated with the individual bin signal, are well distributed as a Gaussian centered at 0 with a width of 0.94. The main sources of this fluctuation are sky noise, statistical fluctuation, the error in shower geometry determination, the error in Cerenkov light subtraction and the error in determining shower maximum. The first three, and particularly the first one, contribute more than 80% to the total fluctuation in the signal. The error from Cerenkov subtraction contributes the least.

The data is binned into age intervals, instead of the original angular bin. We take an average of the normalized shower size within each bin over all events to obtain an average shower development profile. The curve is well determined between $`s`$=0.5 and 1.25. We demonstrate the effect of the Cerenkov light subtraction in FIG 3, by plotting the average shower transition curves before and after the subtraction in the same plot.

## V Tests of trial functions

After calculating the average normalized shower longitudinal development profile, we fit the shape of the profile with several well known trial functions such as the Gaisser-Hillasgh , GreisenGreisen  and Gaussian forms.

In order to describe the characteristics of the shower longitudinal development, i.e. the asymmetric rising-falling shape, Gaisser and Hillas introduced a function which depends on the essential characteristics an air shower: the initial point, $`X_0`$, the shower maximum, $`N_m`$, the shower maximum location, $`X_m`$, and the shower decay length, $`\lambda `$. Except for $`N_m`$ which is in number of shower electrons, all the parameters are in units of atmosphere depth, $`g/cm^2`$. The function describes the change of the shower size with the atmospheric depth as

$`N(X)=N_m\left({\displaystyle \frac{XX_0}{X_mX_0}}\right)^{\frac{X_mX_0}{\lambda }}e^{\frac{X_mX}{\lambda }},`$ (1)
where the atmospheric depth, $`X`$, is in $`g/cm^2`$.

Translating the depth $`X`$ into age $`s`$ and using the normalized shower size $`n=N/N_m`$, eqn. (1) becomes

$`n(s)=\left(1{\displaystyle \frac{1s}{3s}}{\displaystyle \frac{3T_m}{T_mT_0}}\right)^{T_mT_0}e^{3T_m\frac{1s}{3s}},`$ (2)
where $`T_m=X_m/\lambda `$ and $`T_0=X_0/\lambda `$ are the two remaining parameters. Parameter $`T_0`$ is constrained to be less than $`\frac{2s_{min}}{3s_{min}}T_m`$, where $`s_{min}`$ is the lower limit of the data points (about 0.5 as shown in FIG. 3). The comparison of this function with the data is shown in FIG. 4. The two parameters are strongly correlated.

The fit is poor at small and large $`s`$, particularly beyond $`s=1.3`$. While a possible reason for the deviation of those 3$``$4 points may be poor statistics, a detailed systematic error analysis is described in the next section.

The second trial function is the Greisen function which was suggestedGreisen  to describe the development of a pure electromagnetic shower. Depth is therefore expressed in the unit of radiation length, $`L_0=36.66g/cm^2`$. The function is defined by a single parameter, $`y=X_m/L_0`$, and has a form

$`N(T)={\displaystyle \frac{0.31}{\sqrt{y}}}e^{T(1\frac{3}{2}\mathrm{𝑙𝑛}s)},`$ (3)
where $`T=X/L_0`$ refers to the atmospheric depth in radiation lengths in the air. Since the age $`s`$ appears explicitly in the formula, it is straightforward to rewrite it as a function of either $`s`$ or $`T`$. It is apparent after conversion to a function of $`s`$ that the radiation length no longer appears explicitly in the formula. This feature of the Greisen function make it potentially useful for the description of hadronic shower with a single parameter $`y`$. However, the meaning of $`L_0`$ need to be changed. In the case of electromagnetic cascade, Greisen shows that $`y=\mathrm{𝑙𝑛}(E_0/E_c)`$ fits a elongation law with a rate about 60g/cm2 where the $`E_c`$ refers to the electron critical energy. We will see how the rate is changed to fit the case of hadronic shower below. The fit of the Greisen function to our measured average shower is shown in the dashed line in FIG. 4.

In the Fly’s Eye data analysisFEcomp , the poorer detector resolution dominated over the shape of the shower longitudinal development and a symmetric Gaussian in $`X`$ could be used to represent the shower longitudinal development. Using $`\xi =\sigma /X_m`$ as the single parameter, the Gaussian can be normalized by following the rule set in previous sections, i.e.

$`n(s)=exp\left\{{\displaystyle \frac{1}{2\xi ^2}}\left({\displaystyle \frac{1s}{1s/3}}\right)^2\right\}.`$ (4)
The result is also plotted in FIG 4. It is a poor fit to the data with a $`\chi ^2=26`$ per degree of freedom. Thanks to the significantly improved resolution, this experiment provides the first clear observation of an asymmetrical air shower development profile in $`X`$ at 10$`{}_{}{}^{17}`$ 1018eV.

However, based on the result of this experiment, we observe that the average longitudinal development profile of air showers appears quite symmetrical as a function of shower age, $`s`$. A simple symmetrical Gaussian shaped function of age can be used to describe shower profiles. Since the function is centered at 1, it is governed by only one parameter, i.e. the width of the function, denoted as $`\sigma `$. This fit is shown in FIG. 4 also. Hereafter, we call this function the “new Gaussian function” to distinguish it from the one that is symmetric in the depth variable.

We plot the fit results in the form of deviations in FIG. 5. Quantitatively, the $`\chi ^2`$ per degree of freedom for Gaisser-Hillas, Greisen and the “new Gaussian function” are 1.93, 1.87 and 1.79 respectively. They all have essentially the same goodness of fit. More data beyond the region currently covered is required to further differentiate between them. From a practical standpoint however, the Gaisser-Hillas function contains too many parameters, especially $`X_0`$ which is not directly measured in experiments and is consequently poorly defined. These parameters are strongly correlated. The Greisen function does not have this problem because it fits the data with a single parameter. Both functions need a scale to measure the shower development, however. The scale is given by $`\lambda `$ and $`L_0`$, respectively. The scale must be associated with the type of primary particle. This actually reduce the freedom of fitting in the reconstruction of real experimental events. One has to set the scale as a parameter in a complete fitting procedure. The “new function” proposed in this work is more practically useful in real shower reconstruction. As a function of depth, it has the form

$`n(X)=exp\left\{{\displaystyle \frac{2}{\sigma ^2}}\left({\displaystyle \frac{XX_m}{X+2X_m}}\right)^2\right\}.`$ (5)
This function uses a single parameter, $`X_m`$, to locate a shower in the atmosphere, and a single parameter, $`\sigma `$, to indicate the shower width and it is scale-free. The last feature is very useful since the composition is a-priori unknown.

## VI Systematic error analysis

In order to investigate the reasons for deviation from the test functions at both early and late stages of shower development, we apply the above analysis to Monte Carlo simulated events. The generator used here is the full simulation code of the HiRes/MIA detector developed in the composition studyprl . It is driven by a shower generator based on the CORSIKA package. It produces the shower size using the Gaisser-Hillas function with full fluctuations in $`X_m`$, $`N_m`$, $`X_0`$ and $`\lambda `$. The fluctuation in muon density and muon arrival time on the ground is generated according to parameterizations of CORSIKA simulated events. All details of the detector are carefully considered including the sky noise, atmospheric attenuation, “ray-tracing”, phototube response, electronics and triggering. Generated data are passed through the complete calibration and reconstruction procedure just like the real data. The details of the generator and the comparison with the data can be found in the Ref. apj . FIG. 6 shows that the overall quality of the fit is good. This means that the Cerenkov subtracting scheme, shower maximum finding and normalization procedure proposed in this paper work well. The shower longitudinal development function is reconstructed almost exactly except for the few points at the beginning and the end of the shower as is the case for the real data. The deviation turns out to be associated with the correction to tube signals which lie far from the shower-detector-plane. Details of this correction can be found in Ref. Brian . Briefly, a tube can be triggered even if it is not located exactly on the shower-detector plane since the shower has a finite lateral width, the mirror image on the focal plane has a finite spot size and the broadening effect due to the atmospheric scattering. The amount of light reaching the detector can be calculated by correcting the amount of measured light with a response function in the off-plane-angle. This function is determined by performing a careful “ray-tracing” procedure which folds all of the effects in. It turns out to be a Gaussian with a typical width of a few tenths of a degree depending on the shower distance and the shower age. On the other hand, we know that there is an uncertainty in the determination of the shower-detector plane direction of about 0.7 degreesapj . This uncertainty causes a systematic overcorrection in the signals for off plane tubes. This becomes worse for a tube located far from the center of the phototube cluster. This effect can be verified by simply substituting the ideal MC input shower-detector plane for calculating the correction and comparing with the correction using the fitted shower-detector plane. This is shown in FIG. 6. The effect on the beginning and end of the average profile is clear.

One way to handle this systematic error would be to correct it based on this MC simulation result. Another approach is to simply drop those affected points in the fitting procedures. We prefer the later for several reasons. First, it avoids introducing any modeling dependence. Secondly, the points being dropped at age $`>1.25`$, correspond to positions deeper than 1.43$`X_m`$, or about 860$`g/cm^2`$ in the atmosphere, where the uncertainty in atmospheric attenuation is greatest. As a part of the systematic uncertainty study, we varied the atmospheric parametersprl  by one standard deviation from their most probable values in the shower reconstruction. We observed that there is no effect for most of the average profile and that this variation affected only the last few deepest points. This effect is smaller than that produced by the shower-detector plane error. Finally, in practice, dropping those points does not change our result because of the relatively large statistical errors associated with those points as can be seen in comparing the FIG. 4 with FIG. 5, where events associated with the last three points in FIG. 4 are dropped.

We now discuss the systematic error caused by setting the Cerenkov light to be zero in the first angular bin in the field of view of the detector. Based on the previous results, we known that the Gaisser-Hillas function works well in describing the shower longitudinal development. We can use the Gaisser-Hillas function in a least $`\chi `$-square fitting procedure event by event. Since in this case we know the shape of the shower, the Cerenkov light component in all angular bins can be calculated. The fraction of Cerenkov light in each bin is compared to our original method in FIG. 7. It is clear that most of the Cerenkov light is subtracted with our bin signal based Cerenkov light subtracting scheme. Only 7.7% of the Cerenkov light is systematically underestimated. This is a minor effect in calculating the average transition curves as shown in the FIG. 3. The diamonds in the figure are obtained by carrying out the same averaging and normalization procedure as before but with the Cerenkov light fraction in each bin replaced by that determined using the Gaisser-Hillas function. The points are visibly lower which makes the profile appear slightly wider. However, the widening is not large enough to change our conclusion on the shape of the longitudinal development function. For instance, the values of fitted parameter $`T_0`$ and $`T_m`$ are -2.57 and 8.6, respectively, which are consistent with the result listed below.

The other systematic issue is the reliability of our maximum searching by using local fitting with a parabolic function. The uncertainty associated with this method is demonstrated by comparing the maxima with those obtained from the reconstruction procedure using the Gaisser-Hillas profile. FIG. 8 a) shows the comparison between the $`X_m`$’s and FIG. 8 b) shows the difference between the $`N_m`$’s. The other way to illustrate the uncertainty is to compare the resolution functions in $`X_m`$ and $`N_m`$ from simulation for both methods. It turns out the resolution in $`X_m`$ for local parabolic fitting is 16% worse than the standard Gaisser-Hillas approach, and the resolution in $`N_m`$ is comparable in both cases.

## VII Results and conclusions

In order to extract the best values for the parameters, $`T_m`$, $`T_0`$, $`y`$ and $`\sigma `$, in the Gaisser-Hillas, Greisen and the “new Gaussian function”, a well reconstructed data set is used. The parameters are determined to be

$`T_m`$ $`=`$ $`9.7\pm 2.0,`$ (6)
$`T_0`$ $`=`$ $`3.2\pm 2.9,`$ (7)
$`y`$ $`=`$ $`11.11\pm 0.34,`$ (8)
$`\sigma `$ $`=`$ $`0.272\pm 0.002.`$ (9)
A strong correlation between $`T_m`$ and $`T_0`$ is found for the Gaisser-Hillas function. This correlation can be expressed as $`T_00.10T_m^2+0.74T_m0.10`$ and is shown in FIG. 9.

These results imply: 1) the strong correlation between parameters in the Gaisser-Hillas function means that the data can be fit with a simpler form. Alternatively one can fix one of those parameters and apply the correlation between them in the shower reconstruction; 2) The parameter $`\lambda `$ in the Gaisser-Hillas function can be estimated to be $`67\pm 11g/cm^2`$ by using the average $`X_{max}=651.9\pm 3.5g/cm^2`$; 3) the parameter $`T_0`$, or $`X_0`$, is likely to be a negative number in any individual shower fit. It is thus difficult to interpret as the first interaction point of a shower. 4) the parameter $`y`$, defined as $`X_m/L_0`$ in the Greisen function, is close to $`T_m`$, defined as $`X_m/\lambda `$. This implies that the best value of the parameter $`L_0`$, $`58.7\pm 2.0g/cm^2`$, is no longer the radiation length. This is to be expected since it is used to describe a hadronic shower; 5) a symmetrical function in age is a suitable representation of the shower longitudinal profile. This function is depth-scale invariant.

In summary, the shape of the extensive air shower longitudinal development is investigated in the energy range from 1017 to 1018eV with the HiRes/MIA hybrid experiment. The profile is quite symmetrical as a function of the age of the shower. The “new Gaussian function” (5), the Greisen function (3) and the Gaisser-Hillas function (2) describe the shower shape almost equally well, with the $`\chi ^2`$’s of the test 1.79, 1.87 and 1.93 respectively. This is the first direct measurement of the shower average longitudinal development profile at these energies and covering such a wide range of shower age.

## VIII Acknowledgements

We acknowledge the assistance of the command and staff of Dugway Proving Ground. We are greatly grateful for the discussion with James W. Cronin and Paul Sommers. This work is supported by the National Science Foundation under contract No. PHY-93-21949, PHY-93-22298 and the U.S. Department of Energy.


text:
# 𝜎 Decay at Finite Temperature and Density

## Abstract

$`\sigma `$ decay and its relation with chiral phase transition are discussed at finite temperature and density in the framework of the Nambu-Jona-Lasinio model. The decay rate for the process $`\sigma 2\pi `$ to first order in a $`1/N_c`$ expansion is calculated as a function of temperature $`T`$ and baryon density $`n_b`$. In particular, only when the chiral phase transition happens around the tricritical point, the $`\sigma `$ decay results in a non-thermal enhancement of pions in the final state distributions in relativistic heavy ion collisions.

It is generally believed that there are two QCD phase transitions in hot and dense nuclear matter. One of them is related to the deconfinement process in moving from a hadron gas to a quark-gluon plasma, and the other one describes the transition from the chiral symmetry breaking phase to the phase in which it is restored.

At low energies, there is no obvious reason for the presence of the scalar meson $`\sigma `$ in the study of chiral properties. However, there is a definite need for $`\sigma `$ at finite temperature and density in the investigation of chiral symmetry restoration. For instance, the contribution from $`\sigma `$ to thermodynamics can for most case be neglected, because of its heavy mass, while it plays an important role in the region around the critical point of the chiral phase transition, since there, $`m_\sigma =0`$. Because it is still difficult to extract definite information from the lattice simulations with nonzero baryon density, we need QCD models to investigate the phase transitions at finite density. From the study of several models, the chiral phase transition in the chiral limit with massless pions is of second order at high temperature $`T`$, and of first order at high density $`n_b`$. Therefore, there is a tricritical point $`P`$ in the $`Tn_b`$ plane which separates the first- and second-order phase transitions. In the real world with nonzero pion mass, the second-order phase transition becomes a smooth crossover and the tricritical point $`P`$ becomes a critical end-point $`E`$ of a first-order phase transition. In relativistic heavy ion collisions, for the choice of control parameters like colliding energy and impact parameter such that the freeze-out of the system occurs near the point $`E`$, $`\sigma `$ is one of the most numerous species at the freeze-out, since it is nearly massless. During the expansion of the system, the in-medium $`\sigma `$ mass rises towards its vacuum value and eventually exceeds the $`\pi \pi `$ threshold. As the $`\sigma \pi \pi `$ coupling is large, the decay proceeds rapidly. Since this process occurs after freeze-out, the pions generated by it do not have a chance to thermalize. Thus, it is expected that the resulting pion spectrum should have a non-thermal enhancement at low transverse momentum.

In this letter we discuss the decay process $`\sigma 2\pi `$ at finite temperature and density in the framework of the Nambu-Jona-lasinio (NJL) model which is believed to describe the chiral properties well. This was firstly studied at zero density by Hatsuda and Kunihiro. Our motivation here is to see that in which temperature and density region of freeze-out there is a non-thermal enhancement of pions due to the $`\sigma `$ decay, by calculating the decay rate in the whole $`Tn_b`$ plane.

The two-flavor version of the NJL model is defined through the lagrangian density,

$$L_{NJL}=\overline{\psi }(i\gamma ^\mu _\mu m_0)\psi +G[(\overline{\psi }\psi )^2+(\overline{\psi }i\gamma _5\tau \psi )^2],$$
(1)
where only the scalar and pseudoscalar interactions corresponding to $`\sigma `$ and $`\pi `$ mesons, respectively, are considered, $`\psi `$ and $`\overline{\psi }`$ are the quark fields, $`\tau `$ is the $`SU(2)`$ isospin generator, $`G`$ is the coupling constant with dimension $`GeV^2`$, and $`m_0`$ is the current quark mass.

In an expansion in the inverse number of colors, $`1/N_c`$, the zeroth order (Hartree) approximation for quarks together with the first order (RPA) approximation for mesons gives a self-consistent treatment of the quark-meson plasma in the NJL model. To the lowest order, namely $`\left(1/N_c\right)^{1/2}`$, the Feynman diagram for the decay process $`\sigma 2\pi `$ is sketched in Fig.1. The partial decay rate in the rest frame of $`\sigma `$ is given in terms of the Lorentz-invariant matrix element $`M`$ by

$$\frac{d\mathrm{\Gamma }_{\sigma 2\pi }}{d\mathrm{\Omega }}=\frac{1}{32\pi ^2}\frac{|𝐩|}{m_\sigma ^2}|M|^2,$$
(2)
where $`|𝐩|=\sqrt{\frac{m_\sigma ^2}{4}m_\pi ^2}`$ is the pion momentum, $`m_\sigma `$ and $`m_\pi `$ are respectively $`\sigma `$ and $`\pi `$ masses. using an obvious notation for the $`M`$-matrix element, one has

$$M=g_\sigma g_\pi ^2A_{\sigma \pi \pi }$$
(3)
with

$$A_{\sigma \pi \pi }=T\underset{n}{}e^{i\omega _n\eta }\frac{d^3𝐪}{(2\pi )^3}TrS(𝐪,i\omega _n)\mathrm{\Gamma }_\pi S(𝐪+𝐩,i\omega _n+\frac{m_\sigma }{2})\mathrm{\Gamma }_\pi S(𝐪,i\omega _n+m_\sigma ).$$
(4)
Here $`\omega _n=(2n+1)\pi T,n=0,\pm 1,\pm 2,\mathrm{}`$ are fermionic Matsubara frequencies, $`\mathrm{\Gamma }_\pi =i\gamma _5\tau `$ is the pseudoscalar interaction vertex. The quark propagator is denoted by $`S(𝐪,i\omega _n)=\left(\gamma _0(i\omega _n)\gamma 𝐪+m_q\right)/\left((i\omega _n)^2E_q^2\right)`$, with $`E_q^2=m_q^2+𝐪^2`$. In Eq.(4), $`Tr`$ refers to the trace over color, flavor and spinor indices. The pion-quark and sigma-quark coupling strengths in Eq.(3) are determined in the model via

$`g_\pi ^2(T,\mu )={\displaystyle \frac{\mathrm{\Pi }_\pi (k_0,\mathrm{𝟎};T,\mu )}{k_0^2}}|_{k_0^2=m_\pi ^2},`$
$`g_\sigma ^2(T,\mu )={\displaystyle \frac{\mathrm{\Pi }_\sigma (k_0,\mathrm{𝟎};T,\mu )}{k_0^2}}|_{k_0^2=m_\sigma ^2},`$ (5)
where $`\mathrm{\Pi }_\pi `$ and $`\mathrm{\Pi }_\sigma `$ are the standard mesonic polarization functions for $`\pi `$ and $`\sigma `$. The dynamically generated quark mass $`m_q(T,\mu )`$ and the meson masses $`m_\pi (T,\mu )`$ and $`m_\sigma (T,\mu )`$ are calculated using the usual gap equations in the Hartree approximation for quarks and RPA approximation for mesons. After evaluation of the Matsubara sum in Eq.(4), one has

$`A_{\sigma \pi \pi }(T,\mu )=`$ $`4m_qN_cN_f{\displaystyle }{\displaystyle \frac{d^3𝐪}{(2\pi )^3}}{\displaystyle \frac{f_F(E_q\mu )f_F(E_q\mu )}{2E_q}}\times `$ (6)
$`{\displaystyle \frac{8(𝐪𝐩)^2(2m_\sigma ^2+4m_\pi ^2)𝐪𝐩+m_\sigma ^2/22m_\sigma ^2E_q^2}{(m_\sigma ^24E_q^2)\left((m_\pi ^22𝐪𝐩)^2m_\sigma ^2E_q^2\right)}},`$
where $`f_F(x)=\left(1+e^{x/T}\right)^1`$ is the Fermi-Dirac distribution function.

$`\sigma `$ can decay into neutral and charged pions. By considering the exchange contribution and $`0\theta \pi /2`$ for $`\sigma 2\pi _0`$ due to the identical particle effect, the total decay rate can be written as

$`\mathrm{\Gamma }_{\sigma 2\pi }(T,\mu )`$ $`=\mathrm{\Gamma }_{\sigma 2\pi _0}(T,\mu )+\mathrm{\Gamma }_{\sigma \pi _+\pi _{}}(T,\mu )`$ (7)
$`={\displaystyle \frac{3}{8\pi }}{\displaystyle \frac{\sqrt{\frac{m_\sigma ^2}{4}m_\pi ^2}}{m_\sigma ^2}}g_\sigma ^2g_\pi ^4|A_{\sigma \pi \pi }(T,\mu )|^2\left(1+f_B({\displaystyle \frac{m_\sigma }{2}})\right)^2,`$
where we have taken into account the Bose-Einstein statistics in terms of the distribution function $`f_B(x)=\left(e^{x/T}1\right)^1`$ for the final state pions.

From the comparison of $`\sigma \pi \pi `$ coupling strength $`g_{\sigma \pi \pi }`$ defined through the lagrangian $`L_Ig_{\sigma \pi \pi }\sigma \pi \pi `$ and the $`M`$-matrix element (3) in the NJL model, we have

$$g_{\sigma \pi \pi }(T,\mu )=2M(T,\mu )=2g_\sigma g_\pi ^2A(T,\mu ).$$
(8)
Now let us turn to numerical calculations of the coupling strength $`g_{\sigma \pi \pi }`$ and the decay rate $`\mathrm{\Gamma }_{\sigma 2\pi }`$ as functions of temperature $`T`$ and baryon density $`n_b`$. The relation between $`n_b`$ and the chemical potential $`\mu `$ is rather complicated and is also a function of $`T`$,

$$n_b(T,\mu )=\frac{N_cN_f}{3}\frac{d^3𝐪}{(2\pi )^3}\left[\mathrm{tanh}\frac{1}{2T}(E_q+\mu )\mathrm{tanh}\frac{1}{2T}(E_q\mu )\right],$$
(9)
where the factor $`3`$ reflects the fact that $`3`$ quarks make a baryon. In the calculation we use the dynamical quark mass $`m_q=0.32GeV`$, the pion decay constant $`f_\pi =0.093GeV`$ and the pion mass $`m_\pi =0.134GeV`$ in the vacuum ($`T=\mu =0`$) as input to determine the $`3`$ parameters in the NJL model, namely the coupling constant $`G=4.93GeV^2`$, the momentum cutoff $`\mathrm{\Lambda }=0.653GeV`$ and the current quark mass $`m_0=0.005GeV`$.

We first discuss the temperature dependence of $`g_{\sigma \pi \pi }`$ and $`\mathrm{\Gamma }_{\sigma 2\pi }`$ in the case of zero baryon density which corresponds to the central region of relativistic heavy ion collisions. The coupling strength $`g_{\sigma \pi \pi }`$ is about $`2GeV`$ in the vacuum and almost a constant in the temperature region $`T<0.1GeV`$, and then decreases smoothly with increase of the temperature, as shown in Fig.2. The energy condition for the decay process $`\sigma 2\pi `$ to happen is $`m_\sigma >2m_\pi `$ in the rest frame of $`\sigma `$. At the threshold temperature $`T0.19GeV`$ determined by $`m_\sigma (T)=2m_\pi (T)`$, $`g_{\sigma \pi \pi }`$ has a sudden jump from about $`1.2GeV`$ to zero. In a wide temperature region, the decay rate (7) goes up gradually from its vacuum value $`0.09GeV`$ to the maximum $`0.11GeV`$ at $`T0.17GeV`$, and then it drops down rapidlly but continuously to zero at the threshold temperature. The fact that the maximum value of $`\mathrm{\Gamma }`$ is not in the vacuum, but at $`T0.17GeV`$ is beyond our expectation. It originates from the important Bose-Einstein statistics of the final state pions at high temperature. This can be seen clearly in Fig.2 where we show the decay rates with and without consideration of the Bose-Einstein statistics. The rapid decrease of $`\mathrm{\Gamma }_{\sigma 2\pi }`$ at high temperature is due to the behavior of $`g_{\sigma \pi \pi }`$ and the phase-space factor $`\sqrt{m_\sigma ^2/4m_\pi ^2}/m_\sigma ^2`$.

The density effect on $`g_{\sigma \pi \pi }`$ and $`\mathrm{\Gamma }_{\sigma 2\pi }`$ is shown in Fig.3. Different from the temperature effect in Fig.2 where the statistics of the final state pions is a crucial factor, there is no difference between with and without such statistics when we discuss density effect at low temperature (see Fig.2) which corresponds to the region of the first-order chiral phase transition in the $`Tn_b`$ plane. Instead of the statistics, the chiral restoration dominates the behavior of $`g_{\sigma \pi \pi }`$ and $`\mathrm{\Gamma }_{\sigma 2\pi }`$. At $`T=0`$, the mixing phase of the first-order chiral transition is in the density region $`0.33<n_b/n_0<1.56`$ where $`n_0`$ is the baryon density in normal nuclear matter. At the up limit $`n_b/n_0=1.56,m_\sigma <2m_\pi `$, there is no $`\sigma `$ decay, while at the low limit $`n_b/n_0=0.33,m_\sigma >2m_\pi `$, $`\sigma `$ decay begins with maximum rate. With the decrease of the density, the decay rate and the coupling strength drop down from the maximums to their vacuum values, respectively. The fact that the maximum decay rate and the maximum coupling strength local at the low density limit of the mixing phase is directly due to the first-order chiral transition. The situation for $`T=0.02GeV`$ is similar, but the mixing region becomes narrower, $`0.59<n_b/n_0<1.44`$, and the jump of the decay rate (coupling strength) becomes smaller. As the temperature increases further, the density interval of the mixing phase becomes narrower and narrower, and the jump becomes smaller and smaller. Finally when approaching the end-point $`E`$ of the first-order transition at $`T0.03GeV`$, the interval and the jump disappear. The decay begins at the point $`E`$ and proceeds faster and faster with the decrease of the density. The decay rate goes up from zero at the phase transition point $`E`$ to the maximum at low density. While the coupling strength has still a jump at the point $`E`$, the decay rate is continuous.

To see the parameter dependence of the above discussion, we considered another set of parameters, $`G=5.46GeV^2,\mathrm{\Lambda }=0.632GeV`$ and $`m_0=0.0055GeV`$. It makes shifts of the magnitudes of $`g_{\sigma \pi \pi }`$ and $`\mathrm{\Gamma }_{\sigma 2\pi }`$, but does not change their shapes, especially the behavior around the threshold temperature or the chiral phase transition point.

From the above investigation, we see that when $`\sigma `$ decay into $`2\pi `$ happens in low density or low temperature region, the maximum decay rate does not locate in the vacuum, but is very close to the threshold point or just at the first-order chiral transition point. Only when the decay trajectory passes through the region around the end-point $`E`$ of the phase transition, the position of the maximum decay rate is far from the point $`E`$ where the sigma has minimum mass and therefore can be numerously produced. According to this picture, when the system of a relativistic heavy ion collision evolves in the low density or low temperature region, the $`\sigma `$ decay begins with maximum rate immediately after or even simultaneously at its starting point of creation. Since the $`\sigma `$ decay happens already at the freeze-out, there is no obvious non-thermal enhancement of pions. However, when the freeze-out occurs around the tricritical point, the $`\sigma `$ decay happens far from the freeze-out, there is really a non-thermal enhancement of pions which can be considered as a direct signature of chiral phase transition.
acknowledgments:
This work was supported in part by the Natural Science Foundation of China under grant numbers 19845001 and 19925519.

Figure Captions

Fig.1: The Feynman diagram for $`\sigma `$ decay into $`2\pi `$ to the lowest order in $`1/N_c`$ expansion. The solid lines denote quarks and antiquarks, the dashed lines denote mesons.

Fig.2: The total decay rate $`\mathrm{\Gamma }_{\sigma 2\pi }`$ and the coupling strength $`g_{\sigma \pi \pi }`$ as functions of the temperature at $`n_b=0`$. The thick and thin solid lines indicate the decay rates with and without consideration of Bose-Einstein statistics for the final state pions.

Fig.3: The total decay rate $`\mathrm{\Gamma }_{\sigma 2\pi }`$ (solid lines) and the coupling strength $`g_{\sigma \pi \pi }`$ (dashed lines) as functions of the baryon density scaled by $`n_0`$ for three values of the temperature. The black circle and black lozenge indicate the starting- and end-point of the mixing phase at $`T=0`$ and $`T=0.02GeV`$, respectively.