text:
# 1 Introduction

## 1 Introduction

In general, neutral Higgs bosons do not couple directly to massless photons. For instance the standard model Higgs boson couples to photons only through loops of charged particles, i.e., W’s, quarks and leptons, and the branching ratio into $`\gamma \gamma `$ is small ($`10^3`$  for $`\mathrm{m}_\mathrm{H}90`$ GeV/$`\mathrm{c}^2`$). However the $`\gamma \gamma `$ branching fraction of Higgs bosons can be increased with respect to the standard model prediction in each of the following four configurations.

* The direct couplings to fermions are suppressed, as is the case for models with at least two Higgs multiplets , of which one couples only to fermions and the others only to gauge bosons. The physical states with couplings only to gauge bosons are called fermiophobic Higgs bosons.
* The direct couplings to gauge bosons are enhanced with anomalous couplings . These couplings are described in the most general formulation with four effective six-dimensional operators with strength $`\mathrm{f}_i/\mathrm{\Lambda }^2`$, where $`\mathrm{\Lambda }`$ is the scale of the new underlying interaction.
* Couplings to both fermions and bosons are modified, as is the case in the minimal supersymmetric extension of the standard model (MSSM).
* Additional light, charged particles enter the loops that couple Higgs bosons and photons, as is again the case in the MSSM (charginos, squarks, sleptons, charged Higgs bosons).

With some particular choices of parameters, the branching ratio into $`\gamma \gamma `$ may be enhanced in the MSSM, and can reach a value close to $`100\%`$ in models with fermiophobia or with anomalous gauge couplings. It is therefore possible that a Higgs boson has escaped the standard search for the Higgs-strahlung process $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Hf}\overline{\mathrm{f}}`$ with $`\mathrm{H}\mathrm{b}\overline{\mathrm{b}}`$ . In this letter, a complementary search for the Higgs-strahlung process with $`\mathrm{H}\gamma \gamma `$ is described.

The analysis addresses all topologies arising from the $`\mathrm{e}^+\mathrm{e}^{}\mathrm{H}Z^{()}`$ process, as characterized by the charged track multiplicity of the final state: (i) acoplanar photons with missing energy and no charged particles for $`\mathrm{H}\nu \overline{\nu }`$; (ii) photon pairs with exactly two charged particles identified as leptons for $`\mathrm{H}\mathrm{}^+\mathrm{}^{}`$; (iii) photon pairs accompanied with two thin, low multiplicity jets for $`\mathrm{H}\tau ^+\tau ^{}`$, from two to four charged particles; and (iv) photon pairs with a hadronic system for $`\mathrm{Hq}\overline{\mathrm{q}}`$ with at least five charged particles.

The analysis is performed with the data collected with the ALEPH detector from 1991 to 1999 including the Z peak data collected during the LEP 2 period. This sample corresponds to an integrated luminosity of $`672\mathrm{pb}^1`$ at centre-of-mass energies ranging from 88 to 202 GeV. Details are given in Table 1.

After a short description of the detector properties relevant for this search, the common preselection based on photon identification is reviewed in Section 3. The global search strategy is developed in Section 4. The systematic uncertainties affecting the selection efficiency are discussed in Section 5 and the results are given in Section 6.

## 2 The ALEPH detector

The ALEPH detector and its performance are described in Refs. . The tracking detectors, composed of the silicon vertex detector surrounded by the inner tracking chamber and the time projection chamber (TPC), provide efficient reconstruction of charged particles in the angular range $`|\mathrm{cos}\theta |<0.96`$. A charged particle track is called a good track if it is reconstructed with a least four hits in the TPC and if it originated from within a cylinder of length 20 cm and radius 2 cm, coaxial with the beam and centred at the interaction point. A 1.5 T axial magnetic field delivered by a super-conducting solenoidal coil allows a charged particle $`1/\mathrm{p}_{}`$ resolution of $`(6\times 10^45\times 10^3/\mathrm{p}_{})(\mathrm{GeV}/c)^1`$ to be achieved.

The electromagnetic calorimeter (ECAL) is a lead/wire-plane sampling calorimeter covering the angular range $`|\mathrm{cos}\theta |<0.98`$. Anode wire signals provide a measurement of the arrival time of the particles relative to the beam crossing with a resolution better than 15 ns. Cathode pads associated with each wire layer are connected to form projective towers of approximately $`0.9^{}`$ by $`0.9^{}`$ which are read out in three segments in depth. The impact parameter of the photon with respect to the interaction point is estimated from the barycentre of the electromagnetic shower in each segment with a resolution of about 6 cm. A photon candidate is identified using a topological search  for energy deposits in neighbouring electromagnetic calorimeter towers isolated from the extrapolation of any charged particle track to the ECAL. Any photon candidate close to a boundary between ECAL modules or pointing towards an uninstrumented region of the TPC is not considered in the analysis. The energy calibration of the ECAL is obtained from Bhabha events, radiative returns to the Z resonance, $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ and $`\gamma \gamma \mathrm{e}^+\mathrm{e}^{}`$ events. The energy resolution for photons is $`\delta \mathrm{E}/\mathrm{E}=0.25/\sqrt{\mathrm{E}/\mathrm{GeV}}+0.009`$ .

The luminosity monitors (LCAL and SICAL) extend the calorimetric coverage down to small polar angles. The iron return yoke is instrumented with streamer tubes and acts as a hadron calorimeter (HCAL), covering polar angles down to 110 mrad. Surrounding the HCAL are two additional double layers of streamer tubes called muon chambers.

The measurements of the tracking detectors and the calorimeters are combined into objects classified as charged particles, photons and neutral hadrons using the energy flow algorithm described in Ref.. All objects are used to compute the total visible energy $`\mathrm{E}_{\mathrm{vis}}`$ and the missing energy $`\mathrm{E}_{\mathrm{miss}}`$ with a resolution of ($`0.6\sqrt{\mathrm{E}_{\mathrm{vis}}/\mathrm{GeV}}+0.6`$) GeV.

Electron identification is based on the matching between the measured momentum in the tracking system and the energy in the ECAL, the shower profile in the ECAL and the measurement of the specific ionisation energy loss in the TPC. Muons are identified by their characteristic hit patterns in the HCAL and in the muon chambers.

## 3 Event preselection

Signal events are characterized by two isolated, energetic photons, well contained in the apparatus, and in time with the beam crossing.

Photon isolation is ensured by requiring that the total charged energy in a cone of half-angle $`14^{}`$ around the photon direction ($`\mathrm{E}_{14^{}}^\gamma `$) be smaller than 2 GeV, and that the invariant mass between the photon and any charged particle ($`\mathrm{m}_{\gamma ,\mathrm{ch}}`$) be in excess of 1 $`\mathrm{GeV}/c^2`$. Pairs of photons from a $`\pi ^0`$ decay are rejected by the requirement that their invariant mass be larger than $`1\mathrm{GeV}/c^2`$.

Photon centrality is enforced by requiring that the polar angles satisfy $`|\mathrm{cos}\theta _\gamma |<0.9`$. Since the Higgs boson is expected to be produced with a nearly uniform $`\mathrm{cos}\theta `$ distribution, the sum of the two cosines (in absolute value) must be smaller than 1.4, and the cosines of the Higgs boson production and the thrust axis polar angles, $`\theta _{\gamma \gamma }`$ and $`\theta _{\mathrm{thrust}}`$ must be within $`\pm 0.95`$. Events with energy at low polar angles are further vetoed by rejecting events with more than 2 GeV within $`14^{}`$ of the beam axis ($`\mathrm{E}_{14^{}}^{\mathrm{beam}})`$.

If, in a selected event, more than two photons satisfy the above criteria, only the most energetic two photons are considered as originating from a Higgs boson candidate.

Finally, only events in time with the beam crossing are kept. For events with at least two good tracks, a timing to better than 1 ns is ensured by the good track definition. Events with only one good track are rejected. The production time of events with no good tracks is determined as the energy-weighted average of the times reconstructed in all ECAL modules, t0, and is required to be in agreement with the beam crossing time within $`\pm 40`$ ns. For these events, the impact parameter of each identified photon is also required to be less than 25 cm.

The reconstructed particles of events with at least two good tracks are then forced to form four “jets” with the Durham jet clustering algorithm . The consistency of these events with a four-body final-state hypothesis is verified by requiring that each of the four jet energies $`\mathrm{E}_i^{\mathrm{resc}}`$, rescaled to satisfy energy-momentum conservation under the assumption that the jet velocities are perfectly measured, be positive. However, in order to make the reconstructed Higgs boson mass resolution independent of the final state topology, the measured photon energies were used instead of the rescaled energies.

In the following, the neutral electromagnetic energy in a given jet is computed with all neutral objects in the electromagnetic calorimeter, and with the objects in the hadron calorimeter found behind uninstrumented regions of the electromagnetic calorimeter. An electromagnetic jet is defined as a jet with more than $`80\%`$ of electromagnetic energy. In events reconstructed with four jets, i.e., with at least two good tracks, the two jets with the smallest electromagnetic energy fraction are called fermionic jets.

At this level, a signal efficiency of 20 to $`65\%`$ is achieved at all centre-of-mass energies and for any Higgs boson mass above 1 $`\mathrm{GeV}/c^2`$, as estimated with many $`\mathrm{e}^+\mathrm{e}^{}\mathrm{Hf}\overline{\mathrm{f}}`$ simulated event samples produced with the HZHA generator  and processed through the whole detector simulation and event reconstruction chain.

## 4 Event selection

Events passing the above preselection are further classified in the four signal topologies according to their good track multiplicity ($`\mathrm{n}_{\mathrm{ch}}`$).

* Events with no good tracks are classified as $`\mathrm{H}\nu \overline{\nu }`$ candidate events. The main standard model background sources to this final state are (i) $`\mathrm{e}^+\mathrm{e}^{}\nu \overline{\nu }\gamma (\gamma )`$, simulated with the KORALZ package ; and (ii) $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma (\gamma )`$, simulated with the GGG generator . The latter generator does not contain QED contributions of order $`\alpha ^4`$ and above, but they were estimated in Ref.  to be small enough to be considered negligible in the present analysis.
* Events with two good tracks, both positively identified as electrons or muons, are classified as $`\mathrm{H}\mathrm{}^+\mathrm{}^{}`$ (with $`\mathrm{}`$ = $`\mathrm{e}`$ or $`\mu `$) candidate events. The programs BHWIDE  and UNIBAB  on the one hand, and KORALZ on the other, are employed to simulate the main background processes, i.e., $`\mathrm{e}^+\mathrm{e}^{}\mathrm{e}^+\mathrm{e}^{}\gamma (\gamma )`$ and $`\mu ^+\mu ^{}\gamma (\gamma )`$, respectively.
* Other two good track events, and events with up to four good tracks are classified as $`\mathrm{H}\tau ^+\tau ^{}`$ candidate events. Here again, the main background process, $`\mathrm{e}^+\mathrm{e}^{}\tau ^+\tau ^{}\gamma (\gamma )`$ is simulated with KORALZ.
* Finally, events with at least five good tracks are classified as $`\mathrm{Hq}\overline{\mathrm{q}}`$ candidate events. The $`\mathrm{e}^+\mathrm{e}^{}\mathrm{q}\overline{\mathrm{q}}(\gamma )`$ process is simulated with JETSET  for data taken at the Z resonance, and with PYTHIA  and HERWIG  at higher energies. The four-fermion processes WW, ZZ, Zee and $`\mathrm{We}\nu `$ are simulated using PYTHIA .

The selection criteria designed to reduce the contribution of the background processes are summarized in Table 2, and only a brief account is given here.

For each of the topologies, the photons of the $`\mathrm{e}^+\mathrm{e}^{}\mathrm{f}\overline{\mathrm{f}}\gamma (\gamma )`$ background process are emitted by the incoming and the outgoing charged fermions. They are therefore preferentially produced either along the beam axis or along one of the outgoing charged fermion directions. An efficient rejection is achieved by tightening the photon isolation requirements with respect to these directions, by means of the variables introduced in Section 3 and a number of other relevant variables: (i) photon energies $`\mathrm{E}_\gamma `$ and the energy recoiling against the photon pair $`\mathrm{E}_\mathrm{H}^{\mathrm{recoil}}`$, (ii) invariant masses ($`\mathrm{m}_{\gamma ,\mathrm{ch}}`$, $`\mathrm{m}_{\gamma ,\mathrm{f}}`$) and transverse momentum $`\mathrm{p}_{_{\mathrm{ch}}}^\gamma `$ of individual charged particles and fermionic jets with respect to the photons, (iii) transverse momentum of the photon with respect to the closest fermionic jet $`\mathrm{p}__\gamma ^\mathrm{f}`$ and transverse momentum of the fermionic jets with respect to the thrust axis $`\mathrm{p}__\mathrm{f}^{\mathrm{thrust}}`$, (iv) the value $`y_{34}`$ of the Durham $`y_{\mathrm{cut}}`$ transition value between three and four jets and the invariant mass of the two fermionic jets $`\mathrm{m}_{\mathrm{f}\overline{\mathrm{f}}}`$. For the Z peak data, exactly two electomagnetic jets $`\mathrm{N}_{\mathrm{jets}}^{\mathrm{elec}}`$ are required in order to reject events in which both photons are inside the same jet.

Further cuts are made on the cosine of the decay angle of the photons (fermions) in the rest frame of the $`\gamma \gamma `$ ($`\mathrm{f}\overline{\mathrm{f}}`$) system $`\mathrm{cos}\theta _{\gamma \gamma }^{\mathrm{dcy}}`$ ($`\mathrm{cos}\theta _{\mathrm{f}\overline{\mathrm{f}}}^{\mathrm{dcy}}`$). These decay angles are expected to have a flat distribution for signal events, while they are strongly peaked towards small angles for the background. For a high mass Higgs at high energies the signal events are expected to be somewhat spherical in nature. Non-spherical events are rejected by a cut on the sum of the four minimum inter-jet angles $`\theta _{4\mathrm{jets}}^{\mathrm{min}}`$. At high energies the mass recoiling against the photonic system $`\mathrm{m}_{\mathrm{rec}}`$ is required to be consistent with the Z mass. For the $`\mathrm{H}\tau \tau `$ channel at Z peak energies the mass recoiling against the photonic system is required to be consistent with the rescaled mass of the fermionic system $`\mathrm{m}_{\mathrm{ff}}^{\mathrm{resc}}`$. Finally, in the $`\mathrm{H}\nu \overline{\nu }`$ topology, some boost is ensured by requiring that the total produced mass does not saturate the available energy: $`|\mathrm{m}_{\gamma \gamma }+\mathrm{m}_{\mathrm{rec}}\sqrt{s}|>\mathrm{\hspace{0.17em}2}\mathrm{GeV}`$.

The typical efficiencies and expected standard model background for each topology, together with cross-channel contamination, are shown in Table 3.

## 5 Systematic uncertainties

Because no background subtraction is performed to derive the final result, the uncertainty on the background evaluation has not been estimated. The uncertainty on the photon selection efficiency summarized in Table 4 receives a large contribution from the photon identification efficiency. The uncertainties due to photon energy calibration or photon angle resolution affect only the $`\gamma \gamma `$ invariant mass resolution.

To estimate the uncertainty on the photon selection efficiency, the total cross section of the $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma (\gamma )`$ process was measured and compared to its prediction. Two-photon-candidate events are selected as final states with only two identified photons (with the same timing and pointing constraints as those defined in section 3) with a polar angle such that $`|\mathrm{cos}\theta _\gamma |<0.9`$ ($`0.95`$), and an opening angle satisfying $`\mathrm{cos}\alpha _{\gamma \gamma }<0.999`$ ($`0.9999`$), for the Z peak data (at high energy). The $`\mathrm{e}^+\mathrm{e}^{}\mathrm{e}^+\mathrm{e}^{}\gamma \gamma `$ background is rejected by requiring no charged particles in the event and not more than 2 GeV of energy around the beam direction ($`\mathrm{E}_{14^{}}^{\mathrm{beam}}<2GeV`$). The $`\nu \overline{\nu }\gamma \gamma `$ background is reduced to a negligible amount by requiring the invariant mass of the two-photon system to be greater than $`0.75\sqrt{s}`$. The result is displayed in Fig. 1. From this measurement, a conservative relative systematic uncertainty is estimated to be $`3\%`$, including a $`1\%`$ theoretical error. The systematic uncertainty on the photon detection efficiency includes the effects of the cuts on photons used at preselection level as well as that of the integrated luminosity determination, which enters the $`\mathrm{e}^+\mathrm{e}^{}\gamma \gamma (\gamma )`$ cross section measurement.

The effect of the photon energy calibration and angular resolution on the efficiency is estimated with Bhabha events. The relative difference between the values of the electron energy determined independently by the tracking system and by the calorimetry mainly originates from electron Bremsstrahlung. This difference is measured to be $``$ $`1.5\%`$, and it agrees within $`0.2\%`$ with that expected from the simulation. This $`\pm 0.2\%`$ uncertainty is conservatively assigned to the photon energy calibration. It leads to a shift of the measured Higgs boson mass by $`\pm 200`$ $`\mathrm{MeV}/c^2`$, and to a negligible increase of the mass resolution. The uncertainty on the reconstructed photon direction is estimated by comparing the directions determined with the tracking and the calorimetry. It degrades the measured Higgs boson mass resolution by approximately 20 $`\mathrm{MeV}/c^2`$ with a negligible loss of efficiency. These systematic uncertainties are included in the final result by increasing the Higgs boson mass resolution by $`200\mathrm{MeV}/c^2`$ .

As the lepton identification is used only to classify events, the signal efficiency decreases by less than $`0.1\%`$ relative when the lepton identification efficiency is modified by $`10\%`$.

The remaining selection criteria are expected to be largely insensitive to the details of the simulation of the hadronic system. A quantitative estimate of the size of possible discrepancies is performed with the aid of an event reweighting technique. For each selection variable, bin-by-bin correction factors are calculated as the ratio of data to Monte Carlo expectation, evaluated at the preselection level with the cut on $`\mathrm{m}_{\gamma ,\mathrm{ch}}`$ removed. The Monte Carlo signal distribution of the selection variable, obtained when all cuts are applied, is then re-weighted with these correction factors to obtain a new estimate of the efficiency. The difference in efficiencies estimated with this technique are, for almost all variables, at the level of a few parts per mil and are given in Table 4. The largest effect, of $`1.2\%`$, comes from the cut on the isolation of the photon $`\mathrm{m}_{\gamma ,\mathrm{ch}}`$.

Finally, the model dependence of this analysis is estimated with the anomalous coupling model mentioned in Section 1. In that context, all the anomalous coupling parameters $`\mathrm{f}_i/\mathrm{\Lambda }^2`$ are varied independently of each other within $`\pm 100`$ $`\mathrm{TeV}^2`$ using the HZHA generator. The influence of the Higgs boson energy and angular distributions leads to a relative uncertainty on the signal efficiency of at most $`4\%`$.

The decay width of the Higgs boson is in general negligible with respect to the detector resolution (below 1 $`\mathrm{GeV}`$ for $`|\mathrm{f}_i/\mathrm{\Lambda }^2|<100`$ $`\mathrm{TeV}^2`$). In anomalous coupling models with $`|\mathrm{f}_i/\mathrm{\Lambda }^2|`$ above 500 $`\mathrm{TeV}^2`$ the decay width of a heavy Higgs boson becomes larger than the detector invariant mass resolution ($`\mathrm{\Gamma }3`$ $`\mathrm{GeV}`$) and the relative variation of the signal efficiency becomes larger than the previously estimated systematic uncertainties. The final results presented in the next section are therefore valid for models in which the width of the Higgs boson is less than a few GeV.

## 6 Results

The signal efficiencies are displayed in Fig. 2a as a function of the Higgs boson mass hypothesis. Twenty-two events are selected in the data compared with an expectation of $`27.9\pm 1.6`$ events from all standard model background processes. The contributions of the various channels are given in Table 3. The diphoton invariant mass distribution of these selected events is shown for both data and Monte Carlo expectation in Fig. 2b. The low mass background is dominated by events where both photons originated from the same jet, whereas the high mass background is dominated by events where the two photons originated from different jets. The sum of these two contributions has a minimum at around $`8\mathrm{GeV}/c^2`$. No evidence of a resonance decaying to $`\gamma \gamma `$ is observed and a $`95\%`$ confidence level upper limit on the number of signal events at a given diphoton invariant mass is derived following the method described in Ref. . For each LEP energy, the total efficiency and the mass resolution are parametrized as a function of the $`\gamma \gamma `$ invariant mass. The diphoton invariant mass resolution varies linearly from $`1\mathrm{GeV}/c^2`$ to $`3.5\mathrm{GeV}/c^2`$ over the whole Higgs boson mass range. The systematic uncertainties are conservatively taken into account by scaling down the efficiency by $`5.2\%`$ and increasing the mass resolution as described in the previous section.

The present analysis does not apply for Higgs boson masses below 1 $`\mathrm{GeV}/c^2`$. To extend the search to smaller masses, the direct measurement of the Z invisible width with single photon counting described in Ref.  was used. For this measurement, the single photon candidate events were selected as final states with only one cluster in the electromagnetic calorimeter. The $`\mathrm{H}\nu \overline{\nu }`$ events with $`\mathrm{m}_\mathrm{H}<1\mathrm{GeV}/c^2`$ would therefore have been selected with a good efficiency, ranging from 13 to $`45\%`$, thus leading to a sensitivity similar to that of the present analysis.

Assuming a Higgs boson production cross section with the same $`\sqrt{s},\mathrm{m}_\mathrm{H}`$ dependencies as in the standard model, the $`95\%`$ confidence level upper limit on the product branching ratio $`\mathrm{B}(\mathrm{H}\gamma \gamma )\sigma (\mathrm{e}^+\mathrm{e}^{}\mathrm{Hf}\overline{\mathrm{f}})/\sigma ^{\mathrm{SM}}(\mathrm{e}^+\mathrm{e}^{}\mathrm{Hf}\overline{\mathrm{f}})`$ is derived and shown in Fig. 3.

For the case of a Higgs boson produced at the standard model rate, the best upper limit on the branching ratio (4.7$`\times 10^3`$ at $`95\%`$ confidence level) is obtained for Higgs boson masses below 20 $`\mathrm{GeV}/c^2`$. A Higgs boson decaying exclusively to two photons is ruled out up to 109 $`\mathrm{GeV}/c^2`$ at $`95\%`$ confidence level.

A fermiophobic Higgs boson  with no tree-level coupling to fermions is excluded at $`95\%`$ confidence level for any mass up to 100.7 $`\mathrm{GeV}/c^2`$.

The present analysis extends the reach of similar analyses performed by other LEP collaborations .

## 7 Conclusion

With a data sample of 672 $`\mathrm{pb}^1`$ recorded at centre-of-mass energies from 88 GeV to 202 GeV, a search for two photon decays of Higgs bosons produced in association with a fermion pair has been performed in the mass range from 0 up to 120 $`\mathrm{GeV}/c^2`$. No evidence for resonant production of photon pairs has been found. In the framework of Higgs bosons with standard model coupling to gauge bosons, a $`95\%`$ confidence level upper limit on the Higgs boson branching ratio to two photons has been obtained for Higgs boson masses from 0 to 109 $`\mathrm{GeV}/c^2`$. In the fermiophobic model in which the Higgs boson couples exclusively to gauge bosons, a $`95\%`$ confidence level lower limit on the Higgs boson mass has been set at 100.7 $`\mathrm{GeV}/c^2`$.

## 8 Acknowledgements

We wish to thank our colleagues from the CERN accelerator divisions for the successful operation of LEP. We are indebted to the engineers and technicians in all our institutions for their contribution to the excellent performance of ALEPH. Those of us from non-member countries thank CERN for its hospitality.


format: text

text:
# High redshift radio galaxies: Beacons to biased hierarchical galaxy formation within large scale structure

## 1. Some History

Prior to the mid-1990’s, the only high redshift galaxies known were the parent galaxies of powerful radio sources. Starting with Cygnus A at $`z=0.056`$ in the 1950’s, and moving to 3C 295 at $`z=0.461`$ in the 1960’s, then making the jump to 4C 41.17 at $`z=3.8`$ in the 1980’s, powerful radio galaxies were, by far, the highest redshift galaxies known. As such, radio galaxies have provided a unique, although admittedly biased, view of primeval galaxies and galaxy evolution (McCarthy 1993). The last few years has seen the discovery of populations of ‘normal’ galaxies at high redshift, including optically selected UV-dropout galaxies, and submm selected dusty starburst galaxies (Steidel et al. 1998, Hughes et al. 1998). At first glance, these discoveries apparently diminish the importance of high redshift radio galaxies in the general study of primeval galaxies. More correctly, we feel these discoveries have placed radio galaxies in their proper context among high-z galaxy populations. High redshift radio galaxies remain the most (optically) luminous of the known high-z galaxies, and they have extended, high surface brightness line and continuum emission from radio through X-ray wavelengths. Likewise, low frequency radio surveys provide large, relatively unbiased samples. Hence, powerful radio galaxies are still the most easily and best studied of the high-z galaxy populations, and extensive observations across the electromagnetic spectrum have provided information that is fundamental to our understanding of galaxy and large scale structure formation.

While most authors agree that the parent galaxies for luminous Fanaroff-Riley Class II radio galaxies (P$`{}_{178\mathrm{M}\mathrm{H}\mathrm{z}}{}^{}`$ 1035 h-2 erg s-1 Hz-1)111h $`\frac{\mathrm{H}_\mathrm{o}}{100}`$. We use Ho = 50 km s-1 Mpc-1 and qo = 0.5, unless stated otherwise., are giant elliptical galaxies, at least to $`z1.5`$, there has been debate in the recent literature concerning the evolution of the cluster environments of powerful radio galaxies from $`z=0`$ to $`z1.5`$. Hill and Lilly (1991) concluded that there is substantial evolution in the environments of powerful radio galaxies over this redshift range, from small groups at low-z, to Abell class 0 to 1 clusters at $`z>0.5`$. The recent study by Best (2000) of intermediate-z 3C sources seems to support this conclusion, and he uses this fact to argue that intermediate-z radio galaxies cannot be the precusors of the lower redshift sources. McLure and Dunlop (2000) have recently challenged these findings. They find low redshift luminous AGN, both radio loud and radio quiet, typically inhabit Abell class 0 clusters, and hence that the mean environment for powerful radio galaxies doesn’t change dramatically from $`z=0\mathrm{to}1`$. They claim that: ‘powerful AGN do not avoid rich clusters, but rather display a spread in cluster environment which is consistent with being drawn at random from the massive elliptical galaxy population.’ They make the point that the co-moving space density of clusters at low redshift, $`5\times 10^5`$ h3 Mpc-3, is comparable to that of powerful AGN at $`z=2.5`$, with the interesting implication that all cluster dominant elliptical galaxies were active at $`z2.5`$ (West 1994).

The study of the low-to-intermediate redshift evolution of luminous radio galaxies is ham-strung by the paucity of such sources at low-z. Figure 1 shows rest-frame P178MHz versus redshift for a number of samples of radio galaxies. The 3C samples shows the clear Malmquist bias inherent in flux limited samples. Subsequent studies of lower flux density sources have filled-in the distribution at higher redshifts by revealing lower-luminosity sources. However, there is a clear lack of high-power sources at $`z0.5`$. This is simply a consequence of our universe, and cannot be overcome with larger samples.

## 2. The Leiden Sample

Through pain-staking optical identification of sources from low frequency radio samples, a number of groups, most notably those at Leiden, Oxford, and Berkeley, have expanded the size and redshift range over which radio galaxies have been detected. The current record is $`z=5.2`$ (van Breugel et al. 1999a), and there are 150 sources known at $`z>2`$. Many authors have called attention to the interesting parallel between the rapid increase with redshift in the space density of powerful radio galaxies and other luminous AGN, with that of the cosmic star formation rate. Both show more than an order of magnitude increase in co-moving space density from $`z=0`$ to 2 (Dunlop et al. 2000, Boyle and Terlevich 1998, Blain et al. 1999).

Table 1 gives a broad-brush summary of the properties of powerful radio galaxies at $`z1`$ compared with those at $`z>2`$. While all of the statements are defensible, some are based on observations of just a few sources (or even one!), and remain to be verified.

In the near-IR, the host galaxies change from well defined giant elliptical galaxies at low-z with a fairly tight $`Kz`$ relation, to spatially more complex systems at high-z, with a larger scatter in the $`Kz`$ (van Breugel et al. 1999b). Radio galaxies at $`z1`$ also show evidence for old stars, suggesting a formation redshift $`z>3`$ (Dunlop et al. 2000). The blue light shows a general alignment with the radio source axis at all redshifts, but high redshift sources also show blue knots of emission well away from the radio source axis. Also, at low-z the blue emission is typically highly polarized, while at high-z the emission is less polarized, or even unpolarized, and one source shows stellar absorption features (Dey et al. 1997). The implication is that at low-z the blue emission is dominated by scattered light from the obscured nucleus, while at high-z star forming regions make a significant contribution. Diffuse Ly-$`\alpha `$ emission is seen at both low and high-z, but the Ly-$`\alpha `$ halos are typically larger, and more luminous at high-z (Baum and McCarthy 2000). Self-absorption of the Ly-$`\alpha `$ emission is also more prevalent in high redshift sources, indicating large amounts of cooler HI (van Ojik et al. 1997, Debreuck et al. 2000).

The radio source morphologies change from small, irregular systems at high-z to large classical doubles at low-z. Extreme values of Faraday rotation of the polarized radio emission are observed at both low-z and high-z. Pentericci (1999) has argued that the fraction of sources with large rotation measures (RMs) increases with redshift. The origin of large RMs in low-z systems is thought to be magnetized cluster atmospheres, with fields $`10\mu `$G ordered on scales of 10’s of kpc. The location of the Faraday screen in high-z sources remains an open question.

## 3. Dust emission

Powerful radio galaxies at low-z typically have IR luminosities $`<10^{12}`$ L year-1, some of which can be ascribed to synchrotron radiation from the AGN (Golombek et al. 1988). High-z radio galaxies show dramatically different behavior in the rest-frame IR.

Archibald et al. (1999) have shown a rapid evolution with redshift in the rest-frame IR luminosity of radio galaxies, with the median value increasing as $`(1+z)^{3.5}`$ to $`10^{13}`$ L year-1 at $`z=4`$. The spectral energy distributions show that the emission is from warm dust, with implied dust masses $`10^8`$ M, and dust temperatures between 50 and 100 K (Benford et al. 2000). Emission from CO has also been detected from a few high-z radio galaxies, with implied gas masses $`10^{10}`$ M (Papadopoulos et al. 2000).

Perhaps more importantly, imaging at submm and mm wavelengths has shown that the peak of the dust and CO emission is often not centered on the AGN, but can arise in sources distributed at distances of 1 or more from the AGN. The most dramatic example of this is 4C 41.17 at $`z=3.80`$ recently imaged with SCUBA at 850$`\mu `$m by Ivison et al. (2000). They detect six luminous dust emitting sources in an area of about 4 square arcmin centered on the radio galaxy, including emission from the AGN itself. Interestingly, the brightest dust emitting source is located at a distance of 1 from the AGN. This represents an order-of-magnitude over-density of submm sources relative to the field.

A somewhat less dramatic, but equally telling, example of this phenomenon is shown in Figure 2. Dust emission has been detected from the $`z=2.156`$ radio galaxy PKS 1138–262 at 240 GHz using MAMBO (Kreysa et al. 1999) at the IRAM 30m telescope. A series of pointed and imaging observations shows that the dust emission may not peak at the position of the AGN, but peaks about 4′′ to the north-east of the AGN, at the location of a complex of Ly-$`\alpha `$ emitting knots designated BC in Figure 2. The implied star formation rate for this complex, based on the mm emission, is a few hundred M year-1. This is two orders-of-magnitude greater than the star formation rate predicted from the rest-frame UV luminosity (Pentericci et al. 1998), and is characteristic of the recently discovered dust obscured, high-z starburst galaxy population (Rowan-Robinson 1999).

## 4. A proto-cluster around PKS 1138–262

The $`z=2.156`$ radio galaxy PKS 1138–262 was chosen as the target of a pilot project to search for galaxy over-densities around high redshift radio galaxies for a number of reasons: (i) its extremely clumpy optical morphology and large size that resembles simulations of forming massive cluster galaxies (Pentericci et al. 1998), (ii) its extreme RM and distorted radio morphology indicating dense local gas (Carilli et al. 1997), and (iii) detection of extended X-ray emission (Carilli et al. 1998).

Narrow band imaging with the VLT revealed 50 candidate Ly-$`\alpha `$ emitting galaxies within a $`7^{}`$ field centered on PKS1138–262 and in the redshift range of 2.137$`\pm `$0.037 (Kurk et al. 2000), not including the 10 previously known Ly-$`\alpha `$ emitting knots within 8′′ of the radio galaxy itself (Pentericci et al 2000). Follow-up spectroscopy of 43 candidates showed 15 of the sources in the range z = 2.155$`\pm `$0.016. There is no obvious angular sub-structure in the distribution of the Ly-$`\alpha `$ emitting sources, but there is clear redshift sub-structure, with about half the sources being at $`z=2.143\mathrm{to}2.147`$, and the other half at $`z=2.162\mathrm{to}2.166`$. The velocity dispersion of each sub-group is about 400 km s-1, implying a virial mass of $`7\times 10^{13}`$ M for each group.

The over-density of galaxies in the 1138–262 field is comparable to that seen by Steidel et al. (1998) in their UV-dropout selected proto-cluster at $`z=3.090`$ (see the analysis in Pentericci et al. 2000). In terms of large scale structure formation in the context of CDM models with standard cosmologies, Steidel et al. argue that the space density of such structures at high redshift requires either a low value of $`\mathrm{\Omega }_M`$ and/or a large bias parameter, $`b`$, with $`b=2`$ for $`\mathrm{\Omega }_M=0.2`$, or $`b=6`$ for $`\mathrm{\Omega }_M=1`$.

## 5. X-rays from PKS 1138–262

At low-z there is a clear correlation between extreme RM radio galaxies and X-ray emitting cluster atmospheres (Taylor et al. 1994). Based on this correlation, we have begun a search for extended X-ray emission from high-z radio galaxies with extreme RM values. Our pilot study has targeted PKS 1138-262. A long exposure with ROSAT HRI detected X-ray emission from PKS 1138–262 with a rest frame 2 to 10 keV luminosity of $`7\times 10^{44}`$ erg s-1, comparable to that of a $`10^{14}`$ M cluster atmosphere. Unfortunately, the spatial resolution of ROSAT was insufficient to determine the spatial distribution of the X-ray emission, at least on scales comparable to that expected for high surface brightness cluster core emission $`100`$ kpc.

Figure 3 shows the X-ray emission from 1138–262 as seen in a 40 ksec exposure with Chandra. The total X-ray emission (0.2 to 10 keV) shows a point source at the position of the AGN, plus extended emission on a scale of about 20′′. The extended X-ray emission has a softer spectrum than the point source, contributing about 10$`\%`$ of the total luminosity in the 0.2 to 10 keV band, and 25$`\%`$ in the 0.2 to 1 keV band.

The soft X-ray emission (0.2 to 1 keV) shows a general correlation with the radio source in both extent and orientation. However, the correspondence is not exact in terms of individual high surface brightness features, most notably, the eastern radio hot spot is situated about 2′′ further from the nucleus than the local peak in the X-ray emission. Also, the X-ray emission has a broader (transverse) distribution than the radio emission. Note that the relative astrometry is only accurate to about 2′′. However, shifting the images in order to line-up the eastern hot spot would then move the X-ray nucleus and western extension away from the corresponding radio structures. Our analysis of the X-ray emission from 1138–262 has just begun, so the following conclusions are very tentative.

We feel it unlikely that the X-ray emission is synchrotron radiation, given the lack of a one-to-one correspondence with radio features, and the fact that extrapolating the radio spectrum to the X-ray under-predicts the X-ray flux density by more than an order of magnitude. Inverse Compton emission is also unlikely, for similar spatial reasons, and the derived magnetic field is an order of magnitude below the equipartition value. The one caveat is that, if the radio source morphology is significantly different at low radio frequencies relative to the high (rest-frame) frequencies observed thus far, then inverse Compton emission might become a viable option.

One possible explanation is thermal emission from hot gas, which would be consistent with the large RMs observed. The broad-band spectrum is consistent with gas at 5 to 10 keV, although the spectrum varies with position. The correlation with the radio source axis would then suggest gas heated by shocks driven by the expanding radio source. In this case, the lack of emission from a virialized cluster atmosphere might indicate that the ambient medium was at a relatively low temperature ($``$ 1 keV) prior to the on-set of the radio jet. The observation of stronger Ly-$`\alpha `$ emission, and more prevalent Ly-$`\alpha `$ absorption, in high redshift radio galaxies suggests a larger supply of cooler ambient material on scales of 10 to 100 kpc relative to low redshift radio galaxies (van Ojik et al. 1997 Debreuck et al. 2000). We are currently considering the energetics and gas mass involved, to determine if such a model is sensible.

## 6. High redshift radio galaxies: Beacons to biased hierarchical galaxy formation within large scale structure

Based on the observations presented above, we speculate that high redshift radio galaxies are often beacons to highly biased hierarchical galaxy formation within large scale structure (LSS). The LSS will eventually separate from the Hubble flow and collapse to form a rich cluster at $`z1`$. The individual Ly-$`\alpha `$ and submm-selected galaxies will evolve into cluster spheroidal galaxies, including the giant elliptical AGN host galaxy.

Two question raised during the meeting are worthy of mention. One question by Linda Sparkes concerns the energetics: given the 1045 erg s-1 coming from the radio AGN in these systems, could it be that all observational characteristics are dictated by the AGN? While this could be the case for regions within 50 kpc or so of the AGN, it is clearly not the case for most of the submm and Ly-$`\alpha `$ emitting galaxies distributed at distances of arcminutes from the AGN, and, for the most part, not along the radio axis. Also, the implied IR luminosities for the spatially distributed submm sources are comparable to, or larger than, the integrated luminosity of the radio galaxy. Hence, while the extended radio source may appear dramatic on high quality radio images, the bolometric luminosity of the ‘cluster’ as a whole may not be dominated by the radio galaxy. This raises the interesting possibility that the host galaxies of high redshift radio sources may not be pre-ordained dominant cluster galaxies in all cases, but could become normal cluster ellipticals, thereby alleviating the necessity for all cluster dominant galaxies to host AGN at $`z2.5`$ (McClure and Dunlop 2000).

The second question by James Binney concerns timescales. There is a clear correlation with radio source size and redshift, with the median source size dropping from 100 kpc at low-z to 10 kpc at $`z=4`$ (Carilli et al. 1998). The most recent analysis of this phenomenon suggests that it is due to a correlation between radio source size, age, and luminosity: smaller sources are younger and more luminous (Blundell et al. 1999). As with the energtics discussed above, this age discrimination cannot dictate the characteristics of the submm and Ly-$`\alpha `$ sources distributed throughout the (proto-) cluster. Also, timescales for the radio activity are typically 106 to $`10^7`$ years, while those for the star forming galaxies are likely to be an order of magnitude longer.

While high-z radio galaxies are identified via their radio properties, ie. steep spectrum radio sources, current data indicate that powerful radio galaxies and their environments represent ideal laboratories to study the general phenomena of galaxy and large scale structure formation over the widest possible range in redshifts. In essence, the radio source acts like a bolt of lightening illuminating a storm – a dramatic, short timescale event that draws attention to a larger, even more dramatic phenomenon, namely the formation of a cluster of galaxies.

### Acknowledgments.

The National Radio Astronomy Observatory (NRAO) is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc.. We thank Carlos Debreuck for his tireless work on the identifications and spectroscopic redshift measurements of many of the high redshift radio galaxies plotted in Figure 1. The work by W.v.B. at the University of California Lawrence Livermore National Laboratory was performed under the auspices of the US Department of Energy under contract W-7405-ENG-48.

## References

Archibald, E.N. et al. 2000, MNRAS, in press

Baum, S. and McCarthy, P. 2000, AJ, in press

Best, P.N. 2000, MNRAS, in press

Benford, D. et al. 2000, ApJ, 518, 65

Blain, A.W. et al. 1999, MNRAS, 309, 705

Blundell, K.M. and Rawlings, S. 1999, Nature, 339, 330

Boyle, B.J. and Terlevich, R.J. 1998, MNRAS, 293, L49

van Breugel, W. et al. 1999a, in ’The Most Distant Radio Galaxies,’ eds. Rottgering, Best, and Lehnert (KNAW: Amsterdam), p. 49

van Breugel, W. et al. 1999b, ApJ, 518, L61

Carilli, C.L. et al. 1997, ApJ(supp), 109, 1

Carilli, C.L. et al. 1998, ApJ, 496, L57

Debreuck, C. et al. 2000, ApJ, in press

Dey, A. et al. 1997, ApJ, 490, 698

Dunlop, J.S. et al. 2000, MNRAS, in press

Golombek, D., Miley, G.K., and Neugebauer, G. 1988, 95, 26

Hill, G.J. and Lilly, S.J. 1991, ApJ, 367, 1

Hughes, D. et al. 1998, Nature, 394, 241

Ivison, R.J., et al. 2000, ApJ, in press

Kreysa, E. et al. 1999, SPIE, 3357, 319

Kurk, J. et al. 2000, å, 358, L1

McCarthy, P. 1993, ARA&A, 31, 639

McClure, R.J. and Dunlop, J.S. 2000, MNRAS, in press

van Ojik, R. et al. 1997, å, 317, 358

Papadopoulos, P. et al. 2000, ApJ, 528, 626

Pentericci, L. et al. 2000, å(letters), in press

Pentericci, L. et al. 1998, ApJ, 504, 139

Rowan-Robinson, M. 2000, MNRAS, in press

Steidel, C.C. et al. 1998, ApJ, 492, 428

Taylor, G.B., Barton, E. and Ge, J.-P. 1994, AJ, 107, 1942

West, M. 1994, MNRAS, 268, 79


text:
# NRQCD Analysis of Bottomonium Production at the Tevatron

## I Introduction

The NRQCD factorization formalism provides a systematic framework for analyzing the inclusive production of heavy quarkonium . Long-distance effects involving the binding of a heavy quark-antiquark pair into quarkonium are factored into parameters called NRQCD matrix elements. These nonperturbative parameters are universal, so values extracted from one high energy physics experiment can be used to predict the production rate in others. The NRQCD matrix elements scale as definite powers of $`v`$, where $`v`$ is the typical relative velocity of the heavy quark. The NRQCD factorization approach becomes phenomenologically useful upon truncating the expansion in $`v`$ so as to reduce the independent NRQCD matrix elements to a manageable number. The truncation is most reliable for the heaviest quarkonium states, namely the $`b\overline{b}`$ system for which $`v^2`$ is roughly 1/10.

The most abundant source of data on bottomonium production is the Tevatron $`p\overline{p}`$ collider. In Run IA of the Tevatron, the CDF collaboration was able to resolve the individual S-wave bottomonium states $`\mathrm{{\rm Y}}(1S)`$, $`\mathrm{{\rm Y}}(2S)`$, and $`\mathrm{{\rm Y}}(3S)`$ and measure their production cross sections . An analysis of the CDF data within the NRQCD factorization formalism was carried out by Cho and Leibovich . The analysis is complicated by the production of P-wave bottomonium states that subsequently make transitions to S-wave states. Cho and Leibovich found that the CDF data was insufficient to determine all the important NRQCD matrix elements and they had to make educated guesses for some of them.

The CDF collaboration has recently analyzed the data on bottomonium production from Run IB at the Tevatron. In addition to much higher statistics on the cross sections for $`\mathrm{{\rm Y}}(1S)`$, $`\mathrm{{\rm Y}}(2S)`$, and $`\mathrm{{\rm Y}}(3S)`$ , they also have results on the production of the P-wave states $`\chi _b(1P)`$ and $`\chi _b(2P)`$ . The high quality of the new CDF data justifies an updated theoretical analysis, with careful attention to the experimental and theoretical errors.

In this paper, we present a quantitative analysis of the new CDF data on bottomonium production within the NRQCD factorization formalism. The color-singlet NRQCD matrix elements for S-wave states are determined from their electromagnetic decays, while those for P-wave states are estimated from potential models. The color-octet NRQCD matrix elements are determined by fitting the CDF data, taking full account of the feeddown from transitions between bottomonium states. The resulting values of the matrix elements are used to predict the cross sections for the spin-triplet and spin-singlet bottomonium states in Run II of the Tevatron.

## II NRQCD Matrix Elements

The NRQCD factorization approach provides a model-independent framework for analyzing the inclusive production of heavy quarkonium . The factorization formula for the differential cross section for the inclusive production of a bottomonium state $`H`$ of momentum $`P`$ has the schematic form

$$d\sigma [H(P)]=\underset{n}{}d\sigma [b\overline{b}(n,P)]O^H(n),$$
(1)
where the sum extends over both color-singlet and color-octet and over all angular momentum channels for the $`b\overline{b}`$ pair. The $`b\overline{b}`$ cross sections, which are independent of the bottomonium state $`H`$, can be calculated using perturbative QCD. All dependence on the state $`H`$ is factored into parameters $`O^H(n)`$ called NRQCD matrix elements. These phenomenological parameters can be expressed as matrix elements in an effective field theory called nonrelativistic QCD (NRQCD). A nonperturbative analysis of NRQCD reveals how the various matrix elements scale with the typical relative velocity $`v`$ of the heavy quark in quarkonium. Spin symmetry, which is an approximate symmetry of QCD, also gives relations between various matrix elements.

The relative importance of the terms in the factorization formula (1) depends on the size of the $`b\overline{b}`$ cross sections and on the size of the matrix elements. According to the velocity-scaling rules, the most important matrix element for direct $`\mathrm{{\rm Y}}(1S)`$ production is the color-singlet parameter $`O_1^{\mathrm{{\rm Y}}(1S)}(^3S_1)`$. The spin-symmetry relations can be used to reduce the next most important matrix elements to three color-octet parameters: $`O_8^{\mathrm{{\rm Y}}(1S)}(^3S_1)`$, $`O_8^{\mathrm{{\rm Y}}(1S)}(^1S_0)`$, and $`O_8^{\mathrm{{\rm Y}}(1S)}(^3P_0)`$. These color-octet matrix elements are important, because the cross sections for producing color-octet $`b\overline{b}`$ pairs can be much larger than for color-singlet $`b\overline{b}`$ pairs. There are analogous matrix elements that describe the direct production of $`\mathrm{{\rm Y}}(2S)`$ and $`\mathrm{{\rm Y}}(3S)`$. The NRQCD factorization formula (1) for direct $`\mathrm{{\rm Y}}(nS)`$ production reduces to

$`d\sigma [\mathrm{{\rm Y}}(nS)]`$ $`=`$ $`d\sigma [b\overline{b}_1(^3S_1)]O_1^{\mathrm{{\rm Y}}(nS)}(^3S_1)+d\sigma [b\overline{b}_8(^3S_1)]O_8^{\mathrm{{\rm Y}}(nS)}(^3S_1)`$ (3)
$`+d\sigma [b\overline{b}_8(^1S_0)]O_8^{\mathrm{{\rm Y}}(nS)}(^1S_0)+\left({\displaystyle \underset{J}{}}(2J+1)d\sigma [b\overline{b}_8(^3P_J)]\right)O_8^{\mathrm{{\rm Y}}(nS)}(^3P_0).`$
The factor of $`2J+1`$ in the last term comes from using a spin-symmetry relation to eliminate $`O_8^{\mathrm{{\rm Y}}(nS)}(^3P_J)`$ in favor of $`O_8^{\mathrm{{\rm Y}}(nS)}(^3P_0)`$.

The most important matrix elements for the direct production of the P-wave states $`\chi _{bJ}(1P)`$, $`J=0,1,2`$, can be reduced to a color-singlet parameter $`O_1^{\chi _{b0}(1P)}(^3P_0)`$ and a single color-octet parameter $`O_8^{\chi _{b0}(1P)}(^3S_1)`$. There are analogous matrix elements that describe the direct production of $`\chi _{bJ}(2P)`$ and $`\chi _{bJ}(3P)`$. The NRQCD factorization formula (1) for direct $`\chi _{bJ}(nP)`$ production reduces to

$`d\sigma [\chi _{bJ}(nP)]`$ $`=`$ $`d\sigma [b\overline{b}_1(^3P_J)]O_1^{\chi _{bJ}(nP)}(^3P_J)+(2J+1)d\sigma [b\overline{b}_8(^3S_1)]O_8^{\chi _{b0}(nP)}(^3S_1).`$ (4)
In the last term, the factor of $`2J+1`$ comes from using a spin-symmetry relation to eliminate $`O_8^{\chi _{bJ}(nP)}(^3S_1)`$ in favor of $`O_8^{\chi _{b0}(nP)}(^3S_1)`$. We can also use a spin-symmetry relation to replace $`O_1^{\chi _{bJ}(nP)}(^3P_J)`$ in the first term by $`(2J+1)O_1^{\chi _{b0}(nP)}(^3P_0)`$. The matrix elements for $`\mathrm{{\rm Y}}(nS)`$ and $`\chi _{bJ}(nP)`$ enumerated above should be sufficient for a quantitative description of the production of S-wave and P-wave bottomonium states.

The NRQCD factorization formula gives the cross section for the direct production of a given bottomonium state. The cross sections that are most easily measured in experiments are inclusive cross sections that include contributions from the direct production of higher bottomonium states which subsequently decay into the given state. For example, the feeddown from $`\chi _b(1P)`$, $`\mathrm{{\rm Y}}(2S)`$, and $`\chi _b(2P)`$ accounts for roughly 27%, 11%, and 11% of the $`\mathrm{{\rm Y}}(1S)`$ cross section, respectively . Taking into account the feeddown from higher $`\mathrm{{\rm Y}}(mS)`$ and $`\chi _{bJ}(mP)`$ states, the cross section for inclusive $`\mathrm{{\rm Y}}(nS)`$ production can be written

$`d\sigma [\mathrm{{\rm Y}}(nS)]_{\mathrm{inc}}`$ $`=`$ $`d\sigma [b\overline{b}_1(^3S_1)]O_1(^3S_1)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(nS)}+{\displaystyle \underset{J}{}}d\sigma [b\overline{b}_1(^3P_J)]O_1(^3P_J)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(nS)}`$ (7)
$`+d\sigma [b\overline{b}_8(^3S_1)]O_8(^3S_1)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(nS)}+d\sigma [b\overline{b}_8(^1S_0)]O_8(^1S_0)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(nS)}`$
$`+\left({\displaystyle \underset{J}{}}(2J+1)d\sigma [b\overline{b}_8(^3P_J)]\right)O_8(^3P_0)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(nS)},`$
where the “inclusive NRQCD matrix elements” are

$$O[n]_{\mathrm{inc}}^{\mathrm{{\rm Y}}(nS)}=\underset{H}{}B_{H\mathrm{{\rm Y}}(nS)}O^H[n].$$
(8)
The sum over $`H`$ includes $`\mathrm{{\rm Y}}(nS)`$ and all higher bottomonium states that can make transitions to $`\mathrm{{\rm Y}}(nS)`$. The coefficient $`B_{HH^{}}`$ is the inclusive branching fraction for $`H`$ to decay into $`H^{}`$. By convention, we define $`B_{HH}=1`$. The inclusive branching fraction for the observed bottomonium states are collected in Table I. These numbers were obtained by combining the measured branching fractions for the exclusive decays $`\mathrm{{\rm Y}}(nS)\chi _{bJ}(mP)+\gamma `$, $`\chi _{bJ}(nP)\mathrm{{\rm Y}}(mS)+\gamma `$, and $`\mathrm{{\rm Y}}(nS)\mathrm{{\rm Y}}(mS)+\pi \pi `$, with the exception of $`B_{\mathrm{{\rm Y}}(3S)\mathrm{{\rm Y}}(1S)}`$, which is a direct measurement .

In Table I, we have not included the spin-singlet states $`\eta _b(nS)`$ and $`h_b(nP)`$, which have yet to be observed. Transitions between spin-singlet and spin-triplet states are suppressed, because they proceed through magnetic $`\mathrm{\Delta }S=1`$ transitions. The rates for $`\mathrm{\Delta }S=1`$ transitions are suppressed relative to those for $`\mathrm{\Delta }S=0`$ transitions by a factor of $`v^2`$, which is roughly an order of magnitude. The branching fractions for $`\eta _b(2S)`$ into other bottomonium states are further suppressed by its large annihilation width into two gluons. Quantitative estimates of the electromagnetic and hadronic transition rates are given in Refs. . They support the conclusion that the branching fractions for decays of spin-singlet states into spin-triplet states can be neglected.

The tiny branching fractions in Table I for the transition $`\chi _{bJ}(2P)\chi _{bJ^{}}(1P)`$ are the contributions from the double radiative transitions via $`\mathrm{{\rm Y}}(2S)`$. We have not included the contributions from the two-pion decays $`\chi _{bJ}(2P)\chi _{bJ^{}}(1P)+\pi \pi `$, which have not been observed. We can estimate their magnitude by observing that the rates for $`\mathrm{{\rm Y}}(3S)\mathrm{{\rm Y}}(2S)+\pi \pi `$ and $`\mathrm{{\rm Y}}(3S)\mathrm{{\rm Y}}(2S)+\gamma \gamma `$ are equal to within experimental errors. Since the phase space available for the transitions $`\chi _{bJ}(2P)\chi _{bJ^{}}(1P)`$ is similar to that for $`\mathrm{{\rm Y}}(3S)\mathrm{{\rm Y}}(2S)`$, we expect the rate for $`\chi _{bJ}(2P)\chi _{bJ^{}}(1P)+\pi \pi `$ to be comparable to that for $`\chi _{bJ}(2P)\chi _{bJ^{}}(1P)+\gamma \gamma `$. Including the effects of two-pion transitions could increase the branching fraction for $`\chi _{bJ}(2P)\chi _{bJ^{}}(1P)`$ by a factor of 2 or 3, but since the values of $`B_{\chi _{bJ}(2P)\chi _{bJ^{}}(1P)}`$ in Table I are all less than 1.5%, they should still be negligible.

As indicated by the entries “0 ?” in the first row of Table I, we neglect the feeddown from the $`\chi _b(3P)`$ states, which have not been observed. A naive extrapolation from the other entries of the Table suggest that the branching fractions for $`\chi _{b1}(3P)`$ and $`\chi _{b2}(3P)`$ into $`\mathrm{{\rm Y}}(3S)`$ could be about 12%, while their branching fractions into $`\mathrm{{\rm Y}}(1S)`$ could be about 7%. These are small enough that they would not have a significant effect on our analysis. We have also neglected the feeddown from D-wave states.

## III Parton Differential Cross Sections

In hadron collisions, bottomonium with transverse momentum $`p_T`$ of order $`m_b`$ or larger is produced, at leading order in $`\alpha _s`$, by parton fusion processes $`ijb\overline{b}+k`$. The differential cross section for producing a bottomonium state $`H`$ with momentum $`P`$ can be expressed in the schematic form

$$d\sigma [H(P)]_{\mathrm{fusion}}=f_{i/p}f_{j/\overline{p}}d\widehat{\sigma }[ijb\overline{b}(P,n)+k]O^H(n),$$
(9)
where there is an implied sum over the partons $`i,j,k`$ and over the $`b\overline{b}`$ channels $`n`$.

The order-$`\alpha _s^3`$ fusion cross section in (9) gives a good first approximation only if the transverse momentum is not too much larger or too much smaller than $`m_b`$. For $`p_Tm_b`$, the order-$`\alpha _s^3`$ fusion cross section for the channel $`b\overline{b}_8(^3S_1)`$ has the scaling behavior $`d\widehat{\sigma }/dp_T^21/p_T^4`$, while all other channels are suppressed by powers of $`m_b^2/p_T^2`$ at leading order. Parton processes with scaling behavior are called fragmentation processes. The fragmentation contributions to $`b\overline{b}`$ channels other than $`b\overline{b}_8(^3S_1)`$ enter at higher order in $`\alpha _s`$. The order-$`\alpha _s^3`$ fusion cross sections therefore underestimate the $`b\overline{b}`$ cross section in these channels at large $`p_T`$. However, the CDF data on bottomonium production extends only out to $`p_T=20`$ GeV , which is not large enough for fragmentation effects to dominate. In extracting the NRQCD matrix elements from that data, it should therefore be sufficient to use the fusion cross section (9).

The order-$`\alpha _s^3`$ fusion cross section in (9) also fails at small $`p_T`$. For some $`b\overline{b}`$ channels, including $`b\overline{b}_8(^1S_0)`$ and $`b\overline{b}_8(^3P_{0,2})`$, there is an order-$`\alpha _s^2`$ fusion cross section from the parton process $`ijb\overline{b}`$, which produces a $`b\overline{b}`$ pair with $`p_T=0`$. In these channels, the order-$`\alpha _s^3`$ fusion cross sections $`d\sigma /dp_T^2`$ diverge like $`1/p_T^2`$ as $`p_T0`$. The divergence in the integral of the cross section for $`ijb\overline{b}+k`$ is cancelled by the radiative corrections to the cross section for $`ijb\overline{b}`$, so that the cross section integrated over $`p_T`$ is finite order by order in $`\alpha _s`$. In order to obtain a smooth prediction for $`d\sigma /dp_T^2`$ in the small $`p_T`$ region, it is necessary to resum higher order corrections involving soft-gluon radiation. This resummation will have a significant effect on the shape of the $`p_T`$ distribution, and therefore on the values of the NRQCD matrix elements used to fit that distribution. We will avoid the complications due to soft-gluon radiation at small $`p_T`$ by using only the data from $`p_T>8`$ GeV to fit the NRQCD matrix elements.

We proceed to describe each of the factors in the fusion cross section (9) in more detail. We include the contributions from the following combinations of colliding partons: $`ij`$ = $`gg`$, $`gq`$, $`g\overline{q}`$, $`q\overline{q}`$, where $`q=u,d,s,c`$. The parton distributions $`f_{i/p}(x_1,\mu _F)`$ and $`f_{j/\overline{p}}(x_2,\mu _F)`$, which specify the momenta of the colliding partons, depend on a factorization scale $`\mu _F`$. We will consider the CTEQ5L and MRST98LO parton distribution functions. They are both obtained from leading order analyses, and thus can be used consistently with leading order parton cross sections. Explicit expressions for the parton differential cross sections $`d\widehat{\sigma }`$ are given in Ref.  and in Ref. . They are proportional to $`\alpha _s^3(\mu _R)`$, where $`\mu _R`$ is the renormalization scale, and they depend on the mass $`m_b`$ of the bottom quark. As part of the theoretical error, we will allow $`\mu _F`$ and $`\mu _R`$ to vary by factors of 2 from the central values $`\mu _T=(m_b^2+p_T^2)^{1/2}`$. This central value interpolates between half the partonic invariant mass at $`b\overline{b}`$ threshold and half the partonic invariant mass at large $`p_T`$ and central rapidity.

The cross sections also depend on two fundamental QCD parameters: $`\alpha _s`$ and $`m_b`$. We take the QCD coupling constant $`\alpha _s(\mu )`$ to run according to the one-loop formula, with the boundary value appropriate to the parton distribution function that is being used. For CTEQ5L, the coupling constant satisfies $`\alpha _s(M_Z)=0.127`$ and $`\alpha _s(m_b)=0.232`$. For MRST98LO, the coupling constant satisfies $`\alpha _s(M_Z)=0.125`$ and $`\alpha _s(m_b)=0.226`$.

The other fundamental QCD parameter in our calculation is the bottom quark mass $`m_b`$. There have been several recent determinations of $`m_b`$ using sum rules calculated to next-to-next-to-leading order accuracy with nonrelativistic resummation . A useful summary of these results is given by Beneke in Ref. . The value of the pole mass is rather unstable under radiative corrections compared to short-distance definitions of the mass, such as the running mass evaluated at its own scale, $`\overline{m}_b=m_b(\overline{m}_b)`$. Beneke’s best estimate for this mass is $`\overline{m}_b=4.23\pm 0.08`$ GeV. The definition of $`\overline{m}_b`$ is purely mathematical in character and not related to any physical thresholds involving the $`b`$ quark. Two definitions that are also relatively stable under radiative corrections and whose definitions are related to thresholds in the bottomonium system are the 1S mass, which is the perturbative energy of the lowest bound state, and the PS mass, which is the sum of the pole mass and some energy related to the potential between the $`b`$ and $`\overline{b}`$. Beneke’s best estimates for these masses are $`m_{b,1S}=4.77\pm 0.11`$ GeV and $`m_{b,PS}(2\mathrm{GeV})=4.57\pm 0.10`$ GeV. The relation between the two is given by a power series in $`\alpha _s`$:

$$m_{b,PS}(\mu )=m_{b,1S}\frac{4\alpha _s(\mu )}{3\pi }\mu +O(\alpha _s^2),$$
(10)
The difference between Beneke’s values for $`m_{b,1S}`$ and $`m_{b,PS}`$ is mostly accounted for by the order-$`\alpha _s`$ correction. We will choose the 1S mass as our prescription for the $`b`$ quark mass. Beneke’s central value for the PS mass differs by 2 standard deviations from the 1S mass. This difference should not be regarded as an ambiguity in the quark mass, because it could not be eliminated by a more precise determination of $`m_b`$. Instead its effects on the cross section could be decreased by calculating the next-to-leading order radiative correction to the parton cross sections. The uncertainty due to different prescriptions for the quark mass can therefore be regarded as part of the error due to radiative corrections.

## IV Color-singlet Matrix Elements

The color-singlet matrix elements for $`\mathrm{{\rm Y}}(nS)`$ can be determined phenomenologically from its decay rate into a lepton pair. The electronic decay rate of the $`\mathrm{{\rm Y}}(nS)`$, including the QCD radiative correction of order $`\alpha _s`$ and the first relativistic correction of order $`v^2`$, is

$$\mathrm{\Gamma }[\mathrm{{\rm Y}}(nS)e^+e^{}]=\frac{2\pi \alpha ^2e_b^2}{9m_b^2}(1\frac{8}{3}\frac{\alpha _s}{\pi }\frac{1}{3}\frac{M_{\mathrm{{\rm Y}}(nS)}2m_b}{2m_b})^2O_1^{\mathrm{{\rm Y}}(nS)}(^3S_1),$$
(11)
where $`e_b=1/3`$ is the bottom quark charge and $`m_b`$ is the 1S mass. The relativistic correction was first expressed in terms of $`M_{\mathrm{{\rm Y}}(nS)}2m_b`$ by Gremm and Kapustin . The vacuum saturation approximation, which is accurate up to corrections of order $`v^4`$, has been used to express the NRQCD matrix element that enters naturally in annihilation rates in terms of the corresponding production matrix element $`O_1^{\mathrm{{\rm Y}}(nS)}(^3S_1)`$. It has also been used to express the radiative and relativistic correction factor as a square. Setting $`m_b=4.77`$ GeV and $`\alpha _s(m_b)=0.22`$ and using the measured value for the decay rates, we obtain the values for the color-singlet matrix elements in Table II. In addition to the experimental errors in the decay rates, there are theoretical errors from relativistic corrections and from perturbative corrections. As a measure of the relativistic error of order $`v^4`$, we take the square of the largest of the order-$`v^2`$ corrections for the three $`\mathrm{{\rm Y}}(nS)`$ states: $`[(M_{\mathrm{{\rm Y}}(3S)}2m_b)/(3m_b)]^20.2`$%. As a measure of the perturbative error from higher orders in $`\alpha _s`$, we take the square of the order-$`\alpha _s`$ correction in (11): $`[16\alpha _s/(3\pi )]^214`$%. The error bars quoted in Table II are obtained by combining the experimental, relativistic, and perturbative errors in quadrature. The error bars are dominated by the 14% perturbative error, except in the case of the $`\mathrm{{\rm Y}}(3S)`$ for which the experimental error is 16%. The values for $`O_1^{\mathrm{{\rm Y}}(nS)}(^3S_1)`$ in Table II are larger by a factor of 3 than those given in Table I of the first paper in Ref.  because of a normalization error in the Table. This normalization error did not appear in the cross sections and therefore did not affect the results.

There is no data that can be used for phenomenological determinations of the color-singlet matrix elements for the P-wave states. However the color-singlet matrix elements for both the S-wave and P-wave states can be estimated using wavefunctions from potential models. Using the vaccuum-saturation approximation, the color-singlet matrix element for $`\mathrm{{\rm Y}}(nS)`$ can be expressed in terms of its radial wavefunction at the origin, while that for $`\chi _{bJ}(nP)`$ can be expressed in terms of the derivative of its radial wavefunction at the origin:

$`O_1^{\mathrm{{\rm Y}}(nS)}(^3S_1)`$ $``$ $`{\displaystyle \frac{9}{2\pi }}|R_{\mathrm{{\rm Y}}(nS)}(0)|^2,`$ (12)
$`O_1^{\chi _{bJ}(nP)}(^3P_J)`$ $``$ $`(2J+1){\displaystyle \frac{9}{2\pi }}|R_{\chi _b(nP)}^{}(0)|^2.`$ (13)
Eichten and Quigg have tabulated the radial wavefunctions and their derivatives at the origin for 4 potential models that reproduce the observed bottomonium spectrum . As estimates of the color-singlet matrix elements (12) and (13), we take the mean values from the 4 potential models in Ref. . The mean values are tabulated in Table II. The errors are the root-mean-square deviations of the 4 potential-model values. The potential model estimates of $`O_1^{\mathrm{{\rm Y}}(nS)}(^3S_1)`$ are consistent with the phenomenological values, but have larger error bars. This gives us some confidence in the potential-model estimates for $`O_1^{\chi _{b0}(nP)}(^3P_0)`$ in Table II. These values are consistent within errors with those used in the analysis of Ref. .

The most accurate determination of the color-singlet matrix elements for the lowest bottomonium states may eventually come from lattice gauge theory. The corresponding annihilation matrix elements can be readily calculated using lattice simulations of NRQCD . The NRQCD collaboration has calculated the wavefunctions for the lowest S-wave states  and P-wave states  of bottomonium. Inserting these wavefunctions into the expressions in (12) and (13), we obtain the estimates of the color-singlet matrix elements in Table II. The largest errors in the lattice calculations come from matching of lattice NRQCD operators with continuum NRQCD operators and from the omission of dynamical light quarks. Both errors could be as large as 25% in the present calculations. It is premature to quote error bars for the lattice gauge theory results in Table II.

We will adopt the phenomenological values of $`O_1^{\mathrm{{\rm Y}}(nS)}(^3S_1)`$ in Table II and the potential-model values for $`O_1^{\chi _{b0}(nP)}(^3S_1)`$. Using the branching fractions in Table I, we can form the linear combinations that appear in the expressions (7) for the inclusive $`\mathrm{{\rm Y}}(nS)`$ cross sections. These are tabulated in Table III. As indicated by the zeros in the first row, we neglect feeddown from the $`\chi _b(3P)`$ states.

## V Color-octet Matrix Elements

The color-octet NRQCD matrix elements are phenomenological parameters that can only be determined from experimental data. We first extract the inclusive color-octet matrix elements for $`\mathrm{{\rm Y}}(nS)`$ from the CDF measurements of the inclusive $`\mathrm{{\rm Y}}(nS)`$ cross sections. We then extract direct color-octet matrix elements for $`\chi _{bJ}(nP)`$ from the CDF measurements of the fraction of $`\mathrm{{\rm Y}}(1S)`$ coming from $`\chi _b`$’s. This gives us enough information to determine the direct color-octet matrix elements for $`\mathrm{{\rm Y}}(nS)`$ .

### A Inclusive Matrix Elements for S-waves

The inclusive $`\mathrm{{\rm Y}}(nS)`$ cross sections depend linearly on the inclusive matrix elements defined in (8). The inclusive color-singlet matrix elements are given in Table III. We can extract the inclusive color-octet matrix elements from the CDF measurements of the $`\mathrm{{\rm Y}}(nS)`$ cross sections . The differential cross sections integrated over rapidities $`|y|<0.4`$ have been measured out to $`p_T=`$ 20 GeV . The CDF data on $`Bd\sigma /dp_T`$ for $`\mathrm{{\rm Y}}(nS)`$, where $`B`$ is the branching fraction for $`\mathrm{{\rm Y}}(nS)\mu ^+\mu ^{}`$, are shown in Figures 1, 2, and 3 for $`\mathrm{{\rm Y}}(1S)`$, $`\mathrm{{\rm Y}}(2S)`$, and $`\mathrm{{\rm Y}}(3S)`$, respectively. We avoid the problem of carrying out soft-gluon resummation to determine the shapes of the theoretical $`p_T`$ distributions at low $`p_T`$ by using only the data from $`p_T>`$ 8 GeV to fit the color-octet matrix elements. This leaves 5 $`p_T`$ bins for $`\mathrm{{\rm Y}}(1S)`$ and 3 $`p_T`$ bins each for $`\mathrm{{\rm Y}}(2S)`$ and $`\mathrm{{\rm Y}}(3S)`$.

The inclusive $`\mathrm{{\rm Y}}(nS)`$ cross sections depend on the inclusive color-octet matrix elements through the linear combination $`[O_8(^1S_0))+mO_8(^3P_0)/m_b^2]+nO_8(^3S_1))`$, where $`m`$ varies from 4.6 at $`p_T=8`$ GeV to 3.4 at $`p_T=20`$ GeV, while $`n`$ varies from 1.0 at $`p_T=8`$ GeV to 6.3 at $`p_T=20`$ GeV. The parameters $`O_8(^1S_0)`$ and $`O_8(^3P_0)`$ can not be determined independently, because the corresponding parton cross sections have similar dependences on $`p_T`$. We therefore carry out our analysis under the two extreme assumptions that either the $`O_8(^1S_0)`$ term or the $`O_8(^3P_0)`$ term dominates and that the other can be neglected. Assuming both matrix elements are positive, the truth will be somewhere in between the two extremes. We will take the difference between the two extremes as part of the theoretical error.

For a given choice of the parton distributions and the scales $`\mu _F`$ and $`\mu _R`$, we can integrate the $`\mathrm{{\rm Y}}(nS)`$ differential cross section (7) over $`|y|<0.4`$ and over each $`p_T`$ bin. We determine the best fits for the color-octet matrix elements by minimizing the $`\chi ^2`$ associated with the sum over $`p_T`$ bins. Fixing $`\mu _F=\mu _R=\mu _T`$ and assuming that $`O_8(^3P_0)`$ is negligible, we obtain the values of $`O_8(^3S_1)`$ and $`O_8(^1S_0)`$ in the 1st and 3rd columns of Table IV. Repeating the analysis but assuming that $`O_8(^1S_0)`$ is negligible, we obtain the values of $`O_8(^3S_1)`$ and $`O_8(^3P_0)`$ in the 2nd and 4th columns. The first errors in Table IV are extracted from the matrix of second derivatives of the $`\chi ^2`$ function. There is also an error from varying the renormalization and factorization scales $`\mu _F`$ and $`\mu _R`$. These errors are large, but since they are highly correlated, we have separated them out as a second error in Table IV. The upper and lower errors are the shifts in the matrix elements that minimize $`\chi ^2`$ when $`\mu _R`$ and $`\mu _F`$ are changed from the central value $`\mu _T=\sqrt{m_b^2+p_T^2}`$ by multiplicative factors of 2 and 1/2, respectively. The error from varying $`m_b`$ is also highly correlated, but it is smaller and we do not list it separately. It can be taken into account when we use the matrix elements to calculate other observables.

Our fits for $`Bd\sigma /dp_Tdy`$ at $`y=0`$ for inclusive $`\mathrm{{\rm Y}}(1S)`$, $`\mathrm{{\rm Y}}(2S)`$, and $`\mathrm{{\rm Y}}(3S)`$ are compared to the CDF data in Figures 1, 2, and 3. The error bands reflect the statistical uncertainties in the matrix elements. The fits are reasonably good in the region $`p_T>8`$ GeV that we used for fitting. At low $`p_T`$, our fits for $`d\sigma /dp_T`$ behave like $`1/p_T`$, because we have not implemented the effects of soft-gluon radiation on the shape of the $`p_T`$-distribution. The fits therefore diverge from the data below $`p_T=8`$ GeV. For $`\mathrm{{\rm Y}}(2S)`$, the central curve becomes negative at small $`p_T`$ because our fit gives a negative central value for $`O_8(^1S_0)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(2S)}`$ or $`O_8(^3P_0)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(2S)}`$. If we had fit the color-octet matrix elements using the data for $`p_T>4`$ GeV, instead of only the data for $`p_T>8`$ GeV, the central values for $`O_8(^1S_0)_{\mathrm{inc}}`$ or $`O_8(^3P_0)_{\mathrm{inc}}`$ would also have been negative for $`\mathrm{{\rm Y}}(1S)`$ and $`\mathrm{{\rm Y}}(3S)`$.

In Figures 1, 2, and 3, the color-singlet model predictions from order-$`\alpha _s^3`$ fusion processes are shown as dashed lines. At the largest values of $`p_T`$ shown, the predictions fall more than an order of magnitude below the data. The color-singlet model prediction for $`\mathrm{{\rm Y}}(3S)`$ indicates the shape of the $`b\overline{b}_1(^3S_1)`$ cross section. The color-singlet model predictions for $`\mathrm{{\rm Y}}(1S)`$ and $`\mathrm{{\rm Y}}(2S)`$ behave very differently at small $`p_T`$, because they receive contributions from decays of $`\chi _{bJ}(nP)`$. The predictions diverge as $`p_T0`$, because the cross sections $`d\widehat{\sigma }/dp_T`$ for $`b\overline{b}_1(^3P_{0,2})`$ behave like $`1/p_T`$. In order to obtain the correct threshold behavior in these channels, it would be necessary to resum the effects of soft-gluon radiation.

### B  Direct Matrix Elements for P-waves

The color-octet matrix elements $`O_8(^3S_1)`$ for the $`\chi _b`$’s can be determined from CDF measurements of the fractions of $`\mathrm{{\rm Y}}(1S)`$’s that come from $`\chi _b(1P)`$’s and from $`\chi _b(2P)`$’s . The important feeddown decays for $`\chi _{bJ}(1P)`$ and $`\chi _{bJ}(2P)`$ proceed through the $`\mathrm{{\rm Y}}(2S)`$ and $`\mathrm{{\rm Y}}(3S)`$, respectively. The fractions $`F_{\chi _b(nP)}^{\mathrm{{\rm Y}}(1S)}`$ of $`\mathrm{{\rm Y}}(1S)`$’s from $`\chi _b(nP)`$ therefore satisfy

$`F_{\chi _b(1P)}^{\mathrm{{\rm Y}}(1S)}\sigma [\mathrm{{\rm Y}}(1S)]_{\mathrm{inc}}`$ $`=`$ $`{\displaystyle \underset{J}{}}B_{\chi _{bJ}(1P)\mathrm{{\rm Y}}(1S)}\sigma [\chi _{bJ}(1P)]`$ (15)
$`+\left({\displaystyle \underset{J}{}}B_{\mathrm{{\rm Y}}(2S)\chi _{bJ}(1P)}B_{\chi _{bJ}(1P)\mathrm{{\rm Y}}(1S)}\right)\sigma [\mathrm{{\rm Y}}(2S)]_{\mathrm{inc}},`$
$`F_{\chi _b(2P)}^{\mathrm{{\rm Y}}(1S)}\sigma [\mathrm{{\rm Y}}(1S)]_{\mathrm{inc}}`$ $`=`$ $`{\displaystyle \underset{J}{}}B_{\chi _{bJ}(2P)\mathrm{{\rm Y}}(1S)}\sigma [\chi _{bJ}(2P)]`$ (17)
$`+\left({\displaystyle \underset{J}{}}B_{\mathrm{{\rm Y}}(3S)\chi _{bJ}(2P)}B_{\chi _{bJ}(2P)\mathrm{{\rm Y}}(1S)}\right)\sigma [\mathrm{{\rm Y}}(3S)]_{\mathrm{inc}}.`$
Using the branching fractions from Table I, the coefficients of $`\sigma [\mathrm{{\rm Y}}(2S)]_{\mathrm{inc}}`$ and $`\sigma [\mathrm{{\rm Y}}(3S)]_{\mathrm{inc}}`$ in (15) and (17) are $`(3.8\pm 0.7)`$% and $`(1.8\pm 0.2)`$%, respectively.

The CDF result for the fractions of $`\mathrm{{\rm Y}}(1S)`$’s with $`|y|<0.4`$ and $`p_T>8`$ GeV that come from $`\chi _b(1P)`$’s and $`\chi _b(2P)`$’s are

$`F_{\chi _b(1P)}^{\mathrm{{\rm Y}}(1S)}`$ $`=`$ $`(27.1\pm 8.1)\%,`$ (18)
$`F_{\chi _b(2P)}^{\mathrm{{\rm Y}}(1S)}`$ $`=`$ $`(10.5\pm 4.8)\%,`$ (19)
where we have added the statistical and systematic errors in quadrature. The inclusive $`\mathrm{{\rm Y}}(nS)`$ cross sections in (15) and (17) are the cross sections integrated over $`|y|<0.4`$ and $`p_T>8`$ GeV. Using the CDF measurements of these cross sections, (18) and (19) reduce to the following constraints on the cross sections for $`\chi _{bJ}(1P)`$ and $`\chi _{bJ}(2P)`$:

$`{\displaystyle \underset{J}{}}B_{\chi _{bJ}(nP)\mathrm{{\rm Y}}(1S)}\sigma [\chi _{bJ}(nP)]`$ $`=`$ $`(0.85\pm 0.29)\mathrm{nb},n=1,`$ (20)
$`=`$ $`(0.34\pm 0.17)\mathrm{nb},n=2.`$ (21)
The branching fractions and the associated errors are given in Table I.

The theoretical cross sections for $`\sigma [\chi _{bJ}(nP)]`$, $`J=0,1,2`$, are obtained by integrating (4) over the appropriate region of $`y`$ and $`p_T`$. The constraints (20) and (21) are then linear equations for $`O_8^{\chi _{b0}(nP)}(^3S_1)`$, which give the values in Table V. The first error is obtained by setting $`\mu _F=\mu _R=\mu _T`$ and combining in quadrature the experimental error from (20) or (21), the error from the branching fractions in Table I, and the errors from the color-singlet matrix elements in Table II. The second upper and lower errors in Table V are the shifts in the matrix elements when $`\mu _R`$ and $`\mu _F`$ are changed from their central values by multiplicative factors of 2 and 1/2, respectively.

### C Direct Matrix Elements for S-waves

The NRQCD matrix elements in Table IV can be used to calculate the inclusive $`\mathrm{{\rm Y}}(nS)`$ cross sections. To calculate the direct $`\mathrm{{\rm Y}}(nS)`$ cross sections, we must extract direct color-octet matrix elements for the $`\mathrm{{\rm Y}}(nS)`$ states from the inclusive color-octet matrix elements given in Table IV. The linear combinations of matrix elements determined by the inclusive $`\mathrm{{\rm Y}}(1S)`$ cross sections are

$`O_8(^3S_1)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(1S)}`$ $`=`$ $`O_8^{\mathrm{{\rm Y}}(1S)}(^3S_1)+0.311O_8^{\mathrm{{\rm Y}}(2S)}(^3S_1)+0.112O_8^{\mathrm{{\rm Y}}(3S)}(^3S_1)`$ (23)
$`+2.15O_8^{\chi _{b0}(1P)}(^3S_1)+1.08O_8^{\chi _{b0}(2P)}(^3S_1),`$
$`O_8(^1S_0)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(1S)}`$ $`=`$ $`O_8^{\mathrm{{\rm Y}}(1S)}(^1S_0)+0.311O_8^{\mathrm{{\rm Y}}(2S)}(^1S_0)+0.112O_8^{\mathrm{{\rm Y}}(3S)}(^1S_0),`$ (24)
$`O_8(^3P_0)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(1S)}`$ $`=`$ $`O_8^{\mathrm{{\rm Y}}(1S)}(^3P_0)+0.311O_8^{\mathrm{{\rm Y}}(2S)}(^3P_0)+0.112O_8^{\mathrm{{\rm Y}}(3S)}(^3P_0).`$ (25)
The errors on the branching fractions in (23)–(25) have been suppressed, but they are $`(0.311\pm 0.016)`$, $`(0.112\pm 0.005)`$, $`(2.15\pm 0.31)`$, and $`(1.08\pm 0.08)`$. The linear combinations determined by the inclusive $`\mathrm{{\rm Y}}(2S)`$ cross sections are

$`O_8(^3S_1_{\mathrm{inc}}^{\mathrm{{\rm Y}}(2S)}`$ $`=`$ $`O_8^{\mathrm{{\rm Y}}(2S)}(^3S_1)+0.106O_8^{\mathrm{{\rm Y}}(3S)}(^3S_1)+1.49O_8^{\chi _{b0}(2P)}(^3S_1),`$ (26)
$`O_8(^1S_0)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(2S)}`$ $`=`$ $`O_8^{\mathrm{{\rm Y}}(2S)}(^1S_0)+0.106O_8^{\mathrm{{\rm Y}}(3S)}(^1S_0),`$ (27)
$`O_8(^3P_0)_{\mathrm{inc}}^{\mathrm{{\rm Y}}(2S)}`$ $`=`$ $`O_8^{\mathrm{{\rm Y}}(2S)}(^3P_0)+0.106O_8^{\mathrm{{\rm Y}}(3S)}(^3P_0).`$ (28)
The errors on the branching fractions in (26)–(28) are $`(0.106\pm 0.008)`$ and $`(1.49\pm 0.17)`$. Using the color-octet matrix elements for inclusive $`\mathrm{{\rm Y}}(nS)`$ in Table IV and the color-octet matrix elements for direct $`\chi _b(nP)`$ in Table V, we obtain the color-octet matrix elements for direct $`\mathrm{{\rm Y}}(nS)`$ in Table V.

Our analysis gives a negative value consistent with zero for the matrix elements $`O_8^{\mathrm{{\rm Y}}(2S)}(^1S_0)`$ or $`O_8^{\mathrm{{\rm Y}}(2S)}(^3P_0)`$. Our values for $`O_8^{\chi _{b0}(2P)}(^3S_1)`$, $`O_8^{\chi _{b0}(1P)}(^3S_1)`$, and $`O_8^{\mathrm{{\rm Y}}(1S)}(^3S_1)`$ are also consistent with zero given the statistical error. The only direct color-octet matrix elements that differ from zero by two or more statistical error bars are $`O_8^{\mathrm{{\rm Y}}(2S)}(^3S_1)`$ and $`O_8^{\mathrm{{\rm Y}}(1S)}(^1S_0)`$ or $`O_8^{\mathrm{{\rm Y}}(1S)}(^3P_0)`$.

We now compare our values for the matrix elements with those obtained by Cho and Leibovich in their pioneering analysis of bottomonium production at the Tevatron. Their analysis was based on the CDF data sample from Run IA of the Tevatron . To reduce the errors associated with the shape of the $`p_T`$ distribution at small $`p_T`$, they used only the data from $`p_T>3.5`$ GeV in their analysis. The data was insufficient to determine all the matrix elements, so they estimated the matrix elements $`O_8^{\mathrm{{\rm Y}}(nS)}(^3S_1)`$ by applying scaling relations to the corresponding matrix elements in the charmonium sector. Their value for $`O_8^{\mathrm{{\rm Y}}(1S)}(^3S_1)`$ is consistent with ours to within our large error bars, but their value for $`O_8^{\mathrm{{\rm Y}}(2S)}(^3S_1)`$ is smaller than ours by about a factor of 40. They used the CDF data to fit the matrix elements $`O_8^{\chi _b(nP)}(^3S_1)`$ and the linear combinations $`M_5^{\mathrm{{\rm Y}}(nS)}=O_8^{\mathrm{{\rm Y}}(nS)}(^1S_0)+5O_8^{\mathrm{{\rm Y}}(nS)}(^3P_0)/m_b^2`$. Their values for $`O_8^{\chi _b(nP)}(^3S_1)`$ are comparable to ours in magnitude, but they have much smaller error bars. Their values for $`M_5^{\mathrm{{\rm Y}}(nS)}`$ differ from zero by only about one error bar. In our analysis, we included the matrix elements $`O_8^{\mathrm{{\rm Y}}(nS)}(^3S_1)`$ in the list of those to be fit to the CDF data. The much higher quality of the CDF data from Run IB of the Tevatron  allowed us to carry out a reasonable fit using the data restricted to $`p_T>8`$ GeV.

## VI Inclusive Cross Sections for Spin-triplet States

Having determined the most important matrix elements for the production of the spin-triplet bottomonium states, we can use them to calculate the cross sections for these states in other high energy processes. In particular, we can calculate their cross sections in Run II of the Tevatron in which the center-of-mass energy will be increased from 1.8 TeV to 2.0 TeV. In order to cancel the large theoretical errors, such as those from the uncertainties in the matrix elements and from the choice of scale, we normalize the cross sections to that for inclusive $`\mathrm{{\rm Y}}(1S)`$ at $`\sqrt{s}=1.8`$ TeV. For each bottomonium state $`H`$, we define the ratio

$`R^H(\sqrt{s})`$ $`=`$ $`{\displaystyle \frac{\sigma [H;\sqrt{s}]}{\sigma [\mathrm{inclusive}\mathrm{{\rm Y}}(1S);\sqrt{s}=1.8\mathrm{TeV}]}},`$ (29)
where the cross sections are integrated over $`p_T>8`$ GeV and over $`|y|<0.4`$.

To calculate the inclusive cross section for $`\mathrm{{\rm Y}}(nS)`$, we simply use the inclusive color-singlet matrix elements from Table III and the inclusive color-octet matrix elements from Table IV. To calculate the inclusive cross section for $`\chi _{bJ}(nP)`$, we must first compute the direct cross sections for $`\chi _{bJ}(nP)`$ and the higher bottomonium states using the direct color-singlet matrix elements from Table II and the direct color-octet matrix elements from Table V, and then combine them using the inclusive branching fractions in Table I. The resulting ratios $`R^H`$ shown in Table VI are the averages of the 4 values obtained by using either the CTEQ5L or MRST98LO parton distributions and either setting $`O_8(^1S_0)=0`$ or $`O_8(^3P_0)=0`$. The error bars come from combining in quadrature the statistical errors in the matrix elements from the Tables, the error from varying $`\mu `$ by a factor of two from its central value, the difference between using the CTEQ5L and MRST98LO parton distributions, the difference between setting $`O_8(^1S_0)=0`$ and $`O_8(^3P_0)=0`$, and the error from varying $`m_b`$. The error bars in the numerator and denominator of (29) are both large, but they are highly correlated and tend to cancel in the ratio. The largest contributions to the error bars are the statistical errors on the matrix elements, with the exception of $`\mathrm{{\rm Y}}(1S)`$, for which the largest contribution comes from varying the scale. For $`\mathrm{{\rm Y}}(2S)`$ and $`\mathrm{{\rm Y}}(3S)`$, the results in Table VI for $`\sqrt{s}=1.8`$ TeV are consistent with the actual CDF measurements, which give $`0.61\pm 0.12`$ for $`\mathrm{{\rm Y}}(2S)`$ and $`0.29\pm 0.12`$ for $`\mathrm{{\rm Y}}(3S)`$. When the center-of-mass energy is increased from 1.8 TeV to 2.0 TeV, all the cross sections increase by about 16%. The increase depends on $`p_T`$, changing from about 15% at $`p_T=8`$ TeV to about 19% at $`p_T=20`$ TeV.

## VII Direct Cross Sections for Spin-singlet States

Having determined the most important matrix elements for the production of the spin-triplet bottomonium states, we can also use them to calculate the production rate of the spin-singlet states $`\eta _b(nS)`$ and $`h_b(nP)`$. The matrix elements for these states are related to those of the corresponding spin-triplet states by the approximate spin symmetry of NRQCD. Spin symmetry relates the matrix elements for $`\eta _b(nS)`$ to those for $`\mathrm{{\rm Y}}(nS)`$:

$`O_1^{\eta _b(nS)}(^1S_0)`$ $`=`$ $`{\displaystyle \frac{1}{3}}O_1^{\mathrm{{\rm Y}}(nS)}(^3S_1),`$ (30)
$`O_8^{\eta _b(nS)}(^1P_1)`$ $`=`$ $`{\displaystyle \frac{1}{3}}O_8^{\mathrm{{\rm Y}}(nS)}(^3P_0),`$ (31)
$`O_8^{\eta _b(nS)}(^1S_0)`$ $`=`$ $`{\displaystyle \frac{1}{3}}O_8^{\mathrm{{\rm Y}}(nS)}(^3S_1),`$ (32)
$`O_8^{\eta _b(nS)}(^3S_1)`$ $`=`$ $`{\displaystyle \frac{1}{3}}O_8^{\mathrm{{\rm Y}}(nS)}(^1S_0).`$ (33)
The direct matrix elements for $`\eta _b(nS)`$ can therefore be read off from those for $`\mathrm{{\rm Y}}(nS)`$ in Tables II and V. The spin-symmetry relations have been used to calculate the cross sections for producing the $`\eta _c`$ at the Tevatron , in photoproduction and electroproduction, , and at Hera-B . In Refs.  and , the factor of $`\frac{1}{3}`$ which comes from the number of spin states was omitted from the spin-symmetry relations (30)-(33). Thus their estimates of the cross sections may be too large by a factor of 3. Spin symmetry relates the matrix elements for $`h_b(nP)`$ to those for $`\chi _{b0}(nP)`$:

$`O_1^{h_b(nP)}(^1P_1)`$ $`=`$ $`3O_1^{\chi _{b0}(nP)}(^3P_0),`$ (34)
$`O_8^{h_b(nP}(^1S_0)`$ $`=`$ $`3O_8^{\chi _{b0}(nP)}(^3S_1).`$ (35)
The direct matrix elements for $`h_b(nP)`$ can therefore be read off from those for $`\chi _b(nP)`$ in Tables II and V. These relations were first used by Fleming and Mehen to calculate the cross section for the photoproduction of $`h_c`$ . They have also been used to calculate the cross sections for $`h_c`$ at Hera-B .

To calculate the cross sections for the direct production of the spin-singlet states at the Tevatron, we need the appropriate parton differential cross sections $`d\widehat{\sigma }`$. Explicit expressions for most of those that are needed are given in Ref.  and in Ref. . The exception is the cross sections for producing $`b\overline{b}_8(^1P_1)`$, which contributes to $`\eta _b`$ production. The cross sections for the production of $`b\overline{b}_8(^1P_1)`$ in $`q\overline{q}`$, $`gq`$, and $`gg`$ scattering are

$`{\displaystyle \frac{d\sigma }{d\widehat{t}}}\left(q\overline{q}^1P_1^{(8)}g\right)`$ $`=`$ $`{\displaystyle \frac{O_8(^1P_1)}{16\pi \widehat{s}^2}}{\displaystyle \frac{4(4\pi \alpha _s)^3}{9M^3}}{\displaystyle \frac{\widehat{t}^2+\widehat{u}^2}{\widehat{s}(\widehat{s}M^2)^2}},`$ (36)
$`{\displaystyle \frac{d\sigma }{d\widehat{t}}}\left(gq^1P_1^{(8)}q\right)`$ $`=`$ $`{\displaystyle \frac{O_8(^1P_1)}{16\pi \widehat{s}^2}}{\displaystyle \frac{(4\pi \alpha _s)^3}{6M^3}}{\displaystyle \frac{\widehat{s}^2+\widehat{u}^2}{(\widehat{t})(M^2\widehat{t})^2}},`$ (37)
$`{\displaystyle \frac{d\sigma }{d\widehat{t}}}\left(gg^1P_1^{(8)}g\right)`$ $`=`$ $`{\displaystyle \frac{O_8(^1P_1)}{16\pi \widehat{s}^2}}{\displaystyle \frac{(4\pi \alpha _s)^3}{36M^3}}{\displaystyle \frac{1}{z^2\widehat{s}(\widehat{s}M^2)^3(\widehat{s}M^2+z^2)^3}}`$ (46)
$`\{27\widehat{s}z^2(\widehat{s}^84\widehat{s}^6z^2+\widehat{s}^4z^4+\widehat{s}^2z^6+z^8)`$
$`+M^2(27\widehat{s}^{10}243\widehat{s}^8z^2+697\widehat{s}^6z^4665\widehat{s}^4z^6+346\widehat{s}^2z^827z^{10})`$
$`M^4\widehat{s}(135\widehat{s}^8702\widehat{s}^6z^2+1340\widehat{s}^4z^41087\widehat{s}^2z^6+135z^8)`$
$`+M^6(324\widehat{s}^81134\widehat{s}^6z^2+1557\widehat{s}^4z^4698\widehat{s}^2z^6+54z^8)`$
$`M^8\widehat{s}(486\widehat{s}^61091\widehat{s}^4z^2+882\widehat{s}^2z^492z^6)`$
$`+M^{10}(486\widehat{s}^6616\widehat{s}^4z^2+374\widehat{s}^2z^427z^6)`$
$`M^{12}\widehat{s}(324\widehat{s}^4211\widehat{s}^2z^2+38z^4)`$
$`+M^{14}\widehat{s}^2(135\widehat{s}^238z^2)27M^{16}\widehat{s}^3\}.`$
where $`z^2=\widehat{t}\widehat{u}`$.

We proceed to calculate the cross sections for the direct production of the spin-singlet states at the Tevatron at center-of-mass energies $`\sqrt{s}=1.8`$ TeV and 2.0 TeV. To minimize the effect of the highly correlated errors, we calculate the ratio (29) of the direct cross section integrated over $`p_T>8`$ GeV and over $`|y|<0.4`$ to the corresponding inclusive cross section for $`\mathrm{{\rm Y}}(1S)`$ at 1.8 TeV. The resulting predictions are shown in Table VI. The cross sections for the $`h_b(nP)`$ states are small compared to those for $`\mathrm{{\rm Y}}(1S)`$ and they have large error bars. The cross sections for $`\eta _b(nS)`$ are predicted to be several times larger than those for $`\mathrm{{\rm Y}}(nS)`$ and they have reasonably small error bars. These predictions should be fairly reliable, because the cross sections are for the same kinematical region as the data used to extract the matrix elements. When the center-of-mass energy is increased from 1.8 TeV to 2 TeV, all the cross sections increase by about 16%.

We can make a rough estimate of the cross sections integrated over all $`p_T`$ by assuming that the spin-singlet cross sections have the same shape at small $`p_T`$ as the $`\mathrm{{\rm Y}}(1S)`$ cross section. The measured inclusive $`\mathrm{{\rm Y}}(1S)`$ cross section at central rapidity integrated over all $`p_T`$ up to 20 GeV satisfies $`Bd\sigma /dy=690\pm 25`$ pb, where $`B2.5\%`$ is the branching fraction of $`\mathrm{{\rm Y}}(1S)`$ into $`\mu ^+\mu ^{}`$. The cross section integrated only over $`p_T>8`$ GeV satisfies $`Bd\sigma /dy=106\pm 7`$ pb. The ratio of these cross sections is $`6.5\pm 0.5`$. Multiplying the inclusive $`\mathrm{{\rm Y}}(1S)`$ cross section $`d\sigma /dy=28`$ nb by the factor of 6.5 to take into account the small $`p_T`$ region and by the ratio 4.3 from Table VI, we find that the cross section for $`\eta _b(1S)`$ integrated over all $`p_T`$ should be approximately $`d\sigma /dy=800`$ nb.

The cross section for $`\eta _b(1S)`$ indicates that this state must have been produced in abundance in Run I of the Tevatron. However the $`\eta _b(1S)`$ can be observed only if it has a large enough branching fraction into a decay mode that can be triggered upon. One possibility is the double-$`J/\psi `$ decay $`\eta _b(1S)J/\psi +J/\psi `$, followed by the decays $`J/\psi \mu ^+\mu ^{}`$. The decay $`\eta _b(1S)J/\psi +J/\psi `$ has essentially the same kinematics as the decay $`\eta _c\varphi \varphi `$, except that all masses are scaled up by a factor of 3. Thus the branching fraction for $`\eta _b(1S)J/\psi +J/\psi `$ could be as large as that for $`\eta _c\varphi \varphi `$, which is approximately $`7\times 10^3`$. We can obtain a lower bound on the branching fraction by using the fact that in the limit $`m_b\mathrm{}`$ with $`m_c`$ fixed, the branching fraction for $`\eta _b(1S)J/\psi +J/\psi `$ scales like $`1/m_b^4`$ . If $`m_c`$ and $`m_b`$ were both in this scaling region, then the branching fraction into light $`J^{PC}=1^{}`$ mesons would be smaller for $`\eta _b(1S)`$ than for $`\eta _c`$ by a factor of $`(m_c/m_b)^4`$, which is about $`10^2`$. Since we are not deep into this scaling region, the suppression should be smaller than this. Thus the branching fraction for $`\eta _b(1S)J/\psi +J/\psi `$ should be in the range between $`7\times 10^5`$ and $`7\times 10^3`$. Multiplying by the branching fractions of 6% for each of the decays $`J/\psi \mu ^+\mu ^{}`$, our estimate for the branching fraction for $`\eta _b(1S)J/\psi +J/\psi \mu ^+\mu ^{}+\mu ^+\mu ^{}`$ is $`B2.5\times 10^{6\pm 1}`$. The cross section for producing this particular decay mode of $`\eta _b(1S)`$ is therefore $`Bd\sigma /dy2`$ pb, give or take a factor of 10. Multiplying by the rapidity interval 0.8 and by the integrated luminosity of about 100 pb<sup>-1</sup> in Run I of the Tevatron, we obtain between 16 and 1600 produced events.

We must also take into account the acceptances and efficiencies for observing the decays $`J/\psi \mu ^+\mu ^{}`$. These can be estimated using the CDF data on the production of prompt $`J/\psi `$ in Run IA of the Tevatron . Based on the observation of about 22,000 $`J/\psi \mu ^+\mu ^{}`$ candidates with $`p_T>5`$ GeV and pseudorapidity $`|\eta |<0.6`$ in a 15 pb<sup>-1</sup> data sample, they measured the cross section in that region of $`p_T`$ and $`\eta `$ to be $`B\sigma 17`$ nb. We infer that the product of the acceptance and the efficiency is roughly $`ϵ0.09`$. Multiplying the number of events that are produced by $`ϵ^2`$, we get between 0.13 and 13 observed events. Thus this cross section may be large enough to be observed in Run I of the Tevatron. In Run II, the integrated luminosity will be larger by a factor of 20 and there will be significant improvements in the acceptances and efficiencies for observing muons in both the CDF and D0 detectors. Thus the $`\eta _b(1S)`$ should certainly be observed in Run II through the decay $`\eta _bJ/\psi +J/\psi `$.

## VIII Discussion

We have carried out an updated NRQCD analysis of the CDF data on the production of spin-triplet bottomonium states from Run I of the Tevatron. In spite of using only the data from $`p_T>8`$ GeV, we were able to extract all the relevant color-octet matrix elements directly from the data. Only one of the 8 color-octet matrix elements comes out with a negative central value, but several others are also consistent with zero to within errors. In our analysis, we distinguished between the inclusive color-octet matrix elements that can be used to compute inclusive $`\mathrm{{\rm Y}}(nS)`$ cross sections and the direct color-octet matrix elements required to compute direct $`\mathrm{{\rm Y}}(nS)`$ cross sections and, by spin symmetry, direct $`\eta _b(nS)`$ cross sections.

The most serious deficiency in our analysis was our failure to take into account the effects of soft-gluon radiation that are needed to give a smooth $`p_T`$ distribution at small $`p_T`$. This forced us to use only the small fraction of the data from $`p_T>8`$ GeV to fit the color-octet matrix elements, which led to large errors in these matrix elements. If these matrix elements are used to predict bottomonium cross sections in other high energy processes, the predictions will probably have large error bars. Our theoretical cross sections also diverge from the CDF data below $`p_T=8`$ GeV, which gives us another reason to be cautious in applying our matrix elements to other high energy processes.

An analysis that deals properly with the small $`p_T`$ region could take full advantage of the CDF data and therefore determine the color-octet matrix elements more accurately. Such an analysis requires a prescription for combining the leading-order cross sections for $`ijb\overline{b}+k`$ with the next-to-leading order cross sections for $`ijb\overline{b}`$ recently calculated by Petrelli et al.  to get a smooth $`p_T`$ distribution near $`p_T=0`$. The matrix elements extracted from such an analysis should give reliable predictions for observables at the Tevatron that are dominated by low $`p_T`$. They should also give reliable predictions for bottomonium production in other high energy processes.

Our analysis should give reliable predictions for the cross sections of the spin-singlet states $`\eta _b(nS)`$ and $`h_b(nP)`$ at the Tevatron. We find that the direct cross section for $`\eta _b(1S)`$ at $`p_T>8`$ GeV should be greater than the inclusive cross section for $`\mathrm{{\rm Y}}(1S)`$ by a factor of about 4.3. If we assume that the $`p_T`$ distributions have the same shapes at smaller $`p_T`$, we can estimate the cross section for $`\eta _b(1S)`$ integrated over all $`p_T`$. The resulting cross section is large enough that it should be possible to discover the $`\eta _b(1S)`$ at the Tevatron if there is a decay mode with a large enough branching fraction that can be used as a trigger. We argued that the decay $`\eta _bJ/\psi +J/\psi `$ should allow the discovery of the $`\eta _b(1S)`$, if not in the data from Run I, then certainly in Run II of the Tevatron.

###### Acknowledgements.

This work was supported in part by the U.S. Department of Energy under Grants DE-FG02-91-ER40690 and DOE-ER-40682-143 and by the Natural Sciences and Engineering Research Council of Canada. We thank F. Maltoni for valuable comments.