format: text

text:
# Non-adiabatic geometric quantum computation with trapped ions

## Abstract

We propose a non-adiabatic scheme for geometric quantum computation with trapped ions. By making use of the Aharonov-Anandan phase, the proposed scheme not only preserves the globally geometric nature in quantum computation, but also provides the advantage of non-adiabaticity that overcomes the problem of slow evolution in the existing adiabatic schemes. Moreover, the present scheme requires only two atomic levels in each ion, making it an appealing candidate for quantum computation.

Conventionally, the controllable operations in quantum computation (QC) are achieved on the basis of dynamical origins of quantum state evolutions Div95 . In recent years, for the purpose of being fault-tolerant to certain types of computational errors, there are considerable interest in exploiting the possibility of implementing quantum computation by geometrical means, which have been termed as holonomic quantum computation Zan99 ; Pac00 ; Ell01 . Depending on the degenerate property of the eigenspace of the governing Hamiltonian, the holonomy can be either a simple Abelian Berry phase factor Ber84  or a general non-Abelian unitary transformation Wil84 . It has been shown that the universal quantum computation can be implemented in principle by holonomies Zan99 ; Pac00 ; Ell01 . Up to date, several experimental proposals have been suggested for geometric quantum computation, using such as the nuclear magnetic resonance Jon00 , superconducting nanocircuits Fal00 ; Cho01 , trapped ions Dua01 ; Pac01 , and nonlinear optics Pac00a .

The principle of the aforementioned geometric QC is rooted in the adiabatic evolution of quantum system, which may thus imply a slow computing speed. The adiabatically slow evolution may also challenge the sustainment of required coherence in QC. Therefore, geometric QC based on non-adiabatic evolution should be desirable. Very recently, Wang and Keiji suggested to exploit a non-adiabatic evolution to realize geometric QC in NMR system and superconductor nanocircuits Wan01 . Indeed, geometric phase exists in non-adiabatic evolving quantum systems, which is in fact the Aharonov-Anandan (A-A) phase Aha87 . Strictly speaking, the A-A phase depends on certain dynamical quantities such as the rotating angular speed of external (magnetic) field or state vector Bul88 ; Fev92 ; Ni95 . In this sense, the A-A phase differs from the adiabatic Berry phase. However, the dependence of the A-A phase on the angular speed is through the closed path loop depicted by the ending point of state vector, and hence is global in nature that largely retains the geometric sense of A-A phase. Accordingly, quantum computation based on the non-adiabatic A-A phase can be reasonably regarded as a kind of geometric QC.

In this work we propose a scheme for non-adiabatic geometric quantum computation with trapped ions. Besides removing the drawback of the slow adiabatic evolution, the proposed non-adiabatic scheme holds additional merits as follows. Firstly, there is no need to design the reverse evolving path to eliminate dynamical phases that occur in the Berry phase-based, adiabatic geometric QC operations in nondegenerated systems Jon00 ; Fal00 . Secondly, in comparison to the existing fully holonomic QC schemes Cho01 ; Dua01 , the present one does not involve the complicated construction of the degenerate eigenspace of driving Hamiltonian. The above two merits stem from the following observations. In a non-adiabatic quantum evolution, the geometric A-A phase is in general accompanied by a dynamical phase. However, if the evolving path is such designed that along it the state vector is always perpendicular to the driving (magnetic) field, the resulting phase factor after a non-trivial cyclic evolution will be purely geometric. This feature has been exploited by Suter et al in their seminal experiment for demonstrating the A-A phase Sut88 . Finally, only two atomic levels of each ion are needed in our scheme. This merit alone is attractive, since the originally proposed ion-trap QC scheme required three levels Cir95  and the recently proposed holonomic ion-trap QC required four levels Dua01 ; Pac01 . We notice that valuable efforts on improving ion-trap QC protocol by using only two levels have been carried out in dynamic schemes Mon97 ; Chi01 ; Mol99 . In particular, the technique proposed in Ref. Mol99, , which effectively couples the electronic states of a pair of ions by virtually exchanging phonons, is similar to our present one. The major contribution of our work is to perform QC by geometric means. Viewing that the work of Ref. Mol99,  has in fact extended the ion-trap QC scheme to finite temperature, our geometric scheme may also hold to similar regime, although in the following we would restrict our description at zero temperature limit.

Model Description. For quantum logic with trapped ions, we assume that each ion has two relevant internal states $`|0`$ and $`|1`$ with energy separation $`\omega _0`$, and, as usual Cir95 , can be selectively addressed by lasers. Consider, for instance, the $`j`$th ion being exposed to a traveling-wave laser field $`𝐄(𝐳)=𝐄_0\text{cos}(𝐤𝐳\omega _Lt+\varphi )`$ with frequency $`\omega _L`$, wave vector $`𝐤`$, and phase $`\varphi `$. Here $`𝐳=z_0\widehat{𝐳}(a+a^{})`$ is the center of mass coordinate of the ion in terms of the phonon raising (lowering) operator $`a^{}(a)`$ and zero-point spread $`z_0(\mathrm{}/2M\omega )^{1/2}`$, where $`M`$ is the total mass of the ion chain, and $`\omega `$ the phonon frequency. The resulting Hamiltonian reads

$`^{(j)}={\displaystyle \frac{\omega _0}{2}}\sigma _j^z+\stackrel{~}{\omega }_j\left[\sigma _j^+e^{i\eta (a+a^{})i\omega _Lt+i\varphi }+\text{H.c.}\right].`$ (1)
Here, the atomic Pauli operators $`\sigma ^z=|11||00|`$, $`\sigma ^+=|10|`$, and $`\sigma ^{}=|01|`$, are introduced. In Eq. (1) $`\stackrel{~}{\omega }_j`$ is the Rabi frequency, and $`\eta =(𝐤\widehat{𝐳})z_0`$ is the Lamb-Dicke parameter, which accounts for the coupling strength between internal and motional states. In our geometric scheme, the phonon plays no roles in one-bit operation, but it does involve in two-bit logic gate implementation as it effectively couples two ions together.

One-Bit Rotation via Non-adiabatic A-A Phase. As just mentioned, in single bit logic operation, we do not need the phonon-assisted dynamical flipping. This can be the case under situations such as (i) at zero temperature (no phonon excitation), if $`\omega _L<\omega _0`$; and (ii) at finite temperature, but $`\omega _0\omega _Ln\omega `$ ($`n`$ is an integer). In the absence of phonon participation, Eq. (1) can be re-expressed as

$`_R^{(j)}=𝛀_j𝝈_j,`$ (2)
in the rotating frame (with angular velocity $`\omega _L\widehat{𝐞}_z`$). Here notations have been introduced for the effective magnetic field $`𝛀_j\{\stackrel{~}{\omega }_j\text{cos}\varphi ,\stackrel{~}{\omega }_j\text{sin}\varphi ,(\omega _0\omega _L)/2\}`$, and the vector Pauli operator $`𝝈_j\{\sigma _j^x,\sigma _j^y,\sigma _j^z\}`$.

To show how rotation between logic states $`|0`$ and $`|1`$ can be performed by geometric means, we first consider the geometric evolution of the eigenstates of $`\sigma ^y`$, defined by $`\sigma ^y|\pm =\pm |\pm `$, under appropriately designed laser sequences. Specifically, the state $`|+`$, for instance, will complete a cyclic evolution and acquire a non-adiabatic A-A phase under the following operations: (i) Switching on a $`\pi `$-pulse with laser phase $`\varphi =0`$, the state $`|+`$ rotates around an effective magnetic field $`𝛀_1=\{\stackrel{~}{\omega }_j,0,(\omega _0\omega _L)/2\}`$, from $`|+`$ in the $`\widehat{𝐞}_y`$ direction to $`|`$ in the $`\widehat{𝐞}_y`$ direction along the curve ACB on the Bloch sphere; see Fig. 1. (ii) Suddenly changing the laser phase to $`\varphi =\pi `$, after another $`\pi `$-pulse, the state $`|`$ rotates back to $`|+`$ around $`𝛀_2=\{\stackrel{~}{\omega }_j,0,(\omega _0\omega _L)/2\}`$ along curve BDA on the Bloch sphere. According to the A-A phase theory, after the above cyclic evolution, the state $`|+`$ will acquire a geometric phase $`e^{i\gamma }`$, with $`\gamma =4\mathrm{arctan}[2\stackrel{~}{\omega }_j/(\omega _0\omega _L)]`$. Note that during the above operation, the state vector keeps always perpendicular to the effective magnetic field, thus no dynamical phase is accumulated in the evolution. Similarly, the state $`|`$ will acquire A-A phase $`e^{i\gamma }`$ at the same time.

Now consider the evolution of logic states $`|0=\frac{i}{\sqrt{2}}(|+|)`$, and $`|1=\frac{1}{\sqrt{2}}(|++|)`$. After the above operations, they evolve to states

$`|0`$ $``$ $`\mathrm{cos}\gamma |0+\mathrm{sin}\gamma |1,`$
$`|1`$ $``$ $`\mathrm{cos}\gamma |1\mathrm{sin}\gamma |0.`$ (3)
We see here the geometric A-A phase plays a role of rotating the logic states. Particularly, complete state flipping between $`|0`$ and $`|1`$ can take place at $`\gamma =\pi /2`$. Note also that the possible value of $`\gamma `$ ranges from $`0`$ to $`2\pi `$, implying the ability of arbitrary rotation between $`|0`$ and $`|1`$.

The state evolution of the performed qubit described by Eq. (Non-adiabatic geometric quantum computation with trapped ions) is expressed in the rotating frame with frequency $`\omega _L`$, in which other free (not performed) qubits would have relative phase accumulations in the non-resonant case of $`\omega _L\omega _0`$. Conventionally, a more convenient choice is to express states in the interaction picture with respect to $`_0=\frac{\omega _0}{2}_j\sigma _j^z`$ (equivalently, a rotating frame with frequency $`\omega _0`$ around $`z`$-axis). Accordingly, Eq. (Non-adiabatic geometric quantum computation with trapped ions) can be recast in the interaction picture as

$`|0e^{i\omega _D\tau /2}\mathrm{cos}\gamma |0+e^{i\omega _D\tau /2}\mathrm{sin}\gamma |1,`$
$`|1e^{i\omega _D\tau /2}\mathrm{cos}\gamma |1e^{i\omega _D\tau /2}\mathrm{sin}\gamma |0,`$ (4)
where $`\omega _D=\omega _0\omega _L`$ and $`\tau `$ is the total operation time on the performed qubit.

Interestingly, the above state rotation (flipping) is performed by non-resonant pulses via geometrical means. Now we show that by resonant pulses, a pure phase shift gate of single qubit can be geometrically realized. In resonant case, the laser-frequency-associated rotating frame coincides with the interaction picture defined by $`_0`$, in which the ($`j`$th) qubit Hamiltonian reads $`^{(j)}=\stackrel{~}{\omega }_j\left[\sigma _j^+e^{i\varphi }+\text{H.c.}\right]=𝛀_j𝝈_j`$, with $`𝛀_j=\{\stackrel{~}{\omega }_j\text{cos}\varphi ,\stackrel{~}{\omega }_j\text{sin}\varphi ,0\}`$. We see that in the rotating frame the effective magnetic field constantly has zero $`z`$-axis component for arbitrary laser phase $`\varphi `$, i.e., it lies in the $`x`$$`y`$ plane. To realize the single bit phase gate, we first switch on a $`\pi `$-pulse with laser phase at certain value, say, $`\varphi _0`$. The logic state $`|0`$ and $`|1`$ would rotate around the effective magnetic field $`\{\stackrel{~}{\omega }_j\text{cos}\varphi _0,\stackrel{~}{\omega }_j\text{sin}\varphi _0,0\}`$ to $`|1`$ and $`|0`$, respectively. Then, suddenly changing the laser phase to $`\varphi _0`$, after another $`\pi `$-pulse, the logic states rotate around the effective magnetic field $`\{\stackrel{~}{\omega }_j\text{cos}\varphi _0,\stackrel{~}{\omega }_j\text{sin}\varphi _0,0\}`$, and return to the original states $`|0`$ and $`|1`$. Associated with this two step cyclic evolution, the logic states will respectively acquire geometric A-A phases as

$`|0`$ $``$ $`e^{i\stackrel{~}{\gamma }}|0,`$
$`|1`$ $``$ $`e^{i\stackrel{~}{\gamma }}|1,`$ (5)
where $`\stackrel{~}{\gamma }=4\varphi _0`$. With the help of this phase shift operation, the additional phase factor in Eq. (Non-adiabatic geometric quantum computation with trapped ions) can be canceled out by properly choosing $`\varphi _0`$. More importantly, together with this phase shift gate, the qubit state rotation Eq. (Non-adiabatic geometric quantum computation with trapped ions) constitutes a complete logic set for arbitrary single qubit operation.

Two-Bit Gate via Non-adiabatic A-A Phase. To complete the universal gate for quantum computation, non-trivial two-bit gate such as the controlled NOT, or equivalently, the conditional phase shift (CPS) gate, would be required. Below we show how the CPS gate can be implemented via geometric means. Consider two qubits (e.g. the $`j`$th and $`k`$th ones) irradiated by two lasers with frequencies $`\omega _{L,1}`$ and $`\omega _{L,2}`$, and phases $`\varphi _1`$ and $`\varphi _2`$. By setting $`\omega _{L,1}>\omega _0`$ and $`\omega _{L,2}<\omega _0`$, and correspondingly denoting the detunings by $`\delta _1=\omega _{L,1}\omega _0`$ and $`\delta _2=\omega _0\omega _{L,2}`$, the effective coupling between the two-bit states $`|00`$ and $`|11`$ can be established via virtually exchanging phonons, and the resulting two-bit effective Hamiltonian reads Mol99

$`\stackrel{~}{}^{(j,k)}`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{4}{}}}E_m|\stackrel{~}{m}\stackrel{~}{m}|+g_{jk}[e^{i(\omega _{L,1}+\omega _{L,2})t}`$ (6)
$`\times e^{i(\varphi _1+\varphi _2)}\sigma _j^+\sigma _k^++\text{H.c.}].`$
Here notations $`\{|\stackrel{~}{1}=|11,|\stackrel{~}{2}=|00,|\stackrel{~}{3}=|10,|\stackrel{~}{4}=|01\}`$ are introduced for the two-bit computational basis states. Up to the first-order expansion of the Lamb-Dicke parameter $`\eta `$ in Eq. (1), the effective coupling strength can be obtained via second-order perturbation theory as $`g_{jk}=g_jg_k(\frac{1}{\delta _1\omega }\frac{1}{\delta _2+\omega })`$, where $`g_{j(k)}=\stackrel{~}{\omega }_{j(k)}\eta `$ is the one-phonon involved Rabi frequency of single ion transition. The four basis-state energies $`E_m`$ ($`m=1,\mathrm{},4`$) contain also the ac Stark shifts. The effective interaction couples only between $`|\stackrel{~}{1}`$ and $`|\stackrel{~}{2}`$, but leaves $`|\stackrel{~}{3}`$ and $`|\stackrel{~}{4}`$ inactive with respect to the laser operation in study. This effective two state dynamics (e.g. Rabi oscillations) has been demonstrated in Ref. Mol99,  by numerical simulation starting from the original Hamiltonian. As a result, for two-bit operation, one can focus on state evolution in the subspace $`\{|\stackrel{~}{1},|\stackrel{~}{2}\}`$ in which the two-bit Hamiltonian (6) is recast to a similar form as the single qubit Note1

$`\stackrel{~}{}_R^{(j,k)}={\displaystyle \frac{\stackrel{~}{\omega }_D}{2}}\mathrm{\Sigma }^z+g_{jk}\left(e^{i\mathrm{\Phi }}\mathrm{\Sigma }^++e^{i\mathrm{\Phi }}\mathrm{\Sigma }^{}\right),`$ (7)
where $`\mathrm{\Phi }=\varphi _1+\varphi _2`$, and the two-bit Pauli matrices are introduced as $`\mathrm{\Sigma }^z|\stackrel{~}{1}(\stackrel{~}{2})=\pm |\stackrel{~}{1}(\stackrel{~}{2})`$, $`\mathrm{\Sigma }^+|\stackrel{~}{2}=|\stackrel{~}{1}`$, and $`\mathrm{\Sigma }^{}|\stackrel{~}{1}=|\stackrel{~}{2}`$. Note also that the Hamiltonian (7) has been expressed in the rotating frame with respect to $`\mathrm{\Sigma }^z`$ with the rotation frequency $`\omega _{L,1}+\omega _{L,2}`$; thus $`\stackrel{~}{\omega }_D=(E_1E_2)(\omega _{L,1}+\omega _{L,2})`$. Simple comparison of Eq. (7) with Eq. (2) indicates that an arbitrary rotation between $`|\stackrel{~}{1}`$ and $`|\stackrel{~}{2}`$ can be performed geometrically as that in the single bit case.

In the spirit of pulse-sequence operations in dynamic scheme based on an XY spin model Los98 ; Ima99 ; Li02 , we shall in the following show that the CPS gate can be implemented geometrically as

$`𝒰_{\text{CPS}}=`$ $`e^{i\pi /4}e^{i\pi 𝐧_j𝝈_j/3}e^{i\pi 𝐧_k𝝈_k/3}e^{i\pi \sigma _k^x/2}`$ (8)
$`\times U_{jk}(\pi /4)e^{i\pi \sigma _j^y/2}U_{jk}(\pi /4)e^{i\pi \sigma _j^x/2}.`$
Here the vector Pauli operator $`𝝈=(\sigma ^x,\sigma ^y,\sigma ^z)`$, unit vector $`𝐧_j=(1,1,1)/\sqrt{3}`$, and $`𝐧_k=(1,1,1)/\sqrt{3}`$. In the interaction picture with respect to the original free ion Hamiltonian, the two-bit rotation operator $`U_{jk}(\alpha )`$ is defined via $`U_{jk}(\alpha )|01=|01`$, $`U_{jk}(\alpha )|10=|10`$, and

$`U_{jk}(\alpha )|00`$ $`=`$ $`\text{cos}\left({\displaystyle \frac{\alpha }{2}}\right)|00+i\text{sin}\left({\displaystyle \frac{\alpha }{2}}\right)|11,`$
$`U_{jk}(\alpha )|11`$ $`=`$ $`\text{cos}\left({\displaystyle \frac{\alpha }{2}}\right)|11+i\text{sin}\left({\displaystyle \frac{\alpha }{2}}\right)|00.`$ (9)
Since the arbitrary rotation and phase shift gate of single qubit have been built up, the single bit operations in Eq. (8) can be readily implemented via geometric means by properly combining the single bit logic elements. Viewing the similarity between the two-bit reduced Hamiltonian (7) in the subspace $`\{|00,|11\}`$ and the single bit Hamiltonian (2), the two-bit rotation of Eq. (Non-adiabatic geometric quantum computation with trapped ions) can be straightforwardly implemented by the following two-step procedures:

(i) In the two-bit rotating frame with frequency $`\omega _{L,1}+\omega _{L,2}`$ around $`\mathrm{\Sigma }^z`$, performing cyclic evolution for the eigenstates of $`\mathrm{\Sigma }^y`$ by controlling the laser phases $`\mathrm{\Phi }`$ similarly as in the single bit case, one can geometrically rotate the states $`|00`$ and $`|11`$. Expressed in the interaction picture which also corresponds to $`\stackrel{~}{}_0E_1|1111|+E_2|0000|`$, this operation can realize the following state transformation

$`|00`$ $``$ $`e^{i\stackrel{~}{\omega }_D\tau /2}\mathrm{cos}\mathrm{\Gamma }|00+e^{i\stackrel{~}{\omega }_D\tau /2}\mathrm{sin}\mathrm{\Gamma }|11,`$
$`|11`$ $``$ $`e^{i\stackrel{~}{\omega }_D\tau /2}\mathrm{cos}\mathrm{\Gamma }|11e^{i\stackrel{~}{\omega }_D\tau /2}\mathrm{sin}\mathrm{\Gamma }|00.`$ (10)
Here $`\mathrm{\Gamma }`$ is the geometric A-A phase determined by the evolution contour of the two bit state vector, and $`\stackrel{~}{\omega }_D\tau `$ is the detuning-induced phase accumulation.

(ii) Tuning the laser frequencies in resonance with the two ions Note2 , i.e., $`\stackrel{~}{\omega }_D=0`$, a phase-shift gate associating with $`|00`$ and $`|11`$ can be implemented to cancel out the phase factors in Eq. (Non-adiabatic geometric quantum computation with trapped ions). In the resonance case the $`\omega _{L,1}+\omega _{L,2}`$ rotating frame coincides with the interaction picture of $`\stackrel{~}{}_0`$. As the one-bit resonant case, the effective magnetic field corresponding to Eq. (7) now lies in the $`x`$$`y`$ plane since $`\stackrel{~}{\omega }_D=0`$. By successively choosing two different values of the laser phase $`\mathrm{\Phi }`$, one can perform two $`\pi `$ rotations on the states $`|11`$ and $`|00`$ around the effective magnetic fields, and readily generate the A-A geometric phases, $`e^{i\stackrel{~}{\mathrm{\Gamma }}}`$ and $`e^{i\stackrel{~}{\mathrm{\Gamma }}}`$, for states $`|11`$ and $`|00`$, respectively. Now, after a phase-shift operation with $`\stackrel{~}{\mathrm{\Gamma }}=\stackrel{~}{\omega }_D\tau /2\pi /4`$, Eq. (Non-adiabatic geometric quantum computation with trapped ions) becomes

$`|00`$ $``$ $`e^{i\pi /4}\left[\text{cos}\mathrm{\Gamma }|00+i\text{sin}\mathrm{\Gamma }|11\right]`$
$`|11`$ $``$ $`e^{i\pi /4}\left[\text{cos}\mathrm{\Gamma }|11+i\text{sin}\mathrm{\Gamma }|00\right].`$ (11)
This is identical to Eq. (Non-adiabatic geometric quantum computation with trapped ions), except for the additional global phases. Obviously, these global phases will not appear if we first generate a phase shift of $`e^{i\pi /4}`$ on $`|00`$ and $`e^{i\pi /4}`$ on $`|11`$ at the same time, by the phase-shift gate just described, prior to the operation of Eq. (Non-adiabatic geometric quantum computation with trapped ions).

We have thus realized the two-bit gate $`U_{jk}`$ as defined in Eq. (Non-adiabatic geometric quantum computation with trapped ions). Together with the arbitrary one-bit operations (rotation and phase shift) described earlier, we can now readily implement the important CPS gate Eq. (8), whose role is to transform $`|11e^{i\pi }|11`$, while to keep other computational two-bit basis states unchanged.

Conclusion and Discussion. The proposed non-adiabatic geometric QC scheme based on the A-A phases is expected to overcome several drawbacks of the adiabatic schemes, namely, the slow evolution, need of refocusing to eliminate the dynamical phases, and continuous control over many fields to construct non-trivial loops in the parameter space. Viewing that the trapped-ion is one of the best exploited systems for quantum computation, and that the proposed scheme requires a relatively simple atomic level configuration, we suggest, as a first step, to exploit it as an interferometer for principle-proof of the non-adiabatic geometric A-A phase discussed in this work. We believe that the interference associated with the non-adiabatic A-A phase can be readily demonstrated by experiment in ion-trap systems.

As a possible QC architecture, the elementary operation steps in the proposed non-adiabatic geometric scheme are comparable to its dynamic counterparts. Specifically, the time scales for both the one-bit and two-bit geometric operations are about the same as those in the dynamic operations. Concerning the possible fault-tolerance, in the adiabatic case, quantum logic is tolerant to certain types of errors such as the field fluctuations that preserve the loop area in parameter space, i.e., the Berry phase. Similarly, in the non-adiabatic case, the A-A phase is of error tolerance to any fluctuations around the state-evolving path that preserve the path loop area. In principle, there exist many possible driving field deviations that can preserve the state path loop area. However, the most natural and possible errors appear to be random (but small) fluctuations of the laser phase, frequency, and coupling strength to the atomic levels, which equivalently result in fluctuations of the effective magnetic field. The global A-A phase is expected to be largely immune from this type of errors, and, at the same time, the net dynamic phase accumulation is approximately zero due to the cancellation of the positive and negative contributions.

Acknowledgments. Support from the Major State Basic Research Project No. G001CB3095 of China, the Special Fund for ”100 Person Project” from Chinese Academy of Sciences, the Youth Science and Technology ”Qimingxing” Program of Shanghai (No.99QA14036), and the Research Grants Council of the Hong Kong Government is gratefully acknowledged.


text:
# Production of 𝜋⁰ and 𝜂 mesons at large transverse momenta in pp and pBe interactions at 530 and 800 GeV/𝑐

## I INTRODUCTION

The study of inclusive single-hadron production at large transverse momentum ($`p_T`$) has been a useful probe in the development of perturbative quantum chromodynamics (PQCD) . Early in the evolution of the parton model, a departure from an exponential dependence of particle production at lower $`p_T`$ was interpreted in terms of the onset of interactions between pointlike constituents (partons) contained in hadrons. Large $`p_T`$ is a regime where perturbative methods have been applied to QCD to provide quantitative comparisons with data. Such comparisons yield information on the validity of the PQCD description, and on parton distribution functions of hadrons and fragmentation functions of partons.

This paper reports high-precision measurements of the production of $`\pi ^0`$ and $`\eta `$ mesons with large $`p_T`$. The $`\pi ^0`$ production cross sections are compared with next-to-leading order (NLO) PQCD calculations . As illustrated in a previous publication , our data (for both inclusive $`\pi ^0`$ and direct-photon production) are not described satisfactorily by the available NLO PQCD calculations, using standard choices of parameters. Similar discrepancies have been observed  between conventional PQCD calculations and other measurements of high-$`p_T`$ $`\pi ^0`$ and direct-photon cross sections (see also ). The origin of these discrepancies may be attributed to the effects of initial-state soft-gluon radiation. Such radiation generates transverse components of initial-state parton momenta, referred to below as $`k_T`$ . Evidence of significant $`k_T`$ in various processes, and a phenomenological model for incorporating its effect on calculated high-$`p_T`$ cross sections, have been extensively discussed in Refs. ; recent studies of the photoproduction of direct photons at HERA may provide additional insights . The inadequacy of NLO PQCD in describing $`k_T`$-sensitive distributions has been discussed in Ref. . In this paper, we follow the phenomenological prescription of Ref.  when comparing calculations with our $`\pi ^0`$ data. We also present cross sections for $`\eta `$ meson production at large-$`p_T`$. As might have been expected from earlier measurements (see, $`e.g.`$, ), $`\eta `$ production, relative to $`\pi ^0`$ production, shows little dependence on $`p_T`$ or on center of mass rapidity ($`y_{\mathrm{cm}}`$).

## II THE EXPERIMENTAL SETUP

The E706 experiment at Fermilab was designed to measure direct-photon production at high-$`p_T`$, and to investigate the structure of events containing direct photons. The data collection phase of the experiment spanned three fixed-target running periods, and included a relatively low statistics commissioning run in 1987-88 , and primary data runs in 1990 and 1991-92. The results presented here are from data recorded during the 1991-92 run. The E706 apparatus, operated in tandem with the E672 muon spectrometer, constituted the Meson West Spectrometer, displayed schematically in Fig. 1. The experiment used a right-handed Cartesian coordinate system, with the $`Z`$-axis pointed in the nominal beam direction, and the $`Y`$-axis pointed upward. The principal elements of the Meson West Spectrometer are discussed below. A more detailed description of the triggering methodology appears in Ref. .

### A Beamline and target

The Meson West beamline was capable of transporting either 800 GeV/$`c`$ primary protons from the Fermilab Tevatron or secondary beams of either polarity. We report here on results from studies using 800 GeV/$`c`$ primary protons and a 530 GeV/$`c`$ positive secondary beam. The beamline was instrumented with a differential Čerenkov counter  to identify incident pions, kaons, and protons in secondary beams. This helium-filled counter was 43.4 m long and was located $``$100 m upstream of the experiment’s target. Using the Čerenkov counter, the proton fraction at 530 GeV/$`c`$ was determined to be 97% .

A 4.7 m long stack of steel surrounding the beam pipe was placed between the last beamline magnet and the target box (see Fig. 1) to absorb hadrons. A water tank was placed at the downstream end of this hadron shield to absorb low-energy neutrons. During the 1991-92 run, two walls of scintillation counters were located both upstream and downstream of the hadron shield to identify penetrating muons.

The target region during the 1991-92 run consisted of two 0.8 mm thick copper disks of 2.5 cm diameter, located immediately upstream of a liquid hydrogen target, followed by a 2.54 cm long beryllium cylinder of diameter 2.54 cm. The hydrogen target consisted of a 15 cm long mylar flask, supported in an evacuated volume, with beryllium windows at each end (2.5 mm thickness upstream and 2.8 mm thickness downstream) .

### B Charged particle tracking

The spectrometer employed a charged particle tracking system consisting of silicon strip detectors (SSDs) , a dipole analysis magnet, proportional wire chambers (PWCs), and straw tube drift chambers (STDCs) . The SSD system consisted of sixteen planes of silicon wafers, arranged in eight modules. Each module contained two SSD planes, one providing $`X`$-view information and the other instrumenting the $`Y`$-view. Six 3$`\times `$3 cm2 SSD planes were located upstream of the target and used to reconstruct beam tracks. Two hybrid 5$`\times `$5 cm2 SSD planes (25 $`\mu `$m pitch strips in the central 1 cm, 50 $`\mu `$m beyond) were located $``$2 cm downstream of the Be target . These were followed by eight 5$`\times `$5 cm2 SSD planes of 50 $`\mu `$m pitch. The SSDs were instrumented to cover a solid angle of $`\pm 125`$ mr. Figure 2 displays the distribution of reconstructed vertices from a representative sample of 1991-92 data, showing clear separation of the different target elements.

The analysis dipole magnet imparted a transverse momentum impulse of $``$$`450\mathrm{MeV}/c`$ (in the horizontal plane) to singly-charged particles. Downstream track segments were measured by means of four stations of four views ($`XYUV`$) of 2.54 mm pitch PWCs and two stations of eight (4$`X`$4$`Y`$) layers of STDCs with tube diameters 1.03 cm (upstream station) and 1.59 cm (downstream station). The STDC stations were installed prior to the 1990 fixed-target run, and improved the angular resolution of the downstream tracking system (to $``$0.06 mrad) to make it comparable to that of the upstream system.

### C Calorimetry

The central element of the E706 apparatus was the finely-segmented liquid argon electromagnetic calorimeter (EMLAC) used to detect and measure electromagnetic showers . The EMLAC had a cylindrical geometry with an inner radius of 20 cm and an outer radius of 160 cm. It was divided into four mechanically independent quadrants, which were further subdivided electronically to create octants. The calorimeter had 33 longitudinal cells read out in two sections: an 11 cell front section ($``$8.5 radiation lengths) and a 22 cell back section ($``$18 radiation lengths). This front/back split was used for measuring the direction of incidence of showering particles, for discriminating between electromagnetic and hadronic showers and for resolving closely separated electromagnetic showers. The longitudinal cells consisted of 2 mm thick lead cathodes (the first cathode was constructed of aluminum), double-sided copper-clad G-10 radial ($`R`$) anode boards, followed by 2 mm thick lead cathodes and double-sided copper-clad G-10 azimuthal ($`\mathrm{\Phi }`$) anode boards. There were 2.5 mm argon gaps between each of these layers in a cell. The physical layout is illustrated in Fig. 3.

The copper-cladding on the anode boards was cut to form strips. Signals from corresponding strips from all $`R`$ (or $`\mathrm{\Phi }`$) anode boards in the front (or back) section were jumpered together. The copper-cladding on the radial anode boards was cut into concentric strips centered on the nominal beam axis. The width of the strips on the first $`R`$ board was 5.5 mm. The width of the $`R`$ strips on the following $`R`$ boards increased slightly so that the radial geometry was projective relative to the target, which was located 9 m upstream of the EMLAC. The azimuthal readout was subdivided at a radius of 40 cm into inner and outer segments, with each inner $`\mathrm{\Phi }`$ strip subtending an azimuthal angle of $`\pi /192`$ radians, and outer $`\mathrm{\Phi }`$ strips covering $`\pi /384`$ radians. Subdivision of the azimuthal strips in the outer portion of the detector improved both the position and energy resolution for showers reconstructed in this region. It also reduced $`R`$$`\mathrm{\Phi }`$ correlation ambiguities from multiple showers in the same octant of the calorimeter.

Data acquisition and trigger-signal processing for the EMLAC was based upon the FNAL RABBIT system . To achieve tolerable deadtimes, the zero suppression features of this system were used extensively during the experiment’s 1987-88 commissioning run. However, zero suppression limited our ability to understand the effects of out-of-time photon-induced showers and the tails of hadron-induced showers, thereby compromising efforts to understand the detailed response of the detector. Consequently, FASTBUS modules (the ICBM—Intelligent Control and Buffering Module—and the Wolf interface ) were developed by E706 to replace the original, CDF-designed, MX readout controllers. These FASTBUS modules enabled us to eliminate zero suppression during the experiment’s two primary data runs.

The apparatus also included two other calorimeters: a hadronic calorimeter (HALAC) located downstream of the EMLAC in the same cryostat, and a steel and scintillator calorimeter (FCAL), positioned further downstream, to increase coverage in the very forward region. The HALAC had 53 longitudinal cells read out in 2 sections: a 14 cell front section ($``$2 interaction lengths), and a 39 cell back section ($``$6 interaction lengths). Each cell consisted of read-out planes separated by 3 mm argon gaps and a 2.5 cm thick steel plate.

The FCAL acceptance covered the beam hole region of the EMLAC. It was split into three longitudinally similar sections. Each section was composed of alternating layers of 1.9 cm thick steel absorber plates and 4.8 mm thick acrylic scintillator sheets. The distance between steel plates was 6.9 mm. The downstream module contained 32 steel absorber plates and 33 scintillator sheets; the other two modules were comprised of 28 steel absorber plates and 29 scintillator sheets. Together, the three modules constituted $``$10.5 interaction lengths of material.

### D Muon identification

The E672 muon spectrometer, consisting of a toroidal magnet, shielding, scintillators, and proportional wire chambers, was deployed immediately downstream of the FCAL. The combined Meson West Spectrometer was triggered on high-mass muon pairs in order to investigate the hadroproduction of $`J/\psi `$, $`\psi (2s)`$, $`\chi _c`$, and B mesons . E706 and E672 collected data simultaneously and shared trigger logic, the data acquisition system, and event reconstruction programs. Data collected with the dimuon trigger were also used for several technical studies—for example, the $`J/\psi \mu ^+\mu ^{}`$ signal was used to calibrate the momentum scale of the tracking system.

### E Triggering

The E706 trigger selected events yielding high transverse momentum showers in the EMLAC. The selection process involved four stages: beam and interaction definitions, pretrigger requirements, and the final trigger requirements. Beam particles were detected using a hodoscope consisting of three planes (arranged in $`X`$, $`Y`$ and $`U`$ views) of scintillator, located $``$2 m upstream of the target region. Each plane contained 12 scintillator strips. The widths of the strips varied from 1 mm in the central region, to 5 mm along the edges. The edges of individual strips in each plane were overlapped to avoid gaps in acceptance. A beam signal was generated by a coincidence of signals from counters in at least two of the three hodoscope planes. A beam1 signal required that less than two hodoscope planes detected two or more isolated clusters of hits in coincidence with beam. This beam1 requirement rejected events with multiple beam particles incident upon the target. A plane of four scintillation counters (referred to as beam hole counters) arranged to produce a 0.95 cm diameter central hole was located downstream of the beam hodoscope and used to reject interactions initiated by particles in the beam halo.

Two pairs of scintillation counters were mounted on the dipole analysis magnet, one pair upstream and the other downstream of the magnet. Each pair had a central hole that allowed non-interacting beam particles to pass through undetected. An interaction was defined as a coincidence between signals from at least two of these four interaction counters. To minimize potential confusion in the EMLAC due to out-of-time interactions, a filter was used to reject interactions that occurred within 60 ns of one other.

For those interactions that satisfied both the beam1 and interaction definitions, the $`p_T`$ deposited in various regions of the EMLAC was evaluated by weighting the energy signals from the fast outputs of the EMLAC $`R`$ channel amplifiers by $``$$`\mathrm{sin}\theta _i`$, where $`\theta _i`$ is the polar angle that the $`i^{th}`$ strip subtends with respect to the nominal beam axis. The pretrigger hi requirement was satisfied when the $`p_T`$ detected either in the inner 128 $`R`$ channels or the outer $`R`$ channels of any octant was greater than the threshold value of $``$2 GeV/$`c`$. A pretrigger signal was issued only if the signals from a given octant satisfied the pretrigger requirement, there was no evidence in that octant of substantial noise or significant $`p_T`$ attributable to an earlier interaction, and there was no incident beam halo muon detected by the walls of scintillation counters surrounding the hadron shield upstream of the spectrometer. The pretrigger signal latched data from the various subsystems until a final trigger decision was made.

Localized trigger groups were formed for each octant by clustering the $`R`$ channels into 32 groups of 8 channels. Each of the adjacent pairs of 8 channel groups (groups 1 and 2, 2 and 3, $`\mathrm{}`$ , 31 and 32) formed a local group of 16 strips. If the $`p_T`$ detected in any of these groups of 16 was above a specified high (or low) threshold, then a local hi (or local lo) signal was generated for that octant. The single local hi (single local lo) trigger required a local hi (local lo) signal from an octant that also satisfied the pretrigger hi. The local hi (local lo) threshold was $``$3.5 GeV/$`c`$ ($``$2 GeV/$`c`$). The single local lo trigger was prescaled by a factor of $``$200.

In addition to these high-$`p_T`$ triggers, prescaled samples of low-bias beam, interaction, and pretrigger events were also recorded. The prescale factors for these triggers were typically set at $`15^6`$, $`15^5`$, and $`15^3`$, respectively. These low-bias triggers constituted $``$10$`\%`$ of the events recorded.

## III $`\pi ^0`$ AND $`\eta `$ ANALYSIS

The data sample used in this analysis corresponds to an integrated luminosity of 6.8 (1.2) events/pb for 530 GeV/$`c`$ $`\mathrm{p}`$Be ($`\mathrm{p}\mathrm{p}`$) interactions, and 6.5 (1.2) events/pb at 800 GeV/$`c`$. The following subsections describe the analysis procedures and methods used to correct the data for losses resulting from inefficiencies and selection biases.

### A Event reconstruction

The charged-track reconstruction algorithm produced track segments upstream of the dipole magnet using information from the SSDs, and downstream of the magnet using information from the PWCs and STDCs. These track segments were projected to the center of the magnet, and linked to form the final tracks, and the interaction vertex. The charged-track reconstruction and vertex-finding methodology are described in Ref. .

The readout in each EMLAC quadrant consisted of four regions: left $`R`$ and right $`R`$ (radial strips of each octant in that quadrant), and inner $`\mathrm{\Phi }`$ and outer $`\mathrm{\Phi }`$ regions (azimuthal strips divided at $`R=40\mathrm{cm}`$). Strip energies from clusters in each region were fitted to the shape of electromagnetic showers, as determined from Monte Carlo simulations and isolated-shower data. These fits were used to evaluate the positions and energies ($`E_R`$ and $`E_\mathrm{\Phi }`$) of the peaks in each region. Shower positions and energies were obtained by correlating peaks of approximately the same energy in the $`R`$ and $`\mathrm{\Phi }`$ regions within the same half octant (more complex algorithms were used to handle configurations with overlapping showers in either the $`R`$ or $`\mathrm{\Phi }`$ regions). Any differences in photon energy as measured in the $`R`$ and $`\mathrm{\Phi }`$ views reflect the intrinsic resolution properties of the calorimeter, and provide a test of the quality of the Monte Carlo simulations. The EMLAC’s longitudinal segmentation provided discrimination between showers generated by electromagnetically or hadronically interacting particles. For individual showers, the ratio of energy reconstructed in the front section to the sum of energy in the front and back section of the EMLAC (referred to as $`E_{\mathrm{FRONT}}/E_{\mathrm{TOTAL}}`$) also tested the Monte Carlo simulation of longitudinal shower development (see the detector simulation section below). An expanded discussion of the EMLAC reconstruction procedures and performance can be found in Ref. .

### B Data sample selection and corrections

The events contributing to the cross section measurements were required to have a reconstructed vertex within the fiducial volume of the Be or H2 targets . Vertex reconstruction efficiencies were evaluated for each target using a detailed Monte Carlo simulation of the spectrometer (described below). These efficiencies were used to correct for reconstruction losses and for resolution smearing across target fiducial boundaries. The vertex reconstruction efficiency was $`1`$ for the H2 and downstream Be targets and $`0.97`$ for the upstream Be target.

Both $`\pi ^0`$ and $`\eta `$ mesons were reconstructed via their $`\gamma \gamma `$ decay modes. The photons were required to be within the fiducial region of the EMLAC, which excluded areas with reduced sensitivity. In particular, photons incident upon regions of the detector near quadrant boundaries (which abutted steel support plates), the central beam hole, the outer radius of the EMLAC, and octant boundaries were excluded from consideration. A simple ray-tracing Monte Carlo program was employed to determine the correction for the effect of the azimuthal fiducial boundaries .

To reduce backgrounds from hadronic showers, only showers with at least 20$`\%`$ of their shower energy deposited in the front part of EMLAC ($`E_{\mathrm{FRONT}}/E_{\mathrm{TOTAL}}`$$`>`$ 0.2) were accepted as photon candidates. To correct the cross sections for this requirement, as well as for other larger effects including resolution-smearing and reconstruction losses, a full simulation of the showers in the calorimeter (described below) was employed. In addition, $`\gamma \gamma `$ combinations were considered as $`\pi ^0`$ or $`\eta `$ candidates only when the two photons were detected in the same octant (to simplify subsequent analysis of the trigger response), and only those combinations with energy asymmetry \[$`A_{\gamma \gamma }|E_{\gamma _1}E_{\gamma _2}|/(E_{\gamma _1}+E_{\gamma _2})`$\] less than 0.75 were considered (to reduce uncertainties due to low energy photons).

The $`\gamma \gamma `$ invariant mass distribution in the $`\pi ^0`$ and $`\eta `$ regions for photon pairs that satisfied the above requirements is shown in Fig. 4 for representative low $`p_T`$ intervals. A $`\pi ^0`$ candidate was defined as a combination of two photons with invariant mass, $`M_{\gamma \gamma }`$ , in the range $`100\mathrm{MeV}/c^2<M_{\gamma \gamma }<180\mathrm{MeV}/c^2`$. An $`\eta `$ candidate was defined as a two-photon combination in the range $`450\mathrm{MeV}/c^2<M_{\gamma \gamma }<650\mathrm{MeV}/c^2`$. To evaluate production cross sections, combinatorial background in the $`\pi ^0`$ and $`\eta `$ regions was evaluated as follows. Sideband regions were defined to cover an equivalent mass range of the $`\pi ^0`$ and $`\eta `$ peak regions (using the same acceptance criteria as for the peak regions). The $`p_T`$ and rapidity distributions from these side bands were then subtracted from the corresponding distributions within the $`\pi ^0`$ and $`\eta `$ mass ranges to obtain the respective signals. This technique is appropriate as long as the combinatorial background depends approximately linearly dependence on $`M_{\gamma \gamma }`$. At low $`p_T`$ (below $`2`$ GeV/$`c`$), the shape of the combinatorial background in the signal regions is not linear, and a more sophisticated fitting procedure was used to evaluate the background. The $`\gamma \gamma `$ mass distributions were fitted to Gaussians for signal plus second-order polynomials in $`M_{\gamma \gamma }`$ to represent the background. The combinatorial background in the peak regions was then defined using the resultant parameters of the fit, and the signals defined as the differences between the totals and the fitted backgrounds.

For the cross section measurements, the signals have been corrected for losses due to the energy asymmetry cut and the branching fractions  for the $`\gamma \gamma `$ decay modes. The correction for losses due to the conversion of one or both of the photons into $`e^+e^{}`$ pairs was evaluated by projecting each reconstructed photon from the event vertex to the reconstructed position in the EMLAC. The number of radiation lengths of material traversed along the photon path was calculated on the basis of a detailed description of the detector. The photon conversion probability was then evaluated and used to account for losses from conversion.

### C Trigger response

As mentioned previously, the single local hi and single local lo trigger decisions were based upon depositions of $`p_T`$ in the EMLAC within groups of 16 contiguous radial strips. For each $`\pi ^0`$ or $`\eta `$ candidate, a probability to satisfy the trigger was defined, based upon energy deposition in the entire octant, as $`P=1(1p_\mathrm{i})`$, where $`p_\mathrm{i}`$ is the efficiency of the $`\mathrm{i}^{th}`$ trigger group in the octant containing the candidate . The inverse of this probability was applied as a trigger weight to each meson candidate. To avoid excessively large trigger weights, meson candidates with trigger probabilities of $`<`$ 0.1 were excluded from further consideration. The correction for losses from this requirement was determined from Monte Carlo, and absorbed into the reconstruction efficiency.

The efficiency of the pretrigger hi was determined in a manner similar to that used for the local triggers . Additional details on the trigger can be found in Refs. .

The cross sections presented in this paper reflect composite measurements, utilizing a combination of results from the interaction, pretrigger hi, single local lo, and single local hi triggers. The $`p_T`$ spectra (corrected for only prescale factors) from a representative sample of these triggers are shown in Fig. 5. The transition points chosen between the high and low threshold triggers were determined by comparing the fully corrected results from each trigger, and were different for $`\pi ^0`$ and $`\eta `$ mesons and also depended on rapidity.

### D Rejection of beam halo muons

Spurious triggers were produced by muons in the beam halo that radiated energy in the electromagnetic calorimeter in random coincidence with an interaction in the target. Particularly in the outer regions of the EMLAC, this energy was recorded as a high-$`p_T`$ deposition that satisfied the local trigger requirements. This occurred much more frequently in data from the 530 GeV/$`c`$ secondary beam than for the primary 800 GeV/$`c`$ beam, due to the absence of an upstream interaction target in the latter case. To reduce this background, the pretrigger logic relied on signals from the veto walls of scintillator counters to reject events associated with such muons in the beam halo. In the off-line analysis, we employed expanded requirements on the latched veto-wall signals, the direction of reconstructed showers , the shower shape (halo muon-induced showers have a different shape than photon or electron-induced showers), and the total $`p_T`$ (balance) in the event. For the latter, we calculated the net $`p_T`$ of the photons and charged particles which, based upon their initial trajectories, would intercept the EMLAC in the transverse plane within the $`120^{}`$ sector opposite the meson candidate. In interactions which generate a high-$`p_T`$ meson, the ratio of this “away-side” $`p_T`$ ($`P_T^{away}`$) to that of the meson $`p_T`$ should be near unity. However, for events triggered by showers from beam halo muons, $`P_T^{away}/p_T`$ should be near zero, since the interaction in random coincidence with the beam halo muon is typically a soft (low $`p_T`$) interaction. Candidates with $`P_T^{away}/p_T<0.3`$ were considered likely to be due to beam halo muons and were rejected.

To illustrate the effect of the above off-line requirements, the $`\gamma \gamma `$ invariant mass distribution, both before and after application of the rejection criteria, is shown in Fig. 6 for the 530 and 800 GeV/$`c`$ data, for $`p_T>7.0`$ GeV/c. The large muon-induced background at low $`\gamma \gamma `$ mass values in the 530 GeV/$`c`$ data is due to the occasional transverse splitting of the muon-induced showers into two closely-separated photon candidates. This happens because the reconstruction software assumes that the showers originate from the target-region rather than from the beam halo. The 800 GeV/$`c`$ data have very few muon-induced triggers, and is consequently not affected very strongly by these rejection criteria.

The impact of these rejection criteria on the physics signal was checked using more restrictive selection criteria to define a pure sample of $`\gamma \gamma `$ pairs. The fraction of signal lost by the application of each of the muon-rejection requirements determined a correction to the cross section. The product of the correction factors for muon rejection corresponds to an increase of $``$8$`\%`$ in the cross secton at $`p_T`$ = 4 GeV/$`c`$, and $``$10$`\%`$ at $`p_T=7`$ GeV/$`c`$, for the 530 GeV/$`c`$ beam data. (The corrections were smaller for the 800 GeV/$`c`$ data because of cleaner beam conditions.)

### E Detector simulation

The Meson West spectrometer was modeled using a detailed geant  simulation. Because the full simulation of electromagnetic showers requires extensive computing time, we developed a hybrid approach using geant-tracking through the magnetic spectrometer and in the initial stages of the shower development in the calorimeter. We used the standard geant algorithms for tracking particles with energies above 10 MeV, below which we relied on an empirical parameterization for the deposition of energy in the EMLAC . This cutoff was selected to be at the point at which bremsstrahlung still dominates the energy loss in lead, and led to significant improvement in processing speed. In doing this, we took advantage of the steady advances in computational power of the FNAL UNIX farms  to reach the desired level of statistical accuracy.

As inputs to the geant simulation, we employed single particle distributions, reconstructed data events, and two physics event generators: herwig  and pythia . For the analyses described herein, we chose herwig as the principal Monte Carlo event generator based on a better match of particle multiplicities between data and Monte Carlo using the default parameters. Over 5.5 million herwig events were passed through the geant simulation. We weighted the herwig $`\pi ^0`$ and $`\eta `$ spectra (in $`p_T`$ and rapidity) to match our measured results, so that the corrections obtained from the Monte Carlo were based on the data distributions rather than on the behavior of the physics generator.

The calibration of the energy response of the EMLAC was based on the reconstructed masses of $`\pi ^0`$ mesons in the $`\gamma \gamma `$ decay mode . The steeply falling $`p_T`$ spectrum for $`\pi ^0`$ production, combined with the calorimeter’s resolution, produced a small offset ($``$1$`\%`$) in the mean reconstructed photon energies. Using the same calibration procedure in the simulated EMLAC as in the detector, we corrected this offset and minimized any potential biases in the calibration. We also employed the geant Monte Carlo simulation to evaluate the mean correction (as a function of photon energy) for energy deposited in the material upstream of the EMLAC. The impact of detector resolution on the energy scale and on the $`\pi ^0`$ and $`\eta `$ production spectra was incorporated in the overall reconstruction efficiency corrections.

To ensure that the Monte Carlo simulation reproduced the data, a special preprocessor was used to convert geant information into signals and strip energies as measured in the various detectors, and to simulate hardware effects, such as channel noise and gain variations. The generated Monte Carlo events were then processed through the same reconstruction software used for the analysis of data, and thereby provided measures of the inefficiencies and biases of the reconstruction algorithms.

Comparisons between results from the Monte Carlo simulation and the data for the distributions in $`E_RE_\mathrm{\Phi }`$ and in the fraction of energy deposited in the front section of the EMLAC, $`E_{\mathrm{FRONT}}/E_{\mathrm{TOTAL}}`$, are presented in Figs. 7 and 8. The Monte Carlo results are in satisfactory agreement with the data, indicating that the simulation properly treats shower development in the EMLAC. Figures 9 and 10 show the $`\gamma \gamma `$ mass spectra in the $`\pi ^0`$ and $`\eta `$ mass regions for two minimum-$`p_T`$ requirements, and compare these spectra to the Monte Carlo simulation. In addition to giving further evidence that the Monte Carlo provides a good simulation of the resolution of the EMLAC, the agreement in the levels of the combinatorial background indicates that the Monte Carlo also provides reasonable simulation of the underlying event structure. Figure 11 shows a comparison between the Monte Carlo and the data for the sideband-subtracted $`\pi ^0`$ energy asymmetry distribution, and the agreement indicates that the Monte Carlo simulation describes accurately the losses of low-energy photons.

Reconstruction inefficiencies for $`\pi ^0`$ and $`\eta `$ mesons (which satisfied the $`A_{\gamma \gamma }`$, energy asymmetry, requirement) were relatively small over most of the kinematic range. Figure 12 shows the probability for a $`\pi ^0`$ to pass the selection requirements imposed on the Monte Carlo events at 530 GeV/$`c`$, as a function of $`p_T`$, for different rapidity intervals. This probability includes losses due to the reconstruction algorithm, the $`E_{\mathrm{FRONT}}/E_{\mathrm{TOTAL}}`$ requirement, and the 10$`\%`$ minimum trigger probability requirement. The inefficiency at forward rapidities and high-$`p_T`$ is attributable to the increased difficulty in separating two photons from $`\pi ^0`$ decays (coalescence) in this kinematic region.

### F Normalization

Electronic scalers that counted signals from the beam hodoscope, interaction counters, and beam hole counters were used to determine the number of beam particles incident on the target. Other scalers logged the state of the trigger and of components of the data acquisition system. Information from these scalers was used to determine the number of beam particles that traversed the spectrometer when it was ready to record data. This number was corrected for multiple occupancy in the beam hodoscope (beyond that excluded via the beam1 requirement, and for absorption of beam in the target material.

The normalization of the low $`p_T`$ $`\pi ^0`$ cross section was independently verified using events from the prescaled beam and interaction trigger samples. In these samples, the absolute normalization can be obtained just by counting events. For these low $`p_T`$ events, the normalization as determined via the scalers and via event counting techniques agree to better than $`3\%`$.

Based upon the good agreement between results from these independent normalization methods, combined with the stability of the cross section results from various sections of the run, an evaluation of the internal consistency of the scalers, and a detailed analysis of the design, implementation and performance of the trigger, the net systematic uncertainty in the overall normalization is $``$8$`\%`$.

### G Secondary Beam Contamination

The 530 GeV/$`c`$ secondary beam cross section measurements were corrected for the small admixture of pions and kaons present in the beam, estimated as 2.75$`\%`$ $`\pi ^+`$ and 0.5$`\%`$ $`K^+`$, respectively . Although the percentage contamination is small, it’s effect at high $`p_T`$ is enhanced by the presence of two, rather than three, valence quarks in incident mesons. The effect of $`\pi ^+`$ contamination was estimated using our high statistics study of $`\pi ^0`$ and $`\eta `$ meson production by 515 GeV/$`c`$ $`\pi ^{}`$ beam . This is justified because: neutral meson production by $`\pi ^+`$ and $`\pi ^{}`$ beams is expected to be similar, based on isospin arguments and earlier measurements ; the $`\pi ^+`$ component of the nominally 530 GeV/$`c`$ positive secondary beam had a mean momentum of 515 GeV/$`c`$; and the ratio of the measured high $`p_T`$ $`\pi ^0`$ cross section from the Čerenkov-tagged $`\pi ^+`$ component of the positive secondary beam to the corresponding cross section from the 515 GeV/$`c`$ $`\pi ^{}`$ beam was consistent with unity .

The effect of $`K^+`$ contamination was assumed to be half that of the same amount of $`\pi ^+`$ contamination, consistent with earlier, lower energy measurements  and with our own, more statistically limited, data . After correcting for beam contamination, the 530 GeV/$`c`$ cross sections were reduced by $`2\%`$ at low $`p_T`$ and by $`10\%`$ at high $`p_T`$.

### H Summary of systematic uncertainties

The principal contributions to the systematic uncertainty arose from the following sources: calibration of photon energy response $`(5`$$`9\%)`$, $`\pi ^0`$ and $`\eta `$ reconstruction efficiency and detector-resolution unsmearing $`(5\%)`$, the overall normalization $`(8\%)`$, and, for the 530 GeV/$`c`$ secondary beam, the beam contamination $`(0`$$`7\%)`$. A more complete list of systematic uncertainties is presented in Table I. Some of these uncertainties ($`e.g.`$ normalization) are strongly correlated between bins. The systematic uncertainties, combined in quadrature, are quoted with cross sections in the appropriate tables. Combining all systematic uncertainties yields a net uncertainty that varies from $`11\%`$ ($`13\%`$) at low $`p_T`$ for $`\pi ^0`$ ($`\eta `$) mesons to $`15\%`$ at high-$`p_T`$.

The secondary proton beam was determined to have a mean momentum of $`530\pm 2\mathrm{GeV}/c`$ with an estimated halfwidth of $`30\mathrm{GeV}/c`$. This momentum spread introduces a small uncertainty ($`5\%`$) in comparisons of theory with data. (For the 800 GeV primary beam, the momentum bite and the corresponding uncertainty are negligible.)

## IV INCLUSIVE CROSS SECTIONS

### A $`\pi ^0`$ production

The inclusive $`\pi ^0`$ cross sections per nucleon versus $`p_T`$ are shown in Figs. 13 and 14 for 530 and 800 GeV/$`c`$ proton beams, respectively, incident upon beryllium and liquid hydrogen targets. Because of the steeply falling spectra, the data are plotted at abscissa values that correspond to the average values of the cross section in each $`p_T`$ bin (assuming local exponential $`p_T`$ dependence) . These cross sections are also tabulated in Tables II and III. The corresponding cross section measurements, as functions of $`p_T`$ and $`y_{\mathrm{cm}}`$ are reported in Tables IV through VII.

NLO PQCD calculations  are compared to the data in Figs. 15 through 23. In Fig. 15, NLO PQCD results using CTEQ4M parton distribution functions  and BKK fragmentation functions  are compared to the measured inclusive $`\pi ^0`$ cross sections for $`\mathrm{p}`$Be and $`\mathrm{p}\mathrm{p}`$ interactions at 800 GeV/$`c`$. The PQCD calculations for the Be target have been adjusted to account for nuclear effects by the factor $`A^{\alpha 1}`$, where $`A`$ is the atomic number and $`\alpha =1.12(1.08)`$ for 530 (800) GeV/c incident protons. Theoretical results are presented for three values of the factorization scale: $`\mu `$ = $`p_T/2`$, $`p_T`$, and $`2p_T`$. (In these comparisons, the renormalization and fragmentation scales have been set to the value of the factorization scale.) In addition to a large dependence on the choice of scale, the expectations for these choices of $`\mu `$ lie significantly below the data, at both 530 and 800 GeV/$`c`$.

In Fig. 16, NLO calculations using BKK and KKP  fragmentation functions, and $`\mu `$ = $`p_T/2`$, are compared to the $`\pi ^0`$ cross sections for 530 GeV/$`c`$ $`\mathrm{p}`$Be and $`\mathrm{p}\mathrm{p}`$ interactions. The calculations exhibit considerable dependence on the choice of fragmentation function, but both choices predict yields that are significantly lower than the data.

These discrepancies have been interpreted  as arising from additional soft-gluon emission in the initial state that is not included in the NLO calculation, and which results in sizeable parton $`k_T`$ before the hard collision (for a different perspective, see the discussion in Ref. ). Soft-gluon (or $`k_T`$) effects are expected in all hard-scattering processes, such as the inclusive production of jets, high-$`p_T`$ mesons, and direct photons . The Collins-Soper-Sterman resummation formalism  provides a rigorous basis for understanding these radiation effects, and there have been several recent efforts to derive resummation descriptions for the inclusive direct-photon  and dijet cross sections . The calculation of Ref.  for inclusive direct-photon production, which includes the effects of soft-gluon resummation near the kinematic threshold limit $`\left(x_T=2p_T/\sqrt{s}1\right)`$, has a far smaller sensitivity to scale, compared to NLO calculations, and provides cross sections close to those of NLO calculations with a scale of $`\mu =p_T/2`$. Also, for our energies, the calculation of Ref.  using $`k_T`$-resummation, and of Ref. , which treats simultaneously threshold and recoil effects in direct-photon production, yield a substantially larger cross section than the NLO result. However, no such calculations are available for inclusive meson production. In their absence, we use a PQCD-based model that incorporates transverse kinematics of initial-state partons to study the principal consequences of additional $`k_T`$ for high-$`p_T`$ production processes.

Because the unmodified PQCD cross sections fall rapidly with increasing $`p_T`$, the net effect of the “$`k_T`$ smearing” is to increase the expected yield at higher $`p_T`$’s. An exact treatment of the modified parton kinematics has been implemented in a Monte Carlo calculation of the leading-order (LO) cross sections for high-$`p_T`$ particle production , with the $`k_T`$ distribution for each of the incoming partons represented by a Gaussian with one adjustable parameter ($`k_T`$). Unfortunately, no such program is available for NLO calculations, and so we approximate the effect of $`k_T`$ smearing by multiplying the NLO cross sections by the corresponding LO $`k_T`$-enhancement factors. Admittedly, this procedure involves a risk of double-counting since some of the $`k_T`$-enhancement may already be contained in the NLO calculation. However, we expect such double-counting effects to be small.

The $`k_T`$ values used in the calculation of the LO $`k_T`$-enhancement factors are similar to those employed in comparisons of kinematic distributions in data involving production of high-mass $`\gamma \gamma `$, $`\gamma `$$`\pi ^0`$, and $`\pi ^0`$$`\pi ^0`$ systems with calculations relying on the same LO program (see Refs.  for further details). For these comparisons, we used the LO versions of the CTEQ4 distribution and (where appropriate) BKK fragmentation functions, and an average transverse momentum of 0.6 GeV/$`c`$ for the $`\pi ^0`$ mesons relative to the fragmenting parton direction (varying this parameter in the range 0.3 to 0.7 GeV/$`c`$ does not affect our conclusions) .

Comparisons of the $`k_T`$-enhanced calculations with data at 530 GeV/$`c`$ are displayed in Fig. 17, indicating reasonable agreement for the $`k_T`$ values chosen. Similar conclusions can be drawn from comparisons between calculations and data at 800 GeV/$`c`$, as illustrated in Fig. 18. It is interesting to compare the fractional differences between data and the $`k_T`$-enhanced NLO calculations using BKK and KKP fragmentation functions on a linear scale. Such comparisons are shown in Fig. 19 for $`\mathrm{p}`$Be interactions at 530 and 800 GeV/$`c`$. The $`k_T`$-enhanced calculations using the KKP fragmentation functions are seen to reproduce the shape of the $`\pi ^0`$ cross section better than calculations using the earlier BKK fragmentation functions, but small discrepancies are still present. Similar conclusions can be drawn from comparisons of NLO QCD with our data on $`\mathrm{p}\mathrm{p}`$ collisions (Fig. 20), although in this case the fractional differences using the KKP fragmentation functions are systematically greater than zero when we employ the same $`k_T`$ values used in Fig. 19. It is, however, worth noting that the hydrogen nucleus is not typical, having an up to down quark ratio of two rather than a value of less than one as in the case of standard nuclei, such as Be (Cu) for which $`\mathrm{u}/\mathrm{d}`$ = 0.93 (0.94). The theoretical calculations address this difference by calculating the nuclear cross sections using the actual numbers of protons and neutrons present in each case, but this alone does not fully compensate for the influence of the nuclear media surrounding interactions in targets other than hydrogen. As a result, values of $`k_T`$ that yield agreement with nuclear data adjusted by $`A^{\alpha 1}`$ may not be fully appropriate for hydrogen. In particular, raising $`k_T`$ by only 0.05 GeV/$`c`$ drops the points in the lower two quadrants of Fig. 20 by enough ($`0.075`$ at 4.5 GeV/c) to approximately center these distributions on zero.

Figures 21 and 22 show the cross sections for inclusive $`\pi ^0`$ production as functions of rapidity for $`\mathrm{p}`$Be interactions at 530 and 800 GeV/$`c`$, for several intervals in $`p_T`$. The expected peaking at a scattering angle near $`90^{}`$ in the center of mass ($`y_{\mathrm{cm}}`$ = 0) develops slowly as a function of $`p_T`$. The shapes and normalizations of the data are in good agreement with the $`k_T`$-enhanced calculations.

Both theoretical and experimental uncertainties are reduced in the ratio of invariant cross sections for $`\pi ^0`$ production at 800 and 530 GeV/$`c`$, allowing, in principle, a more sensitive test of the calculations. Figure 23 displays this ratio compared to conventional ($`k_T`$=0) and $`k_T`$-enhanced NLO results using KKP fragmentation functions. The $`k_T`$ values used here correspond to those used in Figs. 17 and 18. The energy dependence of the data is accommodated better by the $`k_T`$-enhanced theory. Similar conclusions are obtained for the hydrogen target data (not shown).

The results discussed in this section are not very sensitive to the reasonably well-known parton distribution functions  (quark distributions being of primary importance here). Methods similar to the ones described in this paper have been applied to analyze high-$`p_T`$ hadron spectra in $`\mathrm{p}\mathrm{p}`$, $`\mathrm{p}\mathrm{A}`$, and AA collisions , and $`k_T`$ effects have been found important for describing data on inclusive production of charged mesons.

### B $`\eta `$ production

Cross sections for inclusive $`\eta `$ production are tabulated in Tables VIII through XIII. Theoretical descriptions of $`\eta `$-meson production differ from the $`\pi ^0`$ case primarily because of differences in the fragmentation of partons into the particles of interest. To investigate this aspect, we present $`\eta /\pi ^0`$ relative production rates as functions of $`p_T`$ and $`y_{\mathrm{cm}}`$ (for two $`p_T`$ ranges) in Figs. 24 and 25 — the data average to a value of $`0.45\pm 0.01`$ at 530 GeV/$`c`$, and $`0.42\pm 0.01`$ at 800 GeV/$`c`$ (for 3 $`<p_T<`$ 8 GeV/$`c`$).

## V SUMMARY

The invariant cross sections for $`\pi ^0`$ and $`\eta `$ production have been measured for $`\mathrm{p}\mathrm{p}`$ and $`\mathrm{p}`$Be collisions at 530 and 800 GeV/$`c`$ as functions of $`p_T`$ and $`y_{\mathrm{cm}}`$, over the kinematic range $`1<p_T<10`$ GeV/$`c`$, and 1.5 units in rapidity. Results from $`k_T`$-enhanced NLO QCD calculations are in reasonable agreement with our measured $`\pi ^0`$ cross sections. Employing the recent KKP fragmentation functions in the NLO QCD calculations was found to improve the description of the detailed shape of our $`\pi ^0`$ cross sections (as a function of $`p_T`$) relative to theoretical results obtained using the earlier BKK fragmentation functions. The measured $`\eta /\pi ^0`$ production ratios, which provide information about the relative fragmentation of partons into these mesons, are $`0.45\pm 0.01`$ at 530 GeV/$`c`$ and $`0.42\pm 0.01`$ at 800 GeV/$`c`$, averaged over $`p_T`$ and rapidity.

###### Acknowledgements.

We thank the U. S. Department of Energy, the National Science Foundation, including its Office of International Programs and the Universities Grants Commission of India, for their support of this research. The staff and management of Fermilab are thanked for their efforts in making available the beam and computing facilities that made this work possible. We are also pleased to acknowledge the contributions of our colleagues on Fermilab experiment E672. We acknowledge the contributions of the following colleagues during the upgrade and operation of the Meson West spectrometer: W. Dickerson, E. Pothier from Northeastern University; J. T. Anderson, E. Barsotti, Jr., H. Koecher, P. Madsen, D. Petravick, R. Tokarek, J. Tweed, D. Allspach, J. Urbin, R. L. Schmitt and the cryo crews from Fermi National Accelerator Laboratory; T. Haelen, C. Benson, L. Kuntz, and D. Ruggiero from the University of Rochester; the technical staffs of Michigan State University and Pennsylvania State University for the construction of the straw tubes and of the University of Pittsburgh for the silicon detectors. We would also like to thank the following commissioning run collaborators for their invaluable contributions to the hardware and software infrastructure of the original Meson West spectrometer: G. Alverson, G. Ballocchi, R. Benson, D. Berg, D. Brown, D. Carey, T. Chand, C. Chandlee, S. Easo, W. Faissler, G. Glass, I. Kourbanis, A. Lanaro, C. A. Nelson, Jr., D. Orris, B. M. Rajaram, K. Ruddick, A. Sinanidis, and G. Wu. We thank S. Catani, J.Ph. Guillet, B. Kniehl, J. Owens, G. Sterman, and W. Vogelsang for many helpful discussions and for providing us with their QCD calculations.