text:
# Electron (hole) paramagnetic resonance of spherical CdSe nanocrystals.

## Abstract

A new mechanism of electron paramagnetic resonance in spherical zinc-blende semiconductor nanocrystals, based on the extended orbital motion of electrons in the entire nanocrystal, is presented. Quantum confinement plays a crucial role in making the resonance signal observable. The mechanism remains operative in nanocrystals with uniaxially distorted shape. A theoretical model based on the proposed mechanism is in good quantitative agreement with unusual ODMR spectra observed in nearly spherical CdSe nanocrystals.

Nearly two decades of research have explored the origin and potential applications of unusual opto-electronic properties of semiconducting nanocrystals. In particular, quantum size effects on nanometer length scale can lead to unique features in optical properties (,) as compared to the corresponding properties of the bulk solids. Several reports have recently described the development of a new type of colloidal nanocrystals based on CdSe and CdS core coated with an epitaxial shell of a different semiconductor or by organic ligands (-).

The electronic properties of nanocrystals have been investigated extensively during the last decade (see e.g. ,, and references therein). However, there is still significant uncertainty regarding the orbital and spin dynamics of carriers in the presence of an external magnetic field, under the quantum confined conditions of the nanocrystals. In particular, a direct method for determining the carriers’ effective masses in nanocrystals, similar to the cyclotron resonance technique used routinely in bulk semiconductors, is lacking. As will be shown in this paper, however, application of the Optically Detected Magnetic Resonance (ODMR) technique to nanocrystalline materials can provide the missing information by taking advantage of the quantum confinement imposed on the carriers dynamics.

Our recent ODMR study of CdSe and CdS nanocrystals has revealed unusually broad resonance bands which could not be accounted for in terms of the standard spin Hamiltonian model since such an interpretation would yield unreasonably large g-factors ( i.e. in the range of $`410`$ ). In this paper we present a new mechanism of electron paramagnetic resonance in spherical semiconducting nanocrystals, based on the extended orbital motion of electrons in the entire nanocrystal, which can account for the observed resonance bands. In macroscopic crystals paramagnetic resonance of conduction electrons is strongly damped by the coupling to the lattice degrees of freedom ( typical spin-lattice relaxation times $`T_1`$ are in the range of $`10^{10}`$ to $`10^8`$ sec. ()). The dramatic suppression of the relaxation rates in small nanocrystals ( with $`T_110^610^4`$ sec () ) due to the effect of quantum confinement can make these resonances easily observable.

Nanocrystals of CdSe, capped either with tri-octy-phosphine-oxide (TOPO) or epitaxial layers of CdS with core diameter of about 30$`\AA `$, were prepared according to the procedure described by Peng et al. (). The photoluminescence (PL) and ODMR measurements were carried out by immersing the samples in a cryogenic dewar (at 1.4K) and exciting them with a continuous 457.9 nm $`Ar^+`$ laser. The sample was mounted on a special sample probe at the centre of a High-Q resonance cavity, coupled to a microwave (mw) source ( 10 GHz), and surrounded by a superconducting magnet (B). A detailed description of the experimental system is described in Ref.( ). The ODMR spectra were obtained by measuring a change in luminescence intensity, due to the absorption of a modulated magnetic component of a microwave radiation. The latter was plotted (see Fig.1) versus the external magnetic field, B, yielding a magnetic resonance-like spectrum.

The PL spectra of CdSe (TOPO) and CdSe/CdS core-shell samples consist of an exciton band (centred at 2.15 eV), predominantly tunable with the size of the core, and an additional broad band at lower energies (centered at 1.70 eV). The exciton band in both samples did not show resonance phenomena and therefore, the ODMR spectra were selectively recorded around 1.7 eV. The ODMR spectrum of either CdSe(TOPO) or CdSe/CdS nanocrystals consisted of three broad resonance bands (Fig. 1a), ranging between 200-8500 Gauss. The spectrum of the CdSe/CdS core-shell structure showed additional overlapping narrow resonance bands (see inset in Fig. 1a ) between 3000-5000 Gauss. The broad resonances were mainly pronounced upon an application of kHz mw power modulation, while the narrow resonances in CdSe/CdS were observed under a mw modulation of about 100Hz. This distinction suggests the existence of two types of magnetic resonance events, with different characteristic relaxation times. Furthermore, the appearance of the narrow resonances with g-factors close to 2.000 ( i.e. $`1.988`$ and $`1.845`$ ) in contrast to the broad ones, emphasizes the localized nature of the corresponding states within the forbidden band gap. The chemical identification of the corresponding trapping sites are discussed at length in a separate publication (). They are considered as pure spin resonance transitions of trapped electrons and holes.

The broad ODMR bands are reminiscent of cyclotron resonance data with cyclotron effective mass values in good agreement with the well known conduction and valence band effective masses of CdSe. It is evident, however, that well defined cyclotron orbits cannot exist within the nanocrystal’s interior since the minimal cyclotron radius of an electron under the magnetic field used in such experiments is much larger than the size of a nanocrystal. Thus, the observed magnetic resonances cannot be associated with the field induced diamagnetic currents within the nanocrystals and should be assigned to a paramagnetic effect. They are most probably associated with transitions between initial extended electron (hole) states within the entire space of the nanocrystal and final localized hole (electron) states trapped in defect cites (shown schematically in Fig.2). The magnetic resonance events take place during the initial stages.

Focusing on the magnetic resonance processes we follow an approach exploited in Ref. ( ), where the g-factor for acceptor levels in diamond-type semiconductors was calculated. The spherical symmetry of the electron (hole) wavefunctions at the $`\mathrm{\Gamma }`$-point of the Brillouin zone simplifies the calculation considerably, leaving the crystal field and lattice effects to appear only through effective masses and g-factors. Within this approach the spherical symmetry is broken only at the nanocrystal surface, namely when some of the nanocrystal faces do not follow the elementary Wigner-Seitz polyhedron of the zinc-blende structure. Electron and hole are treated throughout the paper as independent particles. This is justified by the large electronic energy separation in a nanocrystal with respect to the electron-hole Coulomb interaction ().

Near the $`\mathrm{\Gamma }`$-point the conduction band of bulk CdSe is nearly parabolic, with isotropic effective mass tensor. Utilizing the effective mass and g-factor of an electron provided by the $`𝐤𝐩`$ theory, the relevant Hamiltonian is written in the form,

$$\widehat{H}=\frac{1}{2m^{}}\left(𝐩\frac{e}{c}𝐀\right)^2+g^{}\mu _B𝐒𝐁+\widehat{H}_{S.O.}$$
(1)
where $`m^{}`$ is the conduction band effective mass at the $`\mathrm{\Gamma }`$ -point, $`g^{}`$ is the corresponding effective g-factor (i.e. corrected by the crystal field), and $`\widehat{H}_{S.O.}`$ is the spin-orbit (SO) interaction.

Selecting the symmetric gauge form of the vector potential $`𝐀=(0,0,\frac{1}{2}Brsin\theta )`$ to describe a constant and uniform magnetic field $`𝐁`$ aligned along the $`z`$-axis, with spherical polar coordinates $`(r,\theta ,\varphi )`$ , Eq.(1) can be rewritten in the form: $`\widehat{H}=𝐩^2/2m^{}+\widehat{H}_{S.O.}+\beta (L_z+g_SS_z)B`$, where the diamagnetic term, $`(eBrsin\theta )^2/8m^{}c^2`$, is neglected. Here $`\beta =\mu _Bm_0/m^{}`$, $`g_S=g^{}m^{}/m_0`$, $`m_0`$ is the free electron mass, and $`L_z`$ and $`S_z`$ are, respectively, the projections of the orbital and spin angular momentum on the z-axis. In the presence of the spin-orbit interaction the eigenstates correspond to the total angular momentum $`𝐅=𝐋+𝐒`$ , and the corrections to the electron energy due to the paramagnetic term are $`\delta E_{el}=g_{F,el}\beta M_FB`$, where $`g_{F,el}`$ is the electronic Lande factor, and $`M_F`$ is the projection of $`𝐅`$ on the z-axis. Using a standard scheme, we find that $`g_F=1+(g_S1)\alpha `$, where $`\alpha =\left[S(S+1)+F(F+1)L(L+1)\right]/2F(F+1)`$.

The corresponding calculation near the valence band edge is considerably more complicated due to the four-fold degeneracy of the band edge at the $`\mathrm{\Gamma }`$-point. The relevant magnetic resonance transitions involve the Zeeman splitting of a hole near this energy. Since the energies involved are much smaller than the SO split-off energy $`\mathrm{\Delta }`$, the latter can be ignored (). Thus, the motion of holes near the band edge in the presence of the magnetic field can be described by Luttinger Hamiltonian (),

$$\widehat{D}=\frac{1}{m_0}\left\{(\gamma _1+\frac{5}{2}\gamma )\frac{k^2}{2}\gamma (𝐤𝐉)^2+(\kappa \frac{\gamma }{2})\frac{e}{c}𝐉𝐁\right\},$$
(2)
where $`\gamma `$, $`\gamma _1`$, and $`\kappa `$ are Luttinger band parameters, $`𝐤=𝐩\frac{e}{c}𝐀`$ and $`𝐉`$ is the set of $`4\times 4`$ matrices of spin $`3/2`$. In the absence of the magnetic field the $`𝐤𝐩`$ Hamiltonian at the $`\mathrm{\Gamma }`$-point is invariant under rotations so that the total angular momentum $`𝐅=𝐋+𝐉`$ is conserved, while the resulting energy eigenvalues are independent of $`M_F`$ . The perturbation theory with respect to the paramagnetic term yields for the Zeeman splitting of a hole near the band edge $`\delta E_h=g_{F,h}\mu _BM_FB`$, where $`g_{F,h}`$ is the Lande factor of a hole.

The calculation of $`g_{F,h}`$ in the general case is complicated by the fact that in the corresponding unperturbed wavefunctions heavy and light hole components are mixed. A good approximation of $`g_{F,h}`$ can be obtained, however, in the limit $`F1`$ , when the hole radial wavefunctions $`R_{F,l}\left(r\right)`$ () can be approximated by their asymptotic forms for which the motion of the heavy and light holes are completely decoupled. Consequently, the Zeeman energies of the holes can be written in a separable form, $`\delta E^{hh}=[\gamma _1\frac{1}{2}\gamma ]\mu _BM_FB,\delta E^{lh}=[\gamma _1+\frac{1}{2}\gamma ]\mu _BM_FB,`$ where ’lh’ and ’hh’ stand for light and heavy hole, respectively. The relevance of these simple limiting expressions to the present study will be discussed below.

The above analysis relied on the spherical symmetry of the nanocrystals under study. TEM pictures () however, showed that the nanocrystals are slightly elongated spheroids or ellipsoids. The effect of this symmetry breaking on the relevant electron and hole paramagnetic states, was estimated by assuming small deviations from the spherical shape and by exploiting the perturbation theory. A uniaxial shape distortion introduces an energy, $`\widehat{H}_{SD}`$, which is typically much smaller than $`\widehat{H}_{SO}`$ , but much larger than the Zeeman spitting$`\widehat{H}_Z`$ (see below). Thus, the perturbative procedure can be carried out in two stages, within the framework of the distorted wave Born approximation.

We consider the case of an electron near the bottom of the conduction band and use a model of a free particle enclosed in a nearly spherical cavity slightly distorted along some axis. The first order correction to energy due to this distortion is given by $`\delta E_{n,l,F,M_F}^{SD}=ϵE_{n,l,F}^0\left[M_F^2/F(F+1)1/3\right]`$, where $`E_{n,l,F}^0`$ is the corresponding unperturbed energy, $`ϵ=2(ab)/(a+b)`$ is the deformation parameter, assumed to be small, $`a`$ and $`b`$ are the principal semiaxes of the ellipsoid. Note that in the latter formula the z-axis is selected along the symmetry axis of the nanocrystal.

Due to this uniaxial distortion the original $`2F+1`$ fold degenerate multiplet splits into a collection of Kramers’ doublets, each of which having the same value of $`|M_F|`$. For the prolate nanocrystals the lowest lying doublet corresponds to $`M_F=\pm 1/2`$ and the highest one has $`|M_F|=F`$. In the oblate case the picture is reversed.

It is important to note that, due to the finite size of the electron wavefunctions ( i.e. $`R34nm`$ ), the minimal value of the unperturbed energy $`E_{n,l,F}^0\mathrm{}^2R^2/2m^{}`$ is typically much larger than the Zeeman energy $`\beta M_FB`$; consequently, the distortion energy, $`\delta E_{n,l,F,M_F}^{SD}`$ , is significantly larger than $`\beta M_FB`$ .

An external static magnetic field $`𝐁`$, lifts the residual degeneracy of the Kramers’ doublets. When the field is not aligned along the nanocrystal’s symmetry axis, one should use the more general form $`\widehat{H}_Z=g_F\beta 𝐅𝐁`$. The corresponding Zeeman splitting energy of the doublets with $`|M_F|=1/2`$ is not given by the familiar formula $`\delta E_Z=g_F\mu _BBcos\theta _0M_F`$, where $`\theta _0`$ is the angle between $`𝐁`$ and the symmetry axis of the ellipsoid, since $`\widehat{H}_Z`$ mixes the zero-order functions of the doublet with $`M_F=\pm 1/2`$. A slightly more complicated expression, $`\delta E_Z=\pm 1/2g_F\mu _BB[cos^2\theta _0+(F+1/2)^2sin^2\theta _0]^{1/2}`$ , is obtained in this case, which reduces however to the $`\theta _0=0`$ expression, $`\delta E_{el}`$ , at $`F=1/2`$ for any $`\theta _0`$. This dependence of the Zeeman splitting energy on $`F`$ and $`\theta _0`$ requires a careful calaculation of the observable position of the magnetic resonance peak, which takes into account the transition matrix element $`M_{1/21/2}`$. It is easy to show that $`\left|M_{1/21/2}\right|=A(F+\frac{1}{2})\left|\mathrm{cos}2\theta _0\right|/2\left[\mathrm{cos}^2\theta _0+(F+\frac{1}{2})^2\mathrm{sin}^2\theta _0\right]^{1/2}`$ , where $`A`$ is constant, independent of $`\theta _0`$ and $`F`$.

Application of the microwave magnetic field in a plane perpendicular to $`𝐁`$ causes transitions of significant intensity between different Zeeman split sublevels satisfying the selection rule $`\mathrm{\Delta }M_F=\pm 1`$. Considering first only the nanocrystals with symmetry axes parallel to $`𝐁`$ ( i.e. where $`\theta _0=0`$ ), it can be seen that among all the transitions satisfying the selection rule, only those originating from the $`M_F=\pm 1/2`$ doublet can be observed at magnetic field $`B`$ in the experimentally accessible range. This is due, on the one hand, to the large distortion energy $`\delta E_{n,l,F,M_F}^{SD}`$ and, on the other hand, to the the Kramers’ degeneracy of the $`M_F=\pm 1/2`$ states in the distorted nanocrystals. Similar behavior is anticipated for the hole transitions, though the corresponding quantitative analysis of the shape distortion is quite involved ().

The present ODMR spectra was monitored at the PL broad band, centered well below the bandgap energy. Thus it is assigned to transitions between excited electrons (holes) within spin-orbit states at the conduction (valence) band, recombining with holes (electrons) in defect levels, as shown schematically in Fig.2. As mentioned earlier, the spin-lattice relaxation time, $`T_1`$ , of the excited electrons and holes is in the range of $`10^4`$sec, while that of the trapped defect is of two orders of magnitude larger, thus suggesting that the spin flipping at the defect site can be considered stationary at the time scale of the spin flip within the spin-orbit manifold. Therefore, the spin splitting of the defect states was ignored in the model shown in Fig.2.

Let us first consider the ODMR band associated with the conduction electrons. For the excitation energy of the continuous $`Ar^+`$ laser used in the experiment only the ground (’s’-state; $`S=1/2,L=0,F=1/2`$) and the first few excited $`F`$-states can be populated. Although the electrons pumped into an excited state relax quickly to the ground state via nonradiative channels, the observed ODMR band may be considered as arising mainly from $`M_F=1/21/2`$ transitions in such an excited state. Indeed, recent experiments have shown that in CdSe nanocrystals the life time, $`\tau _p`$ , of the first excited state (’p’-state; $`S=1/2,L=1,F=1/2`$ and $`S=1/2,L=1,F=3/2`$) is in the range of $`1200\times 10^{12}`$sec. (, ). The corresponding decay rate, $`\tau _p^1`$,which is much larger than the spin relaxation rate, $`T_1^1`$ , is of the same order of magnitude as the resonance microwave frequency, $`\omega _{res}`$ (see Fig.2). These imply that the spin excitation of electrons occupying the ’p’-state by the mw source can be completed before relaxing into the long-lived ( luminescent ) ’s’-state, where their spin polarization is preserved for time much longer than the characteristic recombination time.

Assuming Lorentzian densities of $`F`$-states centered at $`\delta E_Z`$ with frequency spread corresponding to $`\tau _p^1`$ , and using the above expression for $`M_{1/21/2}`$ , we find , after averaging over the angle $`\theta _0`$ , the resonance lines shown in Fig.(1b) for typical values of $`\tau _p`$. For the selected well known value of the electron effective mass $`m^{}=0.13m_0`$ , the calculated resonance position is in a very good agreement with the lower field resonance band observed in the experimental ODMR spectrum. In contrast , the calculated resonance peak due to spin excitations in the ’s’-state appears at a field about 1.1T , well out of the field range shown in Fig.(1). It is interesting to note that the calculated integral intensity of this resonance is found to be much weaker than that of the ’p’ state described above.

Evaluation of $`g_{F,h}`$ for the low-laying hole states has proved to be quite difficult. The corresponding resonance peak positions can be estimated, however, if we note that for $`g_{F,el}`$ in spherical nanocrystal all the levels of the multiplet with $`L>0`$ ( i.e. with the exception of $`L=0,`$ $`F=1/2`$ ) cluster in the close vicinity of the multiplet’s centre. The latter can be found by taking the mathematical limit $`F\mathrm{}`$ in the formula for $`\delta E_{el}`$. Assuming that the resonance lines of holes follow a similar pattern and using the above formulae for $`\delta E^{hh}`$ and $`\delta E^{lh}`$ together with the valence band parameters reported in ( ) the centres of the respective multiplets are located at $`B_{hh}=h\nu /(\gamma _1\frac{1}{2}\gamma )\mu _B`$, $`B_{lh}=h\nu /(\gamma _1+\frac{1}{2}\gamma )\mu _B`$, yielding $`B_{hh}`$ $`=0.44`$ Tesla and $`B_{lh}`$ $`=0.33`$ Tesla respectively. The value of $`B_{lh}`$ is in close agreement with the maximum of the middle broad band shown in Fig.1. However, $`B_{hh}`$ deviates slightly from the highest field band maximum, which is not surprising in view of the approximate nature of the present calculation.

In conclusion it was shown that the proposed mechanism of electron (hole) paramagnetic resonance in nearly spherical zinc-blende semiconductor nanocrystals accounts well for the unusual features observed in the ODMR spectra of $`CdSe`$ nanoparticles. The observability of the corresponding resonance bands is due to the small size (i.e. in the range of several nanometers) of the crystallites as well as to the Kramers’ degeneracy of the energy levels, which cannot be lifted by any shape distortion of electrostatic origin.

The authors express their deepest gratitude to Prof. G. Kventsel, V. Zhuravlev and R. Guliamov for useful discussions. This research was supported by the Israel Science Foundation, grant no. 36599-12.5


format: text

text:
# Uninfected random walkers in one dimension

## I Introduction

Random walks model a host of phenomena and find application in virtually all areas of physics . In this paper we consider a one-dimensional (1D) system of non-interacting, unbiased diffusing walkers $`A\mathrm{}\mathrm{}A`$, with random initial conditions. Recently, we have studied the persistence properties of this model . The persistence probability $`P(t)`$ has been widely studied in recent years . In the present context $`P(t)`$ is defined as the fraction of sites which have never been visited by a diffusing particle up to time $`t`$. Such sites are termed “persistent sites”. We showed  that, in the case of the $`A\mathrm{}\mathrm{}A`$ processes,

$$P(t)=\mathrm{exp}[(16\rho _0^2Dt/\pi )^{1/2}],$$
(1)
where $`\rho _0`$ is the particle density and $`D`$ is the diffusion constant. The fraction of unvisited sites in this model can equally be thought of as the surviving fraction of immobile particles in the presence of randomly diffusing traps. At $`t=0`$, consider each site to contain a static particle, $`A`$, so that $`\rho _A(0)=1`$. Then for $`t>0`$, let particles of another type, $`B`$, randomly distributed on the lattice with density $`\rho _0`$, diffuse though the system with the $`B`$ particles acting as perfect traps for the $`A`$ particles so that $`A+BB`$. This model is just the well-studied “scavenger” model of Redner and Kang . The fraction of unvisited sites is then just $`\rho _A(t)`$.

The focus of this paper is to explore the relation between $`P(t)`$ and the fraction of walkers $`U(t)`$ that have never encountered another walker up to time $`t`$. To facilitate our discussion, we introduce the following definition. We consider all walkers to be uninfected at time $`t=0`$. The walkers are, however, considered to be mutually infectious, so that for $`t>0`$ any contact between walkers leads to mutual contamination. We therefore define an uninfected walker to be one which has never encountered (occupied the same site as) another walker. We can then refer to $`U(t)`$ as the fraction of uninfected walkers. The fraction of uninfected walkers is, in fact, isomorphic to the survival probability of $`A`$ particles in the reaction $`A+BB`$, with equal mobility for both species . At late times, one expects the asymptotics of this model to be identical to those of the $`A+B\mathrm{}`$ diffusion-limited process with $`\rho _A(0)\rho _B(0)`$. For this process it was proved rigorously by Bramson and Lebowitz  that

$$\rho _A(t)\mathrm{exp}(\lambda t^{1/2})$$
(2)
in one-dimension, with unknown constant $`\lambda `$. Therefore we expect identical asymptotic behavior for $`U(t)`$. The kinetics of the two species annihilation model has received considerable attention over the past two decades and is largely understood . It is well appreciated that in this model, there exists an upper critical dimension $`d_c`$, below which spatial fluctuations in the initial distribution of the reactants play a significant role in the evolution of the density of the particles. This dependence on the microscopic fluctuations invalidates traditional approaches such as the mean-field approximation which yields $`\rho (t)t^1`$. It was first shown by Toussaint and Wilczek  that, for equal initial densities ($`\rho _A(0)=\rho _B(0)`$) of randomly diffusing $`A`$ and $`B`$ particles, the walker density decays anomalously, according to $`\rho _A(t)\sqrt{\rho _A(0)}t^{d/4}`$ for $`dd_c=4`$, assuming that the particles are distributed homogeneously at $`t=0`$. In the case of this model we argued heuristically, on the basis of a toy model  that,

$$P(t)\mathrm{exp}[C(\rho _0^2Dt)^{1/4}]$$
(3)
where $`C`$ is a constant and $`\rho _0=\rho _A(0)+\rho _B(0)`$. We extend our study to include an investigation of the fraction of uninfected walkers in this model in one-dimension.

The fraction of uninfected walkers, $`U(t)`$, is related to the fraction of persistent sites, $`P(t)`$, in the following manner. In evaluating $`P(t)`$, one considers particles with diffusion constant $`D`$ and addresses the probability that a given site $`x`$ has never been visited by any walker up to time $`t`$. Equivalently, one could consider the site $`x`$ to contain a particle with diffusion constant $`D^{}=0`$ and address the probability that the particle remains uninfected (i.e. unvisited). However, if the particle initially at position $`x`$ is given a diffusion constant equal to that of the other particles in the system i.e. $`D^{}=D`$, then the probability that it “persists” up to time $`t`$ is the probability that it remains uninfected. This suggests that the number of uninfected particles is a function of the diffusion constants in the problem, with the limiting case $`D^{}=0`$ corresponding to the standard persistence problem. Clearly, in order to avoid contamination, remaining stationary is a more effective strategy than diffusing. Therefore, we anticipate that the fraction of uninfected walkers decays more rapidly than the fraction of uninfected (unvisited, i.e. persistent) sites. This faster decay of $`U(t)`$ is clearly observed in the simulations.

Given the similar decay of $`P(t)`$ in the $`A\mathrm{}\mathrm{}A`$ and $`A+B\mathrm{}`$ processes, these models make a sensible combination to study together. We discover that in both cases, $`U(t)[P(t)]^\gamma `$ where $`1<\gamma 1.39`$. In fact the relation $`U(t)=k[P(t)]^\gamma `$, with $`k`$ close to unity and $`\gamma =1.39(1)`$, describes the data rather accurately in both models over the whole range of times studied, although there are indications in both models that the asymptotic value of $`\gamma `$ is smaller than $`1.39`$. We discuss the $`A\mathrm{}\mathrm{}A`$ model in Sec. II and follow with the $`A+B\mathrm{}`$ reaction in Sec. III. We conclude in Sec. IV with a discussion and summary of the results.

## II The $`A\mathrm{}\mathrm{}A`$ model

We consider a system of non-interacting, unbiased, diffusing walkers, with random initial conditions. Our numerics are performed on a 1D lattice of size $`L=10^7`$ with lattice constant $`a=1`$ and periodic boundaries. At $`t=0`$, $`N`$ random walkers are randomly distributed on the lattice with a maximum of one particle per site so that $`U(0)=1`$ and $`P(0)=1\rho _0`$, where $`\rho _0=\rho (0)=N/L`$. Clearly, $`P(0)1`$ is a consequence of our lattice description and in the continuum limit $`\rho _00`$, where the mean distance between walkers is much larger than the lattice spacing, $`P(0)1`$. Therefore, in order to place the fraction of persistent sites and the fraction of uninfected walkers on an equal footing in our numerical study, we define,

$$p(t)=\frac{P(t)}{P(0)}$$
(4)
so that both $`U(0)`$ and $`p(0)`$ are normalized to 1. We now use $`p(t)`$ in the description of our numerical results.

Our model is updated using the direct method  i.e. at each computational step, a particle is picked at random and shifted with probability $`D=1/2`$ to a neighboring site, where $`D`$ is the diffusion constant. For $`t>0`$, any site which is visited by a walker becomes nonpersistent for all remaining time $`t`$. Mutual infection occurs when two or more particles simultaneously occupy a single lattice site. For each jump made by a particle, time $`t`$ is increased by $`dt=1/N`$, where $`N`$ is the current number of particles in the system (which is constant in the present model). Our results are averaged over 100 independent runs.

Dimensional analysis demands that, on a continuum, the fraction of unvisited sites is a function of the dimensionless combination of parameters in the problem, i.e. $`P(t)=f(\rho _0^2Dt)`$. Similarly, for the fraction of uninfected walkers, $`U(t)=g(\rho _0^2Dt)`$, where $`g`$ is some other function. However, given that our numerical simulations take place on a discrete lattice, we expect the scaling of $`P(t)`$ and $`U(t)`$ as a function of $`\rho _0^2Dt`$ to become strictly true only in the continuum limit, $`\rho _00`$. In order to approximate this limit as closely as possible, we choose as our densities $`\rho _0=0.007`$, 0.005 and 0.003, so that $`\rho _01`$.

The dimensionless combination of parameters, $`\rho _0^2Dt`$, also defines a natural characteristic time scale in the problem, namely $`t^{}=1/\rho _0^2D`$. Therefore one expects to observe the asymptotic behavior of the system only after a time $`tt^{}`$ has been reached. However, as far as the persistent properties are concerned, Eq. (1) is an exact relation in the continuum limit for all $`t`$. Therefore, we expect to observe $`P(t)`$ approaching asymptopia very quickly in our numerical data. Clearly, in the discrete description of the model used in our simulations, Eq. (1) breaks down at the earliest times since, for very small $`t`$, lattice effects play a dominant role. Eq. (2), however, which holds for $`U(t)`$ is a truly asymptotic relation. Unfortunately, the extremely fast decay of $`U(t)`$ leads to small-number effects at very early times in our simulations and we find the limit $`tt^{}`$ unattainable with good statistics. Therefore, we limit ourselves to studying the regime up to $`t_{max}t^{}/2`$.

The fraction of unvisited sites decays, in the continuum limit, exactly as the stretched exponential form, Eq. (1) . In order to make a direct comparison between the fraction of persistent sites, $`p(t)`$, and the fraction of uninfected walkers, $`U(t)`$, we present these two quantities on a log-linear plot in Fig. 1 , using $`(\rho _0^2Dt)^{1/2}`$ as abscissa.

Fig. 1 shows clearly that the data scale as a function of the dimensionless quantity $`\rho _0^2Dt`$. $`U(t)`$ exhibits behavior characteristic of Eq. (2), with a faster decay rate than for the fraction of persistent sites i.e. $`\lambda >(16\rho _0^2D/\pi )^{1/2}`$. The data seem consistent with

$$U(t)K[p(t)]^\gamma ,$$
(5)
with $`K1`$ and $`\gamma >1`$. The log-log plot of $`U(t)\text{ vs. }p(t)`$ in Fig. 2 indeed shows an almost linear relationship between $`\mathrm{ln}U(t)`$ and $`\mathrm{ln}p(t)`$. In order to evaluate $`\gamma `$ and $`K`$ we performed a linear regression on the data in Fig. 2 in the region indicated by the arrows, thereby avoiding initial transients associated with our lattice description, and the onset of statistical fluctuations between different runs of the data. We summarize our results for $`\gamma `$ and $`K`$ in Table 1.

| $`\rho _0`$ | $`\gamma `$ | $`K`$ |
| --- | --- | --- |
| 0.007 | 1.391(5) | 0.990(3) |
| 0.005 | 1.391(5) | 0.991(3) |
| 0.003 | 1.391(5) | 0.991(3) |

TABLE 1. Numerical values of $`\gamma `$ and $`K`$, where $`U(t)K[p(t)]^\gamma `$, for various $`\rho _0`$, in the $`A\mathrm{}\mathrm{}A`$ model.

The data are sufficiently good to merit closer inspection. To check how far our numerics probe the asymptotic regime, we plot $`(\rho _0^2Dt)^{1/2}\mathrm{ln}p(t)`$ and $`(\rho _0^2Dt)^{1/2}\mathrm{ln}U(t)`$ against $`\mathrm{ln}(\rho _0^2Dt)`$ in Fig. 3. The nature of this plot greatly expands the early time regime. Whereas the data-collapse for different densities is excellent in Fig. 1, a systematic splitting of the curves as a function of $`\rho _0`$ is now clearly evident. This feature is a consequence of the discreteness of the lattice. In fact the behavior of $`p(t)`$ is well described, for all $`t`$, by an exact lattice calculation  which predicts

$$p(t)=\mathrm{exp}[\rho _0af(Dt/a^2)],$$
(6)
where $`a`$ is the lattice spacing. The function $`f(x)`$ has the limiting behavior $`f(x)2x`$ for $`x0`$ and $`f(x)(16x/\pi )^{1/2}`$ for $`x\mathrm{}`$, the latter reproducing the continuum limit result, Eq. (1). Indeed, the structure of Eq. (6) shows that the continuum limit, $`a0`$ at fixed $`t`$, and the asymptotic limit, $`t\mathrm{}`$ at fixed $`a`$ are the same. Eq. (6) also shows that the quantity $`(\rho _0^2Dt)^{1/2}\mathrm{ln}p(t)`$ is, for fixed $`a`$ and $`D`$, a function of $`t`$ only, i.e. independent of $`\rho _0`$. For small $`t`$ it behaves as $`t^{1/2}`$, accounting for the small-$`t`$ behavior in Fig. (3). Plotting against $`\mathrm{ln}(\rho _0^2Dt)`$ displaces the curves for different densities in the horizontal direction, as is clear in the figure. The whole region where the curves are split in Fig. 3 corresponds to the top left-hand corner of Fig. 1. If the data are plotted against $`\mathrm{ln}t`$ in Fig. (3), instead of $`\mathrm{ln}(\rho _0^2Dt)`$, the data for different densities recollapse, but the onset of the continuum limit at late times not as clear.

For large $`t`$, where the curves for the different densities converge in these scaled variables, the persistence data clearly approach the limiting value $`4/\sqrt{\pi }`$, consistent with the continuum limit displayed in Eq. (1). By contrast, Fig. 3 shows that all of our data for $`U(t)`$ is merely preasymptotic. The data reach a maximum and then appear to monotonically decrease in the region where the different curves converge (the continuum limit). Our results are consistent with the recent simulations of Mehra and Grassberger  who probe $`U(t)`$ to much smaller values than achieved here. The slow approach to asymptopia for $`U(t)`$ is not understood and, although our numerics seem roughly consistent with Eq. (2), a numerical determination of $`\lambda `$ appears to be a hard problem. In the light of Fig. 3, Fig. 2 and its associated results (Table 1) are somewhat misleading. Clearly, over the times covered by our simulations, $`\gamma `$ can not be a constant since $`\lambda `$ is not.

In Fig. 4 we plot the effective exponent $`\gamma (t)=\mathrm{ln}U(t)/\mathrm{ln}p(t)`$ directly as a function of time. The data show that $`\gamma `$ decreases monotonically with time even in the “continuum regime” where the curves for different densities converge. For $`t\mathrm{}`$, $`\gamma `$ must tend to a constant if Eq. (2) is correct. We note that since the fraction of uninfected walkers decreases more rapidly than the fraction of persistent sites, $`\gamma >1`$. Therefore, if the trend in the data continues, we can bound $`\gamma `$ by $`1<\gamma <1.39`$: Asymptotically, $`U(t)[P(t)]^\gamma `$ for some fixed $`\gamma `$ where $`1<\gamma <1.39`$.

Let us consider once again the probability of being uninfected of a particle at site $`x`$, with diffusion constant $`D^{}`$, and diffusing particles to its left and right with diffusion constant $`D`$. Then, we suggest that the fraction of uninfected walkers can be expressed as,

$$U(t)\mathrm{exp}[(\rho _0^2Dt)^{1/2}F(D^{}/D)]$$
(7)
for large $`t`$, consistent with the result of Bramson and Lebowitz (Eq. (2)). $`F(x)`$ is some unknown function of the ratio of the diffusion constants in the system. For $`D^{}=0`$, the problem reduces to the familiar study of unvisited sites. Therefore $`F(0)=(16/\pi )^{1/2}`$ and $`U(t)=P(t)`$. The case $`D^{}=D`$ corresponds to equally mobile particles, and the particle at site $`x`$ becomes infected more quickly than if it had remained stationary. In this case $`F(1)=\gamma (16/\pi )^{1/2}`$ giving $`U(t)[P(t)]^\gamma `$ where, according to our numerics, $`1<\gamma <1.39`$. Although we have explored the values $`D^{}=0`$ and $`D^{}=D`$ only, we conjecture that Eq. (7) holds for any ratio of the diffusion constants.

## III The $`A+B\mathrm{}`$ Model

We now address the fraction of uninfected walkers, $`U(t)`$, in the two species annihilation model. Our simulations are performed on a 1D lattice of size $`L=10^7`$ with periodic boundary conditions. At $`t=0`$, exactly equal numbers, $`N_A(0)=N_B(0)`$, of $`A`$ and $`B`$ particles are randomly distributed on the lattice with a maximum of one particle per site, such that at large scales, both densities $`\rho _A(0)`$ and $`\rho _B(0)`$ are initially homogeneous. We define $`\rho _A(0)=\rho _B(0)=N_A(0)/L`$ and $`\rho _0=\rho _A(0)+\rho _B(0)`$. Both species are also given the same diffusion constant $`D_A=D_B=D=1/2`$. Our model is then updated in the same manner as for the $`A\mathrm{}\mathrm{}A`$ model described in Sec. II, i.e. $`dt=1/N(t)`$, where $`N(t)=N_A(t)+N_B(t)`$, is the total current number of particles in the system. Infection occurs when two or more particles occupy a single lattice site, but we impose an instantaneous reaction ($`A+B\mathrm{}`$) so that each lattice contains only one type of particle.

It is well known that in the two species annihilation model in one dimension, there is a effective repulsion between the $`A`$ and $`B`$ particles that favors segregation into single species domains . This domain coarsening leads to early time transients for both $`P(t)`$ and $`U(t)`$. Although using a high initial density serves to accelerate the systems progression into the asymptotic regime, in doing so, not only is the continuum limit obscured, but we also find that the extremely fast decay of both $`P(t)`$ and $`U(t)`$ leads to small-number effects at very early times. Therefore we choose $`\rho _01`$. For consistency, we select our parameters to be the same as those studied in the $`A\mathrm{}\mathrm{}A`$ process, i.e. $`\rho _0=0.007`$, 0.005 and 0.003, and $`\rho _0^2Dt_{max}0.5`$ where $`t_{max}`$ is the maximum number of time steps in the simulations. All the results are averaged over 100 runs.

The format of our analysis is much the same as in the preceding section. We observe that an algebraic relation, $`U(t)[p(t)]^\delta `$, between $`U(t)`$ and $`p(t)`$ also holds approximately in this case over the range of times studied. Surprisingly, $`\delta `$ seems to have a similar value to $`\gamma `$ in the $`A\mathrm{}\mathrm{}A`$ model, i.e. $`1<\delta 1.39`$. For the two-species annihilation process we suggested in our earlier paper, on the basis of a toy model, that the fraction of unvisited sites decays asymptotically according to the stretched exponential form, Eq. (3) . In the manner of Fig. 1, to make a direct comparison between the fraction of persistent sites $`P(t)`$, and the fraction of uninfected walkers $`U(t)`$, we present these two quantities on a log-linear plot in Fig. 5, this time using $`(\rho _0^2Dt)^{1/4}`$ as abscissa, as suggested by Eq. (3). Fig. 5 shows clearly that the data scale excellently as a function of the dimensionless quantity $`\rho _0^2Dt`$. We appreciate, however, that the stretched exponential behavior of Eq. (3) is not particularly well realized. This may be due to the prevalence of large initial transients (we presented more convincing evidence of the exponential behavior in our earlier paper ). In particular, the asymptotic behavior represented in Eq. (3) is only expected for large values of $`\rho _0^2Dt`$, unlike Eq. (1) for $`P(t)`$, which we have shown holds for all $`t`$ on the continuum . Nonetheless, even in the regime we have studied it is clear from Fig. 5 that the fraction of uninfected walkers has a very similar form of decay to the fraction of unvisited sites, but with a faster decay rate. In Fig. 6, therefore, we plot $`\mathrm{ln}U(t)`$ against $`\mathrm{ln}p(t)`$ to see whether there is a simple relationship between $`U(t)`$ and $`p(t)`$.

The approximate linearity of the data in Fig. 6 suggests $`U(t)k[p(t)]^\delta `$. To determine $`\delta `$ and $`k`$, we performed a linear regression on the data in the region indicated by the arrows, thereby avoiding initial transients associated with our lattice description, and the onset of statistical fluctuations between different runs of the data. We summarize our results for $`\delta `$ and $`k`$ in Table 2. Notice that $`\delta `$ and $`k`$ agree remarkably well with $`\gamma `$ and $`K`$ recorded in Table 1 for the $`A\mathrm{}\mathrm{}A`$ process.

| $`\rho _0`$ | $`\delta `$ | $`k`$ |
| --- | --- | --- |
| 0.007 | 1.390(5) | 0.990(3) |
| 0.005 | 1.390(5) | 0.993(3) |
| 0.003 | 1.386(5) | 0.991(3) |

TABLE 2. Numerical values of $`\delta `$ and $`k`$, where $`U(t)k[p(t)]^\delta `$, for various $`\rho _0`$, in the $`A+B\mathrm{}`$ model.

In analogy to Fig. 3, we plot $`(\rho _0^2Dt)^{1/4}\mathrm{ln}p(t)`$ and $`(\rho _0^2Dt)^{1/4}\mathrm{ln}U(t)`$ against $`\mathrm{ln}(\rho _0^2Dt)`$ in Fig. 7.

Early-time lattice effects, where the different densities do not overlap, are evident in Fig. 7. It is possible that all of our data for $`p(t)`$ and $`U(t)`$ are within the preasymptotic regime, since the $`p(t)`$ curve has not yet saturated to the late-time asymptotic values impied by Eq. (3). Note, however, that if the curves in Fig. 7 were asympotically linearly increasing, rather than eventually saturating to a constant, the corresponding forms for $`p(t)`$ and $`U(t)`$ would be $`\mathrm{exp}[C(\rho _0^2Dt)^{1/4}\mathrm{ln}(\rho _0^2Dt)]`$, with $`C`$ a constant (different for $`p`$ and $`U`$). This form differs only by the logarithmic term from that predicted in  on the basis of a toy model. It should be noted that neither our present data nor that presented in  can definitively rule out such logarithmic corrections, and the asympotic form of $`p(t)`$, as well as $`U(t)`$, cannot be regarded as being definitively established. Throughout the preasymptotic regime, however, the relation $`U(t)=K[p(t)]^\delta `$, with $`\delta =1.39(1)`$, holds rather well. In Fig. 8 we plot the effective exponent $`\delta (t)=\mathrm{ln}U(t)/\mathrm{ln}p(t)`$ directly as a function of $`\mathrm{ln}(\rho _0^2Dt)`$.

Fig. 8 is remarkably similar to Fig. 4 for the $`A\mathrm{}\mathrm{}A`$ process. Indeed, except for the earliest times, the corresponding curves differ by less by $`1\%`$ everywhere. In a similar way we can bound the asymptotic $`\delta `$ (assuming that the late-time trend evident in Fig. 8 continues) by $`1<\delta <1.39`$. Since we have argued elsewhere , using a toy model, that the asymptotic behavior of $`P(t)`$ is described by Eq. (3), we conjecture that the asymptotics of $`U(t)`$ are given by

$$U(t)\mathrm{exp}[\delta C(\rho _0^2Dt)^{1/4}],$$
(8)
where $`C`$ is a constant and $`1<\delta <1.39`$. For the $`A\mathrm{}\mathrm{}A`$ model, we expressed \[Eq. (7)\] $`U(t)`$ in terms of the ratio $`D^{}/D`$, where $`D^{}`$ is the diffusion constant of a tagged particle, $`D`$ is the diffusion constant of the other particles, and $`U(t)`$ is the probability that the tagged particle remains uninfected at time $`t`$. The analogous form for the $`A+B\mathrm{}`$ process is naturally,

$$U(t)\mathrm{exp}[(\rho _0^2Dt)^{1/4}G(D^{}/D)]$$
(9)
where $`G(x)`$ is some unknown function. As before, the case $`D^{}=0`$ corresponds to the problem of unvisited sites, so that $`G(0)=C`$ and $`U(t)=P(t)`$. For $`D^{}=D`$, the case studied above, our conjecture implies $`G(1)=\delta C`$. However, we stress that a logarithmic correction of the form discussed above cannot be ruled out by the data. The main motivation for Eq. (9) is the toy model presented in , applied to the case $`D^{}=0`$.

## IV Discussion and Summary

In this paper we have studied two distinct models: (i) a one-dimensional system of non-interacting, unbiased diffusing walkers $`A\mathrm{}\mathrm{}A`$, with random initial conditions and (ii) the $`A+B\mathrm{}`$ diffusion limited process in one dimension with equal initial densities of $`A`$ and $`B`$ particles distributed homogeneously at $`t=0`$. We have reduced the familiar study of persistent sites $`P(t)`$, to a limiting case in the study of uninfected walkers and shown numerically that in both models,

$$U(t)[P(t)]^\gamma ,$$
(10)
where $`\gamma 1.39`$ for the time regimes covered by our data and $`1<\gamma <1.39`$ for asymptotically large time. The common feature of the $`A\mathrm{}\mathrm{}A`$ and $`A+B\mathrm{}`$ processes is that the fraction of unvisited sites decays according to a stretched exponential (possibly with logarithmic corrections for $`A+B\mathrm{}`$). In our earlier paper , we argued heuristically, on the basis of a toy model, that if the density of walkers decays asymptotically as $`t^\alpha `$, then the fraction of unvisited sites $`P(t)`$ decays with a stretched exponential form if $`\alpha <1/2`$. We showed that

$$P(t)\mathrm{exp}[A(\alpha )\rho _0^2Dt]^{(1/2)\alpha },\alpha <1/2$$
(11)
where $`A`$ is some function of $`\alpha `$. This result is consistent with the results for models studied here. Remarkably, Eq. (10) seems to hold for both the $`A\mathrm{}\mathrm{}A`$ and the $`A+B\mathrm{}`$ processes, with approximately the same bounds on $`\gamma `$. Furthermore, the effective time-dependent value of $`\gamma `$ is very similar for both models (see Figs. 4 and 8). An obvious goal for the future is to try to obtain a theoretical understanding of this simple relationship, and test it against other models with $`\alpha <1/2`$.

For both models studied here we expressed the probability that a particle at site $`x`$ is uninfected in terms of its diffusion constant $`D^{}`$ and the diffusion constant $`D`$ of the other particles \[Eqs. (7) and (9)\], noting that $`P(t)`$ corresponds to the limit $`D^{}=0`$. Using Eq. (11) , we can generalize this description to any system of randomly diffusing particles for which the density decays as $`t^\alpha `$ with $`\alpha <1/2`$. We conjecture that, in this case, the fraction of uninfected walkers can be expressed, asymptotically, as

$$U(t)\mathrm{exp}[(\rho _0^2Dt)^{(1/2)\alpha }R(D^{}/D,\alpha )],$$
(12)
with $`U(D^{}=0)=P(t)`$ and $`R`$ some unknown function. Clearly, there is much scope to test the validity of this form since we have only considered $`\alpha =0,1/4`$ and $`D^{}=0,D`$. Note that for systems with $`\alpha >1/2`$, our toy model predicts  that $`P(t)`$ approaches a non-zero constant, while $`U(t)`$ certainly behaves differently since it is bounded above by $`\rho (t)`$ and must therefore vanish for large $`t`$. More interesting, however, is the borderline case where the particle density falls off like $`t^{1/2}`$. In this case, our toy model predicts that $`P(t)`$  decays as a power-law, $`P(t)t^\theta `$. A system where such behavior is observed is the well studied $`q`$-state Potts model . The results of our study of uninfected walkers in this model follows in a separate paper .

## V ACKNOWLEDGMENT

We thank Peter Grassberger for useful comments and pointing out references , and Richard Blythe for useful discussions. This work was supported by EPSRC (UK).


text:
# Summary of session on jets and inclusive hadron production

## 1 Introduction

Measurements of jets and hadron production at both LEP and HERA provide detailed investigations of many facets of Quantum Chromodynamics (QCD). They provide constraints on the parton densities in both the photon (LEP and HERA) and proton (HERA) in a complementary way to the inclusive measurements of the structure functions, $`F_2^\gamma `$ and $`F_2^\mathrm{p}`$. The cross sections are also sensitive to the dynamics of the hard sub-processes and perturbative QCD (pQCD) calculations. The transition of partons into hadrons can also be measured, thereby bridging the gap between pQCD calculations in terms of partons and the hadrons seen in detectors.

A lot of precise data exist in jet and hadron production from LEP and HERA, some of which are here discussed. For some results all of the available data from LEP has been analysed and at HERA the interpretation of many measurements is becoming limited by large theoretical uncertainties. To ensure that the data has a legacy, global fits to parton densities and Monte Carlo (MC) tuning need to be done incorporating much of this data. This will improve future measurements at both HERA and other colliders, in particular the LHC where very precise knowledge of the proton is needed.

## 2 Summary of results

### 2.1 Hadron production

Measurements of hadron production testing both the fragmentation  and production mechanism  were presented. At LEP, the L3 collaboration have measured $`\pi ^0`$ and $`K_s^0`$ production and compared them with predictions from NLO calculations, MC models and an exponential fall-off as shown in Fig. 1. The cross sections as a function of the transverse momentum of the hadron, $`p_t`$, are described by an exponential at low $`p_t`$ and by a power law function, $`Ap_t^B`$, at values larger than 1.5 GeV. The value of the power B is compatible with 4 for both $`\pi ^0`$ and $`K_s^0`$ mesons as expected in $`22`$ hard parton scattering . Both the MC and NLO predictions give qualitatively similar results. At low $`p_t`$ they describe the data well, but fail at high $`p_t`$ in the cross section for $`\pi ^0`$ production. It would be of interest to see the distribution to higher values for the $`K_s^0`$ cross section, but this will be unattainable at LEP.

At HERA large samples of $`K_s^0`$ mesons exist, extending up to transverse momenta of $`8`$ GeV, which could verify the preference of the power law also for this meson. The production of $`K_s^0`$ and $`\mathrm{\Lambda }`$ hadrons has been considered at HERA and their relative rates compared to MC predictions. When the MC predictions are normalised to the measured $`K_s^0`$ production cross section, predictions from Herwig overestimate and from Pythia underestimate, the $`\mathrm{\Lambda }`$ productions cross section.

### 2.2 Prompt photon production

The description of prompt photon production by QCD at hadron colliders has long been an issue of interest. The current status of comparison between data and theory is shown in Fig. 2 in which clear discrepancies are seen over a large range in energy . In particular, at the Tevatron, the data from both D0 and CDF is not well described, particularly at the lower transverse energies.

Explanations for this discrepancy in prompt photon production vary, for example the need for resummed calculations or extra $`k_T`$ in the colliding hadron being invoked. In this light, ZEUS have made an extraction of the $`k_T`$ in the proton using a sample of events with a prompt photon and an approximately back-to-back jet. This was extracted by considering variables sensitive to $`k_T`$ and fitting the Pythia MC predictions to these distributions. The value found was consistent with a trend of increasing $`k_T`$ with increasing centre-of-mass energy. However, a recent NLO calculation  of prompt photon and jet production can describe the $`k_T`$-sensitive variables, with no extra $`k_T`$ required. In fact, calculations from the same group  can describe most of the hadron-hadron measurements shown in Fig. 2. Those data not described have found to be inconsistent with other measurements .

### 2.3 Jet production

Jet production is being extensively studied in $`\gamma \gamma `$ collisions at LEP  and $`ep`$ collisions at HERA . The measurements investigate the detailed dynamics of QCD, the parton densities in the photon and proton and allow measurements of the strong coupling constant, $`\alpha _s`$.

The dijet production cross section in $`\gamma \gamma `$ collisions is sensitive to the gluon content of the photon at LO. To try and constrain this, the cross section has been measured by both the DELPHI and OPAL collaborations ; the measurements are consistent with each other.

Using a luminosity of $`600`$ pb<sup>-1</sup>, some results from OPAL are shown in Fig. 4 compared to prediction from NLO and the Pythia MC. In the region of high-$`x_\gamma ^\pm `$, the NLO and Pythia MC prediction describe the data well. For the case of single resolved events at low-$`x_\gamma ^\pm `$, the Pythia prediction underestimates the data by about $`20\%`$. The cross section in $`x_\gamma `$ is very sensitive to hadronisation effects and needs to be studied further. The accurate data can be used in future fits to constrain the parton and in particular the gluon density in the photon.

To make dijet measurements, the cuts on the jets have to be chosen carefully so as not to be infra-red sensitive. This has been extensively studied at HERA. The cuts on the two jets have to be sufficiently different to allow for soft-gluon emission. This region can be seen in Fig. 5 at high values of $`E_{T,2}^{\mathrm{jet},\mathrm{cut}}`$, where the NLO prediction has very small uncertainties and overshoots the measured cross section. At low values of $`E_{T,2}^{\mathrm{jet},\mathrm{cut}}`$, the data and NLO have similar shapes for high-$`x_\gamma ^{\mathrm{obs}}`$ and very different shapes for low-$`x_\gamma ^{\mathrm{obs}}`$. Although the NLO and data agree within the uncertainties, the inclusion of higher orders would require a significant change in shape. This could also indicate problems with the photon parton densities, and higher order calculations would help further constrain them. It should be noted that the predictions from the Herwig MC describe the shape of the data well. It is normalised to the total cross section for the full $`x_\gamma ^{\mathrm{obs}}`$ region and $`E_{T,2}^{\mathrm{jet},\mathrm{cut}}=`$11 GeV.

## 3 Discussion

Some measurements of jet and hadron production have been discussed. Calculations to NLO QCD do not give a complete description of all the measurements. With the large data samples now collected at LEP and HERA, these precise measurements could be used to constrain parton density functions and tune MC’s. Important theoretical work is needed in producing higher order calculations and/or NLO calculations with parton showers  and hadronisation. Both of these projects have started and are essential if theory is to keep up with the increasing precision of the data.


text:
# A Simple Model for the BFKL-DGLAP Transition

## 1 Introduction

When $`k_{}`$ is large and $`1/x`$ is limited we are in the DGLAP regime, and when $`1/x`$ is large and $`k_{}`$ is limited we are in the BFKL regime. An essential question is then: What is large? Where is the boundary between the regimes, and what is the behaviour in the transition region? We will here present a simple interpolating model, which can illuminate these questions. The model is relevant for small $`x`$, and has a smooth transition between large $`k_{}`$, where ordered chains dominate, and small $`k_{}`$, where non-ordered chains are most important. We will first discuss the case of a fixed coupling, and after that the case with a running coupling. The results are obtained in collaboration with Gabriela Miu, and a more extensive presentation is given in ref. 1.

## 2 Fixed coupling

For large values of $`k_{}`$ (in the “DGLAP region”) the non-integrated structure function $`(x,k_{}^2)`$ is dominated by contributions from chains which are ordered in $`k_{}`$. For small $`x`$ the result is a product of factors $`\frac{3\alpha _s}{\pi }\frac{dx_i}{x_i}\frac{dk_{,i}^2}{k_{,i}^2}`$ :

$`G{\displaystyle \frac{1}{\overline{\alpha }}}(x,k_{}^2){\displaystyle \underset{N}{}}{\displaystyle \stackrel{N}{}\overline{\alpha }dl_i\theta (l_{i+1}l_i)d\kappa _i\theta (\kappa _{i+1}\kappa _i)}`$
$`\mathrm{where}\overline{\alpha }{\displaystyle \frac{3\alpha _s}{\pi }},l\mathrm{ln}(1/x)\mathrm{and}\kappa ln(k_{}^2)`$ (1)
Integration over $`\kappa _i`$ with the $`\theta `$-functions gives the phase space for $`N`$ ordered values $`\kappa _i`$. The result is $`\kappa ^N/N!`$. The integrations over $`l_i`$ give a similar result, and thus we obtain the wellknown double log result

$$\frac{1}{\overline{\alpha }}(x,k_{}^2)\underset{N}{}\overline{\alpha }^N\frac{l^N}{N!}\frac{\kappa ^N}{N!}=I_0(2\sqrt{\overline{\alpha }l\kappa })$$
(2)
For smaller $`k_{}`$, in the BFKL region, also non-ordered chains contribute, and the result is a powerlike increase $`x^\lambda `$ for small $`x`$-values .

The CCFM  model interpolates smoothly between the DGLAP and BFKL regimes. The Linked Dipole Chain (LDC) model  is a reformulation and generalization of the CCFM model, in which gluons in some parts of phase space are treated as final state radiation, instead of as initial state radiation as in the CCFM model. In the LDC model the possibility to “go down” in $`k_{}`$, from $`\kappa ^{}`$ to a smaller value $`\kappa `$, is suppressed by a factor $`\mathrm{exp}(\kappa \kappa ^{})`$. The effective allowed distance for downward steps is therefore limited to about one unit in $`\kappa `$. If this quantity is called $`\delta `$, the result is that the phase space factor $`\frac{\kappa ^N}{N!}`$ is replaced by $`\frac{(\kappa +N\delta )^N}{N!}`$. Thus we obtain

$$\frac{1}{\overline{\alpha }}(x,k_{}^2)\underset{N}{}\frac{(\overline{\alpha }l)^N}{N!}\frac{(\kappa +N\delta )^N}{N!}$$
(3)
When $`\kappa `$ is very large, this expression approaches eq. (2). When $`\kappa `$ is small we find, using Sterling’s formula, that $`(x,k_{}^2)_N(\overline{\alpha }l\delta e)^N/N!=\mathrm{exp}(\lambda l)=x^\lambda `$, with $`\lambda =\overline{\alpha }\delta e`$. For $`\delta =1`$ this gives $`\lambda =e\overline{\alpha }=2.72\overline{\alpha }`$, which is very close to the leading order result for the BFKL equation, $`\lambda =4\mathrm{ln}2\overline{\alpha }=2.77\overline{\alpha }`$. Thus eq. (3) interpolates smoothly between the DGLAP and BFKL regimes. It is also possible to see that the transition occurs for a fixed ratio between $`\kappa `$ and $`l`$, $`\kappa /le\overline{\alpha }`$.

For small $`x`$ we can also go one step further, and expand eq. (3) in powers of $`\kappa /l`$. After some algebra we then find

$$\frac{1}{\overline{\alpha }}\frac{1}{\sqrt{2\pi \lambda l}}\mathrm{exp}\left[\lambda l+\frac{\kappa }{\delta }\frac{\kappa ^2}{2\delta ^2\lambda l}\right].$$
(4)
We here recognize the Gaussian distribution with a width proportional to $`\lambda l`$, corresponding to a random walk in $`\kappa `$, typical for the BFKL equation .

The result of a numerical evaluation of eq. (3) is presented in fig. 1a together with the DGLAP and BFKL approximations in eqs. (2) and (4). We note in particular that $``$ grows monotonically with $`\kappa `$. The Gaussian behaviour is only obtained in an expansion of $`\mathrm{ln}`$ to second order in powers of $`\kappa /l`$.

## 3 Running Coupling

For a running coupling a factor $`\overline{\alpha }d\kappa `$ in eq. (1) is replaced by $`\alpha _0d\kappa /\kappa =\alpha _0du`$, where $`\overline{\alpha }\alpha _0/\kappa `$ and $`u\mathrm{ln}\kappa =\mathrm{ln}\mathrm{ln}k_{}^2`$. In the large $`k_{}`$ region (the “DGLAP region”) the result is therefore similar to eq. (2), but with $`\overline{\alpha }\kappa `$ replaced by $`\alpha _0\mathrm{ln}\kappa `$. For small $`x`$-values we note that the allowed effective size of downward steps, which is still one unit in $`\kappa `$, is a sizeable interval in $`u`$ for small $`\kappa `$, but a very small interval in $`u`$ for larger $`\kappa `$. This is the reason for our earlier experience , which showed that a typical chain contains two parts. In the first part the $`k_{}`$-values are relatively small, and it is therefore easy to go up and down in $`k_{}`$, and non-ordered $`k_{}`$-values are important. The second part is an ordered (DGLAP-type) chain, where $`k_{}`$ increases towards its final value.

Let us study a chain with $`N`$ links, out of which $`Nk`$ correspond to the first part with small non-ordered $`k_{}`$, and the remaining $`k`$ links belong to the second part with increasing $`k_{}`$. Assume that the effective space for each $`u_i`$ in the soft part is given by a quantity $`\mathrm{\Delta }`$. ($`\mathrm{\Delta }`$ is sensitive to the soft $`k_{}`$-cutoff, and an estimate is presented in ref. 1.) Then the total weight for this part becomes $`\mathrm{\Delta }^{Nk}`$. For the $`k`$ links in the second, ordered, part the phase space becomes as in the DGLAP case $`u^k/k!`$. Thus the total result is ($`G/\overline{\alpha }=\kappa /\alpha _0`$)

$$G\underset{N}{}\frac{(\alpha _0l)^N}{N!}\underset{k=0}{\overset{N}{}}\frac{u^k}{k!}\mathrm{\Delta }^{Nk}=\underset{N}{}\frac{(\alpha _0l\mathrm{\Delta })^N}{N!}\underset{k=0}{\overset{N}{}}\frac{(u/\mathrm{\Delta })^k}{k!};u\mathrm{ln}\kappa .$$
(5)
This simple model also interpolates smoothly between the DGLAP and BFKL regions. For large $`u`$-values the last term in the sum over $`k`$ dominates, which gives the “DGLAP” result

$$G(\alpha _0lu)^N/(N!)^2=I_0(2\sqrt{\alpha _0lu})(16\pi ^2\alpha _0lu)^{1/4}\mathrm{exp}(2\sqrt{\alpha _0lu}).$$
(6)
For small $`u`$-values the sum over $`k`$ gives approximately $`\mathrm{exp}u/\mathrm{\Delta }`$, and thus

$$G\mathrm{exp}(\alpha _0\mathrm{\Delta }l)\mathrm{exp}(u/\mathrm{\Delta })=x^\lambda \kappa ^{\lambda /\alpha _0},\mathrm{with}\lambda =\alpha _0\mathrm{\Delta }.$$
(7)
The transition between the regimes occurs now for a fixed ratio between $`u=\mathrm{ln}\mathrm{ln}k_{}^2`$ and $`l=\mathrm{ln}(1/x)`$, $`u/l\lambda ^2/\alpha _0`$, which is of order 0.1 if $`\lambda `$ is around 0.3.

The result of a numerical evaluation is illustrated in fig. 1b. The DGLAP and BFKL limits in eqs. (6) and (7) are also presented, and we can see how the model interpolates between these limits.

## 4 Conclusion

The simple models in eqs. (3) and (5), for fixed and running couplings respectively, interpolate smoothly between large and small $`k_{}`$-values. They contain the essential features of the CCFM and LDC models, and can therefore give an intuitive picture of the transition between these two regimes. The transition occurs for a fixed ratio between $`\mathrm{ln}k_{}^2`$ and $`\mathrm{ln}1/x`$ for fixed coupling, and between $`\mathrm{ln}\mathrm{ln}k_{}^2`$ and $`\mathrm{ln}1/x`$ for a running coupling. A more extensive discussion, including effects of non-leading terms in $`\mathrm{ln}1/x`$ and the properties of Laplace transforms, can be found in ref. 1.


text:
# The Low-Energy Spectral Response of the ACIS CCDs

## 1. Introduction

The Advanced CCD Imaging Spectrometer (ACIS) instrument team, at both Penn State and MIT, and the Chandra X-ray Center (CXC) have been working to improve the quality of the flight calibration of the ACIS CCDs. The CXC has recently released a new response matrix for the S3 CCD at -120 C in the CALDB 2.7 and CIAO 2.1.3 release (see asc.harvard.edu/cal/Links/Acis/acis/Cal\_prods/
matrix/matrix.html). The new matrix includes a more sophisticated representation of the spectral redistribution function and an improved model of the gain at energies below 1 keV. The major difference between the CALDB 2.7 matrices and previous matrices is at energies below 1 keV, with little or no difference at energies above 1.5 keV. There has been no change to the quantum efficiency model of the detector.

The ACIS instrument team has released software (SW) to correct for the charge-transfer inefficiency (CTI) of both the Frontside-illuminated (FI) and Backside-illuminated (BI) CCDs. The physical model of the CCD is described in Townsley et al. (2001a) and the details of the CTI correction are given in Townsley et al. (2000) and Townsley et al. (2001b). The CTI correction SW (version 1.37) and accompanying response matrices (30JUL01) are available at “www.astro.psu.edu/users/townsley/cti/” or from the contributed SW page at the CXC web site “asc.harvard.edu/cgi-gen/cont-soft/soft-list.cgi”. The CTI correction for the BI CCD is small but is rather large for the FI CCDs, and makes a dramatic difference in spectral fits with the FI CCDs.

In this paper, we will focus on the low energy response as described by the response matrices listed above. Specifically, we will use strong lines in the spectrum of the bright, SMC supernova remnant (SNR) 1E0102.2-7219 in the region between 0.5-1.5 keV to test the new response matrices. The spectrum of 1E0102.2-7219 has been well-characterized using the gratings on XMM (see Rasmussen et al. 2001 and Sasaki et al. 2001) and also on Chandra (see Canizares et al. 2001, Flanagan et al. 2001, and Davis et al. 2001). The spectrum is dominated by strong lines from O VII, O VIII, Ne IX, Ne X, Mg XI, and Mg XII. The high resolution gratings spectra show that there is little or no Fe emission in this spectrum. The lack of Fe in the 1E0102.2-7219 spectrum greatly simplifies the modeling of the spectrum, since the Ne lines are not blended with the Fe L lines.

## 2. Data

We used six datasets for the S3(BI) analysis and five datasets for the I3(FI) analysis. The datasets are listed in Table 1 with the date of observation, position on the CCD, and exposure time. The six datasets for S3 sample four regions on the CCD, since two of the observations were executed at the same location on the chip but at different times to check for temporally-varying performance. The five datasets on I3 sampled five different locations in row number on the CCD in order to characterize the effect of CTI on low-energy X-rays. All of these observations are calibration observations and are available in the Chandra public archive.

The spectra were extracted and response matrices and auxiliary response files generated using the standard CIAO tools dmextract, mkrmf, mkarf using CIAO 2.1.3 and CALDB 2.7 for the S3 data. For the I3 data and also for the S3 dataset OBSID 1311, the spectra were extracted using dmextract but specifying the binning of PI channels described in the recipe on the CTI correction SW page. The different channel types are labeled in Figure 1 as “ PHA” and “PI” for the channel types produced by the CIAO SW and “PPI” for the channel types produced by the Penn State CTI-correction SW. We used the response matrices for I3 and S3 appropriate for CTI-corrected data based on the position of the source on the chip (see the CTI-correction web pages for more details). We used an annulus for the spatial extraction region in order to maximize the contribution of the shocked ejecta and to minimize the contribution of the primary blastwave of the remnant (see Hughes et al. 2001 and Gaetz et al. 2001). Extracting spectra from this annulus should maximize the contrast between the lines and the underlying continuum.

## 3. Spectral Analysis

### 3.1. Spectral Model

We developed a spectral model (based on the high-resolution gratings data) which consisted of two components for absorption (one for the galactic contribution and one for the SMC contribution with variable abundances), a bremsstrahlung component for the continuum, and Gaussians for the lines. The energies of the Gaussians were set to their known values and not allowed to vary during the fitting process. The widths of the lines were set to zero, so that the width is completely determined by the spectral redistribution function of the matrix. The HETG data indicate that there are measurable Doppler shifts in these lines on the order of a few eV, but these are small compared to the resolution of the detector ($`75150`$ eV). Only the normalizations of the lines were allowed to vary. In this manner, this model was used to verify that the lines were modeled at the correct energy in the matrices and the shapes of the lines were consistent with the widths contained in the matrices. Our model contains a total of 24 lines determined from the gratings spectra, but in most of the fits three or four of the lines are fitted with a zero normalization.

### 3.2. S3(BI) Spectral Fits

The spectral model described above was fit to all of the datasets in Table 1. The values of the reduced $`\chi ^2`$ ranged from 1.1 to 2.1. Table 2 and Figure 1 summarize the fitted values for the $`N_\mathrm{H}`$ and $`kT`$. OBSID 1311 was fitted in both PHA and PI channel space with the CALDB 2.7 matrices and the resulting values agreed well with each other. OBSID 1311 was also fitted after applying the CTI correction. The fitted value of $`kT`$ was consistent with the value derived from the CXC matrix fits, but the fitted value of the $`N_\mathrm{H}`$ was discrepant at the 90% confidence limit (CL). OBSIDs 1531, 1308, 1530, 141, and 1702 were fitted in PHA space and compared to OBSID 1311. All datasets agreed within the 90% CL except for OBSID 1530. OBSID 1308 agreed within the 90% CL, but had the largest error bars by far of any dataset. Both OBSID 1530 and 1308 warrant further investigation for instrumental effects which may have compromised the data.

Figure 2 displays the fit to OBSID 1311 in PHA space. The matrix and model represent the data extremely well from 700 eV to 1600 eV. Note how well the Ne lines at $`900`$ eV and $`1000`$ eV are fitted by this matrix and model. There are data from 1.6 keV to 5.0 keV which we have chosen to exclude from the plot since there are no strong lines in that region, but which are included in the fit. Only below 700 eV is there any evidence of systematic effects in the residuals. The fits to the other data sets and to the CTI-corrected data are qualitatively similar. These data confirm that the CALDB 2.7 matrix has correctly incorporated the gain and spectral redistribution function of the BI(S3) CCD and that data from four different locations on the S3 detector produce similar fitted results.

### 3.3. I3(FI) Spectral Fits

We fit the five datasets on the I3 CCD which span the regions from low row numbers where the CTI effects are minimized to high row numbers where the CTI effects are largest. The results are summarized for OBSIDs 440, 439, 136, 140, and 420 in Table 2 and Figure 1. The fitted values for the $`N_\mathrm{H}`$ and $`kT`$ are all consistent with each other at the 90% CL, except for the $`N_\mathrm{H}`$ for OBSIDs 420 and 439. The consistency of the fitted results demonstrates the value and accuracy of the CTI correction. The spectra for OBSIDs 420, 136, and 440 are shown in Figures 3, 4,and 5. The quality of the fits is good with the reduced $`\chi ^2`$ values ranging from 1.1 to 1.9. Only around 650 eV is there any evidence of systematic effects in the residuals. It is clear from the these plots that the spectral resolution of the detector changes dramatically with row number. The gain of the detectors is also changing, but this effect has also been corrected by the CTI correction SW. The data from OBSID 420 show that the spectral resolution of I3 at low row numbers is significantly better than S3 (see Virani et al. these proceedings). Note how cleanly the O and Ne lines are separated in Figure 3 and compare to the S3 spectrum in Figure 2. In fact, OBSID 136 shows that I3 has better spectral resolution than S3 halfway across the CCD. It is only in the top half of the FI CCDs where the resolution is worse than S3 after correcting for CTI as shown in Figure 5.

## 4. Conclusions

We have used the soft, line-dominated spectrum of the SMC SNR 1E0102.2-7219 to evaluate the response matrices of the BI and FI CCDs at low energies. We have demonstrated that the new BI(S3) matrices released in August 2001 in CALDB 2.7 and in the Penn State CTI correction SW are an accurate representation of the spectral response of the S3 CCD from 0.7 to 1.6 keV. We have used several datasets of 1E0102.2-7219 on S3 to verify that the new matrix yields consistent fitted parameters for different locations on the CCD. We have applied the Penn State CTI correction SW to five datasets on the I3(FI) CCD and demonstrated that the SW corrects for the spatial dependence of the gain and the spatial dependence of the redistribution function due to CTI. The fits indicate that the response matrix is an accurate representation of the spectral response of the I3 CCD. We derive fitted parameters which are consistent with each other at the 90% CL, even though the spectral resolution of the I3 CCD is changing rapidly with position. In fact, after CTI correction the FI CCDs have superior spectral resolution compared to S3 across the bottom half of the CCDs. Chandra observers should be aware of this performance and use it to their advantage if possible.

## 5. Acknowledgments

We thank all of our colleagues at the CXC, PSU, and MIT who have helped with the flight calibration of the ACIS instrument. We thank all of the engineers, technicians, and scientists who have made the Chandra X-ray Observatory such a success. PPP, RJE, and SNV acknowledge support for this work from NASA contract NAS8-39703. LKT and PSB acknowledge support for this work from NASA contract NAS8-38252.

## References

Canizares, C.R., Flanagan, K.A., Davis, D.S., Dewey, D. Houck, J.C., Schattenburg, M.L. 2001, in “Young Supernova Remnants” (11th Astrophysics Conference in Maryland), S. S. Holt & U. Hwang (eds), AIP, New York (2001)

Davis, D.S., Flanagan, K.A., Houck, J.C., Canizares, C.R., Allen, G.E., Schulz, N.S., Dewey, D., Schattenburg, M.L. 2001, in ”Young Supernova Remnants” (11th Astrophysics Conference in Maryland), S. S. Holt & U. Hwang (eds), AIP, New York (2001)

Flanagan, K.A., Canizares, C.R., Davis, D.S., Dewey, D. Houck, J.C., Markert, T.H., Schattenburg, M.L. 2001, in the ASP Conference Series for X-ray Astronomy 2000, Palermo Italy, eds. Giacconi, Stella & Serio

Gaetz, T. J., Butt, Y. M., Edgar, R. J., Eriksen, K. A., Plucinsky, P. P., Schlegel, E. M., & Smith, R. K. 2000, ApJL, 534, L47

Hughes, J. P., Rakowski, C. E., & Decourchelle, A. 2000, ApJL, 543, L61

Rasmussen, A. P., Behar, E., Kahn, S. M., den Herder, J. W., & van der Heyden, K. 2001, A&A, 365, L231

Sasaki, M., Stadlbauer, T. F. X., Haberl, F., Filipović, M. D., & Bennie, P. J. 2001, A&A, 365, L237

Townsley, L. K., Broos, P. S., Garmire, G. P., & Nousek, J. A. 2000, ApJL, 534, L139

Townsley, L. K., Broos, P. S., Chartas, G., Moskalenko, E., Nousek, J.A., & Pavlov, G.G., 2001, NIM, accepted, A

Townsley, L. K., Broos, P. S., Nousek, J.A., & Garmire, G.P., 2001, NIM, accepted, B

Virani, S.N., Plucinsky, P.P., Grant, C.E., & LaMarr, B. 2001, these proceedings


text:
# The production of proton and lepton fluxes in near Earth orbit

## 1 Introduction

Cosmic rays approaching the Earth interact with the atmosphere resulting in a substantial flux of secondary particles. The knowledge of composition, intensity and energy spectra of these particles is of considerable interest, e.g. for the evaluation of background radiation for satellites and the estimate of the atmospheric neutrino production for neutrino oscillation experiments .

The AMS measurements in near earth orbit  have allowed, for the first time, to gather accurate information on the intensity, energy spectra and geographical origin of charged particle fluxes at energies below the geomagnetic cutoff over a wide range of latitudes and at almost all longitudes. The under cutoff component of proton fluxes at equatorial latitudes has revealed an unexpected intensity of up to $`50\%`$ of the primary proton flux, a positron to electron flux ratio has been found in the undercutoff component which largely exceeds the cosmic one, differences in residence times and geographical origins have been reconstructed for positively and negatively charged particles.

A robust interpretation of these and many other characteristics of the undercutoff fluxes in terms of secondary particles produced in atmosphere requires an accurate description of both the interaction processes at their origin and of the geomagnetic field effects. Recently, different interpretations of the AMS measurements have been proposed  based on Monte Carlo simulations using different approaches on both the generation technique and the interaction model.

In this work, we report results from a fully 3D Monte Carlo simulation based on FLUKA 2000  for the description of cosmic ray interactions with the atmosphere. The key features of our analysis are an efficient generation technique for the incoming proton flux and a true microscopic, theory driven treatment of the interaction processes opposite to empirical parametrization of accelerator data. As a first attempt the contribution of $`He`$ and the heavier nuclei, representing ($`9\%`$)  of the all nuclei cosmic flux, is neglected.

In the following section we give a detailed description of the basic ingredients of this simulation, the generation technique and the interaction model. In section 3 we present the results on both protons and leptons and the comparison with AMS measurements. In section 4 we propose our conclusions.

## 2 The model

An isotropic flux of protons is uniformly generated on a geocentric spherical surface with a radius of 1.07 Earth radii ($`500Km`$ a.s.l.) in the kinetic energy range $`0.1170GeV`$.

We took the functional form suggested in  to describe the proton energy spectrum, the spectral index and the solar modulation parameter are extracted from a fit to the AMS data .

The magnetic field in the Earth’s proximity includes two components: the Earth’s magnetic field, calculated using a 10 harmonics IGRF  implementation, and the external magnetic field, calculated using the Tsyganenko Model222The external magnetic field is calculated only for distances greater than 2 Earth’s radii ($`E_R`$) from the Earth’s center . Its contribution to the total magnetic field is $`<1\%`$ at smaller distances and therefore can be safely neglected. . To account for the geomagnetic effects, for each primary proton we back-trace an antiproton of the same energy until one of the following conditions is satisfied:

1. the particle reaches the distance of $`10E_R`$ from the Earth’s center.
2. the particle touches again the production sphere.
3. neither 1 or 2 is satisfied before a time limit is reached.

If condition 1 is satisfied the particle is on an allowed trajectory, while if condition 2 is satisfied the particle is on a forbidden one. Condition 3 arises for only a small fraction of the events $`O(10^6)`$.

Particles on allowed trajectories are propagated forward and can reach the Earth’s atmosphere. The atmosphere around the Earth is simulated up to $`120Km`$ a.s.l. using 60 concentric layers of homogeneous density and chemical composition. Data on density and chemical composition are taken from the standard MSIS model . The Earth is modeled as a solid sphere which absorbs each particle reaching its surface.

### 2.1  The generation technique

The ideal approach in the generation of the primary cosmic rays spectra would be to start with an isotropic distribution of particles at a great distance (typically $`10E_R`$) from the Earth where the geomagnetic field introduces negligible distortions on the interstellar flux. However, this computational method is intrinsically inefficient since most of the particles are generated with trajectories which will not reach the Earth environment. Kinematic cuts can be applied in order to improve the selection efficiency at generation, however they tend to introduce a bias particles with low rigidity.

A good alternative to this approach is the backtracing method , adopted in the present analysis as outlined in the previous section. In the following, we will shortly discuss the validity of the technique and report the results of a comparison of the two methods. We recall that this method was applied for the first time in ref.  for the generation of atmospheric neutrino fluxes.

Let us consider first the effects of the geomagnetic field on an incoming flux of charged particles in the absence of a solid Earth. For the discussion, we start with an isotropic flux of monoenergetic333The realistic case of an energy spectrum can be treated just as a superposition of monenergetic cases protons at large distance, i.e. at infinity, from the origin of a geocentrical reference system. In this scenario, a negligible fraction of particles, with very particular initial kinematic parameters, will follow complicated paths and remains confined at a given distance from the origin (semi-bounded trajectories); for all practical purposes this sample can be ignored. Most of the particles will follow unbounded trajectories, reaching again infinity after being deflected by the magnetic field.

Unbounded trajectories cross a spherical surface centered in the field source only an even number of times, as shown in Fig.1: we call legs the trajectory parts connecting the spherical surface to infinity and loops the parts of the trajectory starting and ending inside the spherical surface.

Since each trajectory can be followed in both directions and no source or sink of particles is contained within the surface, the incoming and outgoing fluxes are the same. However, the presence of the magnetic field breaks the isotropy of the flux “near” the field source, so for a given location there is a flux dependence due to the direction.

Applying the Liouville Theorem, under the hypothesis of isotropy at infinity, it is straightforward to prove  that the proton flux in a random point is the same as at infinity along a set of directions (allowed directions), and zero along all the others (forbidden directions).

The pattern of the allowed and forbidden directions depends on both the rigidity and the location and is known as the geomagnetic cutoff.

With the introduction of a solid Earth, all the trajectories that are crossing the Earth are broken in two or more pieces (Fig.2): the legs become one-way trajectories and the loops disappear.

The presence of the Earth modifies the flux which exits from the surrounding spherical surface, since particles are absorbed by the Earth, while it has only a minimal effect on the incoming flux which is modified only by the absence of certain loops. To generate the flux of particles reaching the Earth’s atmosphere, it is sufficient to follow the particles along the allowed trajectories corresponding to the legs, taking care to avoid double or multiple counting.

To respect this prescription we reject all trajectories that are back-traced to the production sphere, this allow us to correctly consider the cases like the one shown in Fig.3.

We point out that an important difference with respect to the application in the neutrino flux calculation of  is that for the former, the generation sphere coincided with the Earth’s surface, and therefore the forbidden trajectories included those which touched again the Earth (plus those who remained trapped for a long time). In that case there are no problems of double counting.

To check the validity of our technique we made a test comparing the results of the inefficient generation technique at 10 Earth’s radii distance from the Earth’s center with the backtracing technique described in this paper.

Fig. 4 shows this comparison for several cha-racteristic distributions, the agreement between the two methods is good.

### 2.2  The interaction model

We use the software package FLUKA 2000  to transport the particles and describe their interactions with Earth’s atmosphere. The setup of this simulation is derived from the one used in . This package contains a tridimensional description of both electro-magnetic and hadronic interactions. This code is benchmarked against a wide set of data and is already used in many applications, ranging from low energy nuclear physics to high energy accelerator and cosmic ray physics. For this reason we have preferred this model with respect to the use of “ad hoc” parametrizations of particle production in the energy range of our interest.

In FLUKA hadronic interactions are treated in a theory-driven approach, and the models and their implementations are guided and checked using experimental data. Hadron-nucleon interaction models are based on resonance production and decay below an energy of few $`GeV`$ and on the Dual Parton Model above. The extension from hadron to hadron-nucleus interactions is done in the framework of a generalized intra-nuclear cascade approach including the Gribov-Glauber multi-collision mechanism for higher energies followed by equilibrium processes: evaporation, fission, Fermi break-up and $`\gamma `$ de-excitation. The parameters of the models embedded in the FLUKA package are fixed only by comparing expectations with data from accelerator experiments.

In fig 5 a) we show the map of the primary proton interaction points in geographical coordinates. The distribution reflects the influence of the geomagnetic cutoff. Fig 5 b) shows the interactions altitude profile, the solid histogram is for $`E_k<30GeV`$ while the dashed one is for $`E_k>30GeV`$. The mean interaction altitude depends weakly on the energy.

The cosmic proton impinging in the atmosphere are doing elastic scattering in the 24% of the events and inelastic interactions in the remaining 76%, in tab.1 we show some characteristic of the inelastic interactions as simulated by FLUKA 2000.

## 3 Comparison with the AMS data

To compare with the AMS data, we define a detection boundary corresponding to a spherical surface at the AMS orbit altitude ($`400\mathrm{Km}`$ a.s.l). We record each particle that crosses the detection boundary within the AMS field-of-view, defined as a cone with a 32$`^\mathrm{o}`$ aperture with respect to the local zenith or nadir directions.

To obtain the absolute normalization, we take into account the field-of-view, the corresponding AMS acceptance, and an Equivalent Time Exposure (E.T.E.) corresponding to the number of the generated primary protons.

Our results are based on a sample of $`610^6`$ primary protons generated in the kinetic energy range of $`0.1170GeV`$, which corresponds to $`\mathrm{4\; 10}^{12}s`$ (E.T.E).

### 3.1 Protons

In Fig.6, we show the comparison between the fluxes obtained with the simulation and the measured AMS downgoing proton flux  in nine bins of geomagnetic latitude ($`\theta _M`$) . Fig.7 shows the same comparison for the upgoing proton flux in four selected bins of $`\theta _M`$.

As seen in Fig.6, the simulation well reproduces at all latitudes the high energy part of the spectrum and the falloff in the primary spectrum due to the geomagnetic cutoff, thus validating the general approach used for the generation and detection, as well as the tracing technique.

A good agreement among data and simulation is also found in the under-cutoff part of the spectra. The small and systematic deficit which can be seen in the secondary component of the simulated fluxes is of the same order of the expected contribution from the interaction of cosmic He and heavier nuclei.

This flux is due to the secondaries produced in the atmosphere and that spiral along the geomagnetic field lines up to the detection altitude. Therefore it is sensitive to specific aspects of the interaction model and to the accuracy of the particle transport algorithm.

A correct quantitative prediction of this part of the spectra depends on the quality of low energy nucleon production both in terms of yield and energy distribution. This is in part due to the fragmentation of the target nucleus, and depends on the details of the nuclear physics algorithms describing excitation and break up.

From the analysis of the motion of the secondary protons from their production up to their detection, it can be pointed out that a fraction of the observed flux is due to a multiple counting of the same particles. Within the formalism of adiabatic invariants , it is seen that charged particles trapped in the geomagnetic field, i.e. the undercutoff protons, move along drift shells which can be associated with a characteristic residence time444The mean time after which a particle is absorbed into the atmosphere. In our case it represents the effective life time of the particle.  that depends on the fraction of the shell located inside the Earth’s atmosphere. Thus, particles moving along long-lived shells have a large probability to cross many times a geocentered spherical detector, while those moving along short-lived shells typically cross the detector only once.

The drift shells crossing the AMS orbit, at an altitude of 400 km, are in general short-lived, however in the equatorial region the long-lived shells are present as well .

In the following, we will indicate as the real proton flux that one obtained by counting only once each particle crossing the detector: its intensity is indicated by the dashed distributions in figs. 6 and 7.

A quite relevant effect can be seen in the equatorial region: there the AMS measurement indicates an important secondary proton flux while the real number of protons crossing the detector is more than one order of magnitude lower. At high geomagnetic latitudes, the solid and dashed lines tend to merge. The effect becomes negligible for $`\theta _M>0.8`$.

This can be better seen in Fig.8, where the integral primary proton flux seen by AMS is shown as a function of geomagnetic latitude. The intensities of the real and measured undercutoff fluxes are reported in the same plot for comparison and their ratio with the primary component shown in Fig.9. A minor contribution from the undercutoff proton component can be therefore expected in the atmospheric shower development and neutrino production.

In Fig.10, the residence time is plotted versus the kinetic energy of the trapped secondary protons for $`|\theta _M|<0.4`$. In the scatter plot it is possible to distinguish the populations corresponding to long-lived and short-lived shells similar to those shown in  for leptons.

Fig.11 shows the distribution of trapped secondary proton end points for $`|\theta _M|<0.4`$, Fig.11a is for a lifetimes smaller than 0.3 s., while Fig.11b is for a lifetimes greater than 0.3 s.. The end point distribution agrees with the location of the intersections of the drift shells with the atmosphere as experimentally verified by , and discussed in .

### 3.2 Electrons and positrons

In Fig. 12 we show a comparison of the simulated undercutoff electron and positron downgoing fluxes with the corresponding AMS measured fluxes .

We remind that the AMS positron measurement is restricted at energies below few GeVs, with a dependence of the maximum energy on the geomagnetic latitude which reflects the increasing proton background with $`\theta _m`$.

A comparison of data and simulation in the high energy part of the electron spectra is not possible, since the cosmic electrons have not been used as an input in the current work. However, their contribution to the cosmic rays reaching the atmosphere is $`O(10^2)`$ leading to a negligible effect of in the generation of the undercutoff fluxes.

The simulation well reproduces the general behavior of the undercutoff part of the spectra in terms of shape and intensity; a similar agreement is observed for the upgoing lepton spectra (not shown). The real lepton fluxes, corresponding to the real proton flux described earlier, are shown with the dashed line distribution in Fig. 12. As in the case of protons, a large effect from multiple crossing is present going toward the equatorial region, more pronounced for the positron component.

As for the undercutoff protons, we would have expected a systematic deficit in the simulated electron fluxes coming from the missing contribution of helium and heavier nuclei to the CR fluxes. Subcutoff $`e^\pm `$ are mainly (97$`\%`$) coming from decays of pions produced in the proton collisions with the atmospheric nuclei: charged pions contribute through the $`\pi \mu e`$ chain, while $`\pi ^0`$ through $`\pi ^0\gamma \gamma `$ with subsequent e.m. showers. The relative contribution of charged pions to the subcutoff electrons (positron) fluxes at AMS altitude in our simulation is found to be $`37\%`$ ($`47\%`$), while the remaining $`60\%`$ ($`50\%`$) appears to come from $`\pi ^0`$ production. This point deserves some considerations. The level of agreement in the comparison of data and predictions for $`e^\pm `$ fluxes turns out to be an important benchmark for the interaction model in view of a discussion on particle production in atmosphere, since is strictly linked to the meson production (mostly pions at this energy). This work complements other studies oriented to the validation of the FLUKA model in terms of particle yields. In particular, the quality of $`\pi ^\pm `$ generation in our interaction model has been already checked in  through the comparison with muon fluxes measurements at different depths in atmosphere. The muons are from the charged pions decay chain and experimental data  are well reproduced by the simulation. In the case of $`e^\pm `$ also $`\pi ^0`$’s become relevant. Usually, when parametrized interaction models are used, like in the works of ref. , the yield of $`\pi ^0`$ is fixed assuming a priori an exact charge symmetry in pion production. In practice, this is also made necessary by the large errors that affect the scarce existing experimental data on neutral pion production. Instead, in the case of a microscopic interaction model like FLUKA, there is no constraint of this type, and the balance of $`\pi ^0`$ vs. $`\pi ^\pm `$ automatically emerges from the feature of the model. Recently, it has been pointed out how in FLUKA there exists a significant violation of charge symmetry: $`\pi ^0`$’s are in general more ($``$20%) of the average of $`\pi ^+`$ and $`\pi ^{}`$. Technically this symmetry violation emerges in the hadronization of color strings and this normally occurs also in other codes, like JETSET or PITHYA. We cannot enter here in a detailed discussion of this point, and we limit ourselves to say that there are reasons to believe that, at least for laboratory energies below 100 GeV, the acceptable value of charge asymmetry should be lower than that resulting from the present version of the code. However, the comparisons of predictions to data discussed in this work, combined with the mentioned work on atmospheric muons, already tell that the predicted fraction of $`\pi ^0`$ cannot be significantly wrong, although no definitive quantitative conclusion can be extracted, since the nuclear component has not been yet introduced in the primary spectrum.

In Fig.13b-c we show the integrated positron and electron downgoing fluxes for the kinetic energy range $`0.21.5GeV`$ as a function of $`\theta _M`$. Their ratio is shown in Fig.13a. One of the most remarkable features of the AMS measurement is the large value of this ratio, when compared to the natural cosmic value, and its latitude dependence. In Fig. 12, the contribution from primary protons with $`E_k>30GeV`$ to the electron and positron fluxes is illustrated by the filled area. We can notice that in the equatorial region, the electrons are produced essentially by primary protons with $`E_k>30GeV`$ , while for the positrons lower energy protons contribute as well. This distinction disappears at higher latitudes, where positron and electrons are produced by the protons in the same energy range.

This behavior reflects the East-West asymmetry of the geomagnetic cutoff on primary protons and the larger probabilities of escape from atmosphere for secondary electrons(positrons) generated by Westward(Eastward) moving protons .

Positrons are preferentially injected on drift shells reaching the AMS altitude by eastward moving protons, which experience a lower rigidity cutoff than westward moving ones. This mechanism is more effective at the equator, where the cutoff is larger and its asymmetry maximal, resulting in the excess of undercutoff positrons from low energy protons as indicated by our simulation. The cutoff mechanism becomes irrelevant at high latitudes, where any difference in positron and electron production should be given instead by different $`\pi ^+/\pi ^{}`$ production. Nor the data neither our simulation indicate, within their uncertainties, a relevant charge asymmetry from this source.

## 4 Conclusions

The interactions of cosmic ray protons with the Earth’s atmosphere have been investigated by means of a fully 3D Monte Carlo program.

The proton, electron and positron undercutoff flux intensities measured by AMS, as well as their energy spectra, have been correctly reproduced by our simulation

Geomagnetic effects, and in particular the east-west asymmetry in the cosmic protons rigidity cutoff, have been confirmed as the mechanism responsible for the measured excess of the positron component.

The main features of the geographical origin and residence time distributions for both protons and leptons have been replicated and the effect of multiple crossing of the detector by spiraling secondaries in the geomagnetic field briefly discussed.

Our results indicate that the intensity of the undercutoff proton flux, when the multiple counting is taken into account, never exceeds a 10% of the cosmic proton flux, representing a negligible source for atmospheric production of secondaries. However, this aspect will be object of further and more refined study in the future.

The analysis on the possible strategies to generate the cosmic rays incoming flux has shown the validity of a backtracing approach as an accurate and highly efficient technique.

This work provides also additional way to validate the features of the adopted particle production model. In particular, the study of $`e^\pm `$ fluxes has revealed to be to be an interesting instrument to check the meson production in primary interactions, and the results are satisfactory. We believe that our simulation, validated by the high statistic measurements of AMS, can be used to assess the radiation environment in near Earth orbit, and represents a valuable tool for a more accurate calculation of particle fluxes in atmosphere.

This work has been partially supported by the Italian Space Agency (ASI) under contract ARS-98/47.