{'fraction_non_alphanumeric': 0.03479987813272249, 'fraction_numerical': 0.015515509867864275, 'mean_word_length': 4.633231629799136, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quarkonium Production at High-Energy Colliders

The theoretical description of heavy quarkonium production at high-energy p-pbar and e-p colliders is reviewed. Predictions based on non-relativistic QCD factorisation are confronted with recent charmonium and bottomonium data from the Tevatron and HERA. Potential shortcomings of the present theoretical analyses are discussed, and the prospects for quarkonium physics at the upgraded Tevatron and HERA colliders and at the LHC are summarised.


Introduction
The production of heavy quarkonium states at high-energy colliders has been the subject of considerable interest during the past few years. A wealth of new experimental results has become available, some of which revealed dramatic shortcomings of earlier quarkonium models. In theory, progress has been made on the factorisation between the short-distance physics of heavy-quark creation and the long-distance physics of bound state formation. The colour-singlet model [1,2] has been superseded by a consistent and rigorous framework, based on the use of non-relativistic QCD (NRQCD) [3], an effective field theory that includes the so-called colour-octet mechanisms. On the other hand, the colour-evaporation model [4][5][6] of the early days of quarkonium physics has been revived [7][8][9][10].
However, despite the recent theoretical and experimental developments the range of applicability of the different approaches is still subject to debate, as is the quantitative verification of factorisation. Because the quarkonium mass is still not very large with respect to the QCD scale, in particular for the charmonium system, nonfactorisable corrections [11][12][13] may not be suppressed enough, if the quarkonium is not part of an isolated jet, and the expansions in NRQCD may not converge very well. In this situation, a global analysis of various processes is mandatory to assess the importance of different quarkonium production mechanisms, as well as the limitations of a particular theoretical framework.
Much of the recent interest and theoretical development in quarkonium physics has been stimulated by experimental data from high-energy pp and ep colliders. The purpose of this review is to summarise the theoretical progress in understanding heavy quarkonium production, and to confront the predictions with recent charmonium and bottomonium data from the Tevatron and HERA. Forthcoming experimental results and theoretical advances, for example in the calculation of higher-order corrections, will sharpen the present picture and allow more quantitative tests of the theory. Therefore, this review will mostly focus on the generic features of the various quarkonium production processes and highlight the prospects for quarkonium physics at future collider experiments. Quarkonium production has also been studied in various other processes like e + e − annihilation, Z decays, hadronic collisions at fixed-target experiments and B decays, and several excellent reviews have been presented in the literature which cover these aspects, see e.g. [14][15][16][17][18].
The NRQCD factorisation approach to quarkonium production is introduced in Section 2 and compared to earlier quarkonium models like the colour-singlet model and the colour-evaporation model. Section 3 summarises the NRQCD predictions for charmonium production and polarisation at hadron colliders. The predictions are confronted with recent experimental data from the Tevatron, and potential shortcomings of the theoretical analysis are discussed. Charmonium production at HERA is reviewed in Section 4, and the prospects for quarkonium physics at the upgraded HERA collider are summarised. Section 5 concludes with an outlook on future work that will be required to test the theoretical framework more rigorously.

Quarkonium Production Mechanism
The production of heavy quarkonium at high-energy colliders involves two distinct scales. The creation of a heavy quark pair is a short-distance process on scales of the order 1/m Q or smaller and can be calculated in perturbation theory. The subsequent non-perturbative transition from the intermediate QQ pair to a physical quarkonium, on the other hand, involves long-distance scales of the order of the quarkonium size 1/(m Q v) or larger. Provided that these two scales are well separated, 1/(m Q v) ≫ 1/m Q , the formation of the bound state should be insensitive to the details of the heavy quark creation process, which is essentially local on the scale of the quarkonium size. It is thus intuitive to expect that long-distance and shortdistance physics in quarkonium production can be separated such that binding effects factorise into universal non-perturbative parameters.

Non-Relativistic QCD
The factorisation approach to quarkonium production [3] provides a systematic framework for separating physics at long-distance and short-distance scales. It is based on the use of non-relativistic QCD (NRQCD) [19], an effective field theory in which the heavy quark and antiquark are treated non-relativistically. A non-relativistic approach to quarkonium physics is motivated by the success of potential models in describing static properties of charmonia and bottomonia. The effective field theory framework of NRQCD allows a rigorous analysis of quarkonium physics which can be used to systematically improve earlier model approaches.
The most general effective Lagrangian of NRQCD takes the form where ψ and χ are two-component spinor fields describing the heavy quark and antiquark, respectively, and L light is the usual relativistic QCD Lagrangian for gluons and light quarks. The covariant derivative is D µ = ∂ µ + igA µ , where A µ = (A 0 , A) is the SU(3) gauge field and g is the QCD coupling. The term δL in (1) includes all possible operators consistent with the symmetries of QCD and reproduces the relativistic effects of the full theory. The leading corrections in δL are bilinear in the quark (or antiquark) field: where E i = G 0i and B i = 1 2 ǫ ijk G jk are the chromoelectric and chromomagnetic components of the gluon field-strength tensor G µν . The spin symmetry of the minimal NRQCD Lagrangian is broken by the two terms in (2) that contain the Pauli matrix σ. The creation and annihilation of heavy quark pairs is taken into account by local four-fermion operators in δL: The factors K i , K ′ i are products of a spin and a colour matrix and may also contain polynomials in the spatial derivate D. The operators in δL 4−fermion annihilate a QQ pair in a colour and angular momentum state specified by K ′ i and create a QQ pair at the same space-time point in a state specified by K i . The dimensionless coefficients c i , d i in the effective Lagrangian are sensitive only to short distances of order 1/m Q or smaller and can be determined by matching scattering amplitudes in NRQCD with the corresponding amplitudes in full QCD.
NRQCD calculations can, in principle, reproduce QCD results for quarkonium observables to any desired accuracy by adding correction terms to the effective Lagrangian δL and tuning the couplings to appropriate values.

NRQCD factorisation
Within the framework of NRQCD, the cross section for producing a quarkonium state H can be expressed as a sum of terms, each of which factors into a short-distance coefficient and a long-distance matrix element [3]: The sum includes all colour and angular momentum states of the QQ pair, denoted collectively by n. The short-distance coefficients dσ are proportional to the cross sections for producing a QQ pair in the state n and with small relative momentum. They can be calculated perturbatively in the strong coupling α s . 2 The non-perturbative transition probabilities from the QQ state n into the quarkonium H are given by vacuum expectation values of local four-fermion operators in NRQCD. Their generic form is where the sum is over the quarkonium polarisation λ and any number of light hadrons X in the final state. The factors K n , K ′ n can contain products of colour and spin matrices as well as covariant derivatives. The colour and angular momentum quantum numbers n of the intermediate QQ pair need not equal those of the physical quarkonium H. Soft gluons with energies of order m Q v or smaller can be emitted before the bound state is formed and change the colour and spin of the heavy quark pair. The effects of these soft gluons are included in the long-distance matrix elements The derivation of the factorisation formula (4) relies on the general space-time structure of heavy quark production to separate the short-distance effects of the hard momentum scale ∼ m Q from the long-distance effects of the soft momentum scale ∼ Λ QCD . This step, which is referred to as topological factorisation, uses standard factorisation methods of perturbative QCD [20]. The additional step required to derive the factorisation formula for the production of non-relativistic bound states involves separating the effect of the QQ relative momentum ∼ m Q v from the hard momentum scale m Q . This can be achieved by expanding the hard scattering part of the amplitude as a Taylor series in terms of the heavy quark relative momentum. The long-distance factors generated by the Taylor expansion can be identified with the vacuum expectation values of local four-fermion operators in NRQCD. In the resulting expression (4) all effects of the momentum scale m Q v or smaller are contained in the long-distance matrix elements O H [n] . They have to be calculated using nonperturbative methods or determined from experimental data. The short-distance coefficients dσ, on the other hand, involve only momentum scales of order m Q or larger and therefore have a perturbative expansion in powers of α s (m Q ). A general prescription for calculating the short-distance coefficients dσ can be found in [21].
The factorisation of quarkonium cross sections (4) has not been proven with complete rigour. Factorisation implies, for example, that soft gluons connecting the quarkonium and remnant jet parts of the scattering amplitude cancel for inclusive processes, but no explicit proof has been presented so far. Although it should be pos- 2 In the case of hadronic production, the short-distance cross sections have to be convoluted with parton distribution functions. 3 The dependence of the long-distance matrix elements and short-distance coefficients on the NRQCD factorisation scale has not been indicated and cancels in their product (4).
sible to extend the existing factorisation theorems of perturbative QCD to inclusive quarkonium processes, the applicability of the NRQCD factorisation approach has not yet been established theoretically in all circumstances.
The general factorisation formula (4) contains an infinite number of terms. It provides a viable phenomenological framework only if the relative magnitude of the various non-perturbative transition probabilities can be determined systematically. NRQCD predicts a hierarchy among the matrix elements O H [n] in terms of their dependence on the intrinsic heavy-quark velocity v. Using NRQCD power counting rules, the general expression (4) can be organised into an expansion in powers of v.
To any given order in v only a finite number of terms contribute to the production cross section.

Power counting
The importance of the various long-distance matrix elements O H [n] can be assessed using power counting rules. They depend upon the relative size of the three different low-energy scales in a non-relativistic bound state: the typical heavy quark threemomentum m Q v (in the quarkonium rest frame), the typical kinetic energy m Q v 2 and Λ QCD . 4 The standard NRQCD power counting (velocity scaling) rules [22] have been derived from basic qualitative properties of non-relativistic bound states and from analysing the equations of motion for the NRQCD quantum field operators. For example, the number of heavy quarks in the quarkonium, d 3 r ψ † ψ ≈ 1 implies ψ ∼ (m Q v) 3/2 since the quarkonium state can be localised within a region r ∼ 1/(m Q v). Similarly, the kinetic energy operator d 3 r ψ † D 2 /(2m Q ) ψ has an expectation value of m Q v 2 , and so the spatial part of the covariant derivative has to scale like D ∼ m Q v. From the equation of motion for the vector potential one finds g E ∼ m 2 Q v 3 and g B ∼ m 2 Q v 4 . According to these rules, the terms in δL bilinear (2) are suppressed by O(v 2 ) relative to those of the leading NRQCD Lagrangian. 5 The power counting rules for the long-distance matrix elements O H [n] in (4) can be derived by considering the Fock state decomposition of a quarkonium state 4 The power counting scheme for effective theories describing heavy-light mesons is more transparent. Because in a heavy-light meson the three-momentum and the energy of the heavy quark are both of order Λ QCD , the expansion parameter is Λ QCD /m Q and the importance of the operators in the effective Lagrangian is determined by their dimensions. 5 The contributions from δL 4−fermion (3) contain extra suppression factors O(α 2 s ) from the shortdistance coefficients d i .
|H in Coulomb gauge, The dominant component |QQ comprises a heavy quark pair in a colour-singlet state and with angular momentum quantum numbers 2S+1 L J that are consistent with the quantum numbers of the physical quarkonium. The higher Fock states, such as |QQg , contain dynamical gluons or light qq pairs. The heavy quark pair can be in either a colour-singlet or a colour-octet state with spin S = 0, 1 and angular momentum L = 0, 1, 2, . . .
Such a hierarchy is likely to be realised for bottomonium states where m Q ≈ 5 GeV and v 2 ≈ 0.1. It is, however, not obvious a priori that it can be applied to the charmonium system where m Q ≈ 1.5 GeV and v 2 ≈ 0.25, so that m Q v ≈ 750 MeV. Depending on the relation between the low energy scales m Q v, m Q v 2 and Λ QCD , alternative power counting schemes may be more appropriate for charmonium. If, for example, m Q v ≈ Λ QCD ≫ m Q v 2 the typical momentum of dynamical gluons would be ∼ Λ QCD , and chromomagnetic transitions may no longer be suppressed with respect to chromoelectric transitions [16]. Alternative power counting rules can have important consequences for quarkonium phenomenology [23][24][25], as discussed in more detail in the subsequent sections.

Colour-octet processes
The NRQCD formalism implies that so-called colour-octet processes associated with higher Fock state components of the quarkonium wave function must contribute to the cross section. Heavy quark pairs that are produced at short distances in a colouroctet state can evolve into a physical quarkonium through radiation of soft gluons at late times in the production process, when the quark pair has already expanded to the quarkonium size. Such a possibility is ignored in the colour-singlet model [1,2], where one assumes that only heavy quark pairs produced in the dominant Fock state form a physical quarkonium.
The most profound theoretical evidence for the incompleteness of the coloursinglet model comes from the presence of infrared divergences in the production cross sections and decay rates of P -wave quarkonium. Consider for example the production of χ J states. Within the colour-singlet model, χ J is described by a colour-singlet QQ pair with angular momentum quantum numbers 3 P J . The cross section is assumed to factorise into a short-distance coefficient, which describes the production of a QQ [1, 3 P J ] pair, and a single long-distance factor which contains all the nonperturbative dynamics of the bound state formation: The long-distance factor is related to the derivative of the radial wave function at the origin and can be determined from potential model calculations or from the χ J decay width. At next-to-leading order, the short-distance cross sections for production (or decay) of P -wave quarkonia exhibit a logarithmic infrared divergence due to soft gluon emission, which can not be factored into |R ′ (0)| 2 . In phenomenological applications of the colour-singlet model, an infrared cutoff has been introduced and identified with one of the low-energy scales present in non-relativistic bound states. It is evident, however, that the presence of infrared divergencies implies a failure of the simple factorisation assumption upon which the colour-singlet model is based.
The NRQCD approach provides a natural solution to this problem [26]. According to the NRQCD power counting rules, two terms contribute to the production cross section of P -wave states at leading order in the velocity expansion. For χ J production one has The matrix element O χ J [1, 3 P J ] can be related to the derivative of the radial wave function at the origin, and the first term in (8) corresponds to the expression of the colour-singlet model. The second term represents a contribution to the cross section from a colour-octet mechanism. The short-distance factor is the cross section for producing a QQ pair in a colour-octet 3 S 1 state, and the corresponding long-distance matrix element O χ J [8, 3 S 1 ] describes the probability for such a QQ [8, 3 S 1 ] pair to form a χ J quarkonium. The colour-octet 3 S 1 state can be reached from the dominant χ J Fock state through a single chromoelectric transition. According to the NRQCD power counting rules, the O χ J [1, 3 P J ] and O χ J [8, 3 S 1 ] matrix elements scale with the same power of velocity ∼ v 5 and must both be included in a consistent theoretical analysis. Other terms in the general factorisation formula are suppressed by relative order v 2 or more. In the NRQCD analysis, the infrared divergence in the coloursinglet short-distance cross section dσ(QQ [1, 3 P J ] + X) at next-to-leading order is cancelled by a matching infrared singularity from the radiative corrections to the colour-octet matrix element O χ J [8, 3 S 1 ] . The inclusion of the colour-octet process proportional to O χ J [8, 3 S 1 ] is crucial to remove the infrared divergence from the production cross section. 6 Colour-octet contributions are needed for a consistent description of P -wave quarkonia, but they can be even more important phenomenologically for S-wave states like J/ψ or ψ(2S). According to the power counting rules, all colour-octet matrix elements for the production (or decay) of S-wave quarkonia are suppressed by powers of the velocity compared to the colour-singlet contribution. The leading long-distance matrix element for J/ψ production, for example, involves the operator O J/ψ [1, 3 S 1 ], which annihilates and creates a cc pair with the same colour and angular momentum quantum numbers as in the dominant Fock state of the J/ψ. From the power counting rules it follows that the corresponding vacuum expectation value O J/ψ [1, 3 S 1 ] scales as v 3 . According to the spin and angular momentum selection rules, the cc pair can reach the colour-octet states cc [ 11 or smaller. In the non-relativistic limit v → 0, the NRQCD description of S-wave quarkonium production thus reduces to the colour-singlet model. However, as discussed in detail in the subsequent sections, colour-octet processes can become significant, and even dominant, if the short-distance cross section for producing QQ in a colour-octet state is enhanced.
NRQCD power counting rules predict a hierarchy among the various transition probabilities O H [n] . Truncating the expansion (4) at a given order in the heavyquark velocity, only a finite number of terms contribute to the production cross section. If the operator O H [n] creates and annihilates a QQ pair in a colour-singlet state with the angular momentum quantum numbers of the dominant quarkonium Fock state, its vacuum expectation value O H [n] can be related to the quarkonium wave function at the origin and determined from the leptonic decay width or potential model calculations. The number of independent non-perturbative matrix elements can be reduced further by using the approximate heavy-quark spin symmetry of the NRQCD Lagrangian [3]. The calculation of the remaining NRQCD matrix elements that contribute at a given order in velocity expansion requires non-perturbative methods. In the absence of lattice estimates 7 , the transition probabilities O H [n] have to be treated as phenomenological parameters and determined from experimental data, see Section 3.2.
Colour-octet processes are taken into account in the so-called colour-evaporation model [4][5][6], where the cross section for a specific quarkonium state H is given as a universal fraction f H of the inclusive QQ production cross section integrated up to the open heavy quark threshold: The model assumes that any QQ pair, produced with small relative momentum, can form a physical quarkonium by emission of soft gluons. No constraints are imposed on the colour and angular momentum configuration. The cross section is hence dominated by QQ pairs produced at short distances with vanishing angular momentum and in either a colour-singlet or colour-octet state.
The colour-evaporation model resembles some features of the NRQCD factorisation approach, in particular due to the inclusion of colour-octet processes. However, the colour-evaporation model is not a rigorous theory, and there is no limit in which it could systematically reproduce a full QCD calculation. Moreover, it turns out that the assumption of a single universal long-distance factor is too restrictive from a phenomenological point of view. It implies, for example, a universal σ(χ c )/σ(J/ψ) ratio, in conflict with experimental observations on charmonium production in hadronhadron and photon-hadron collisions.

Quarkonium Production at Hadron Colliders
The first hadron collider measurements of inclusive J/ψ, ψ(2S) and χ c cross sections from the UA1 experiment [30] at the CERN SppS and from the CDF collaboration [31] at the Fermilab Tevatron could not separate charmonium produced directly in a hard-scattering reaction from charmonium produced in weak decays of B-hadrons. The early data showed some discrepancies when compared to the coloursinglet model, which predicts B-decays to be the dominant source of charmonium at hadron colliders [32]. The substantial uncertainties in the theoretical calculation, however, made it difficult to draw definite conclusions [33]. A rigorous test of the colour-singlet model became possible when CDF presented a measurement of the direct ψ(2S) cross section [34][35][36], where the contributions from B-decays had been removed using microvertex-detection. It turned out that only approximately 25% of ψ(2S) at the Tevatron come from the decay of B-hadrons. Surprisingly, the direct ψ(2S) cross section at p t (ψ) ∼ > 5 GeV proved to be orders of magnitudes larger than the leading-order colour-singlet model prediction [37,38]. A similarly striking discrepancy was observed in the direct J/ψ cross section after contributions from both B-decays and radiative χ c decays had been removed. 8 The cross section for quarkonium hadroproduction is given by where f i/p and f j/p denote the parton densities, and ij = gg, gq, gq, qq. In the coloursinglet model, the leading partonic subprocess g + g → QQ [1, 3 S 1 ] + g occurs at order α 3 s [37,38]. In retrospect, the shortcoming of the colour-singlet model to describe direct J/ψ and ψ(2S) hadroproduction can be understood by examining a typical leading-order Feynman diagram, Figure 1(a). At large transverse momentum, the two internal quark propagators are off-shell by ∼ p 2 t so that the parton differential cross section scales like dσ/dp 2 t ∼ 1/p 8 t , as indicated in the figure. On the other hand, when p t ≫ 2m c the quarkonium mass can be considered small and the inclusive charmonium cross section should scale like any other single-particle inclusive cross section ∼ 1/p 4 t . The dominant production mechanism for charmonium at sufficiently large p t must thus be via fragmentation [39], the production of a parton with large p t which subsequently decays into charmonium and other partons. A typical fragmentation contribution to colour-singlet J/ψ production is shown in Figure 1(b). While the fragmentation contributions are of higher order in α s compared to the fusion process Figure 1(a), they are enhanced by a power p 4 t /(2m c ) 4 at large p t and can thus overtake the fusion contribution at p t ≫ 2m c . When colour-singlet fragmentation is included, the p t dependence of the theoretical prediction is in agreement with the Tevatron data but the normalisation is still underestimated by about an order of magnitude [40][41][42], indicating that an additional fragmentation contribution is still missing. It is now generally believed that gluon fragmentation into colour-octet 3 S 1 charm quark pairs [43,44], as shown in Figure 1(c), is the dominant source of J/ψ and ψ(2S) at large p t at the Tevatron. The probability of forming a J/ψ particle from a pointlike cc pair in a colour-octet 3 S 1 state is given by the NRQCD matrix element O J/ψ [8, 3 S 1 ] which is suppressed by v 4 relative to the non-perturbative factor of the leading coloursinglet term. However, this suppression is more than compensated by the gain in two powers of α s /π in the short-distance cross section for producing colour-octet 3 S 1 charm quark pairs, as compared to colour-singlet fragmentation. At O(v 4 ) in the velocity expansion, two additional colour-octet channels have to be included, Figure 1(d), which do not have a fragmentation interpretation at order α 3 s but which become significant at moderate p t ∼ 2m c [45,46]. The importance of the cc[8, 1 S 0 ] and cc[8, 3 P J ] contributions can not be estimated from naive power counting in α s and v alone, but rather follows from the dominance of t-channel gluon exchange, forbidden in the leading-order colour-singlet cross section. The leading-order subchannels which contribute to J/ψ and ψ(2S) hadroproduction at order v 4 are thus Explicit expressions for the parton cross sections can be found in [46] and [47].

Charmonium cross sections
The different contributions to the J/ψ transverse momentum distribution are compared to the CDF data [36] in Figure 2. As mentioned above, the colour-singlet model at lowest order in α s fails dramatically when confronted with the experimental results. When colour-singlet fragmentation is included, the prediction increases by more than an order of magnitude at large p t , but it still falls below the data by a factor of ∼ 30. The CDF results on J/ψ production can be explained by including the leading colour-octet contributions cc [8, 1 S 0 ], cc [8, 3 S 1 ] and cc [8, 3 P J ], and adjusting the corresponding non-perturbative parameters to fit the data. 9 Numerically one finds the non-perturbative matrix elements to be of O(10 −2 GeV 3 ), see Table 1, consistent with the v 4 suppression expected from the velocity scaling rules. Similar conclusions can be drawn for ψ(2S) production, see Figure 3. Table 1: NRQCD matrix elements for charmonium production. The colour-singlet matrix elements are taken from the potential model calculation of [48,49]. The colour-octet matrix elements have been extracted from the CDF data [36], The errors quoted are statistical only. Parameters: CTEQ5L parton distribution functions [50], renormalisation and factorisation scale µ = p 2 t + 4m 2 c and m c = 1.5 GeV. The Altarelli-Parisi evolution has been included for the cc[8, 3 S 1 ] fragmentation contribution. See [51] for further details.
Colour-octet processes need to be included in a theoretically consistent description of χ c production. According to Eq.(8) both cc [1, 3 P J ] and cc [8, 3 S 1 ] intermediate states contribute to the χ c cross section at leading order in the velocity expansion. While the discrepancy between the colour-singlet prediction and the experimental results [52] is less dramatic than in the case of S-wave charmonia, the inclusion of the 3 S 1 colour-octet process (with O χc [8, 3 S 1 ] adjusted to the value listed in Table 1) significantly improves the theoretical description. The comparison with recent data is presented in Figure 4. The measured ratio of production cross sections σ(χ c2 )/σ(χ c1 ) = 0.96 ± 0.27(stat) ± 0.11(sys) [53] is also in good agreement with the NRQCD expectation of 1.1 ± 0.2 [54].
The analysis of the CDF cross section data alone, although very encouraging, does not provide a conclusive test of NRQCD factorisation because free parameters have to be introduced to fit the data. However, if factorisation holds the non-perturbative total colour-octet 1 S 0 + 3 P J colour-octet 3 S 1 LO colour-singlet colour-singlet frag.

Figure 2:
Colour-singlet and colour-octet contributions to direct J/ψ production in pp → J/ψ + X at the Tevatron ( √ s = 1.8 TeV, pseudorapidity cut |η| < 0.6)) compared to experimental data from CDF [36]. Parameters: CTEQ5L parton distribution functions [50]; factorisation and renormalisation scale µ = p 2 t + 4m 2 c ; m c = 1.5 GeV. The leading logarithms (α s ln p 2 t /(2m c ) 2 ) n have been summed by solving the Altarelli-Parisi evolution equations for the gluon fragmentation function. NRQCD matrix elements as specified in Table 1. matrix elements, Table 1, are universal and can be used to make predictions for various processes and observables. The polarisation signature of J/ψ and ψ(2S) at large transverse momentum is one of the single most crucial tests of the NRQCD approach and will be discussed in some detail in Section 3.3. Ultimately, a global analysis of different production processes will be needed to assess the universality of the colour-octet contributions. Such a program, however, requires a careful discussion of the theoretical uncertainties associated with the determination of the non-perturbative matrix elements from experimental data. For charmonium hadroproduction the uncertainties are mainly associated with the calculation of the short-distance cross sections and, as argued in Section 3.2, have not yet been reliably quantified.
The J/ψ and ψ(2S) transverse momentum distributions at the Tevatron can also be described by the colour-evaporation model [7][8][9] or the related soft-colour interaction model [10]. Both approaches allow colour-octet cc pairs produced in gluon     Figure 2 for χ c production. Experimental data from [52]. fragmentation to hadronise into charmonium and thus resemble the NRQCD cross section prediction at large p t . The colour-evaporation model, however, assumes unsuppressed gluon emission from the cc pair during hadronisation, which randomises spin and colour, and consequently predicts unpolarised quarkonium. The NRQCD approach, in contrast, incorporates the spin-symmetry of QCD in the heavy quark limit, which implies transverse polarisation of direct J/ψ and ψ(2S) at large p t . The polarisation prediction of NRQCD will be discussed in detail in Section 3.3.

Determination of NRQCD matrix elements
The numerical values for the non-perturbative matrix elements obtained from the fit to the Tevatron data are sensitive to any effects that modify the shape of the charmonium p t distribution. This is true in particular for the O ψ [8, 1 S 0 ] and O ψ [8, 3 P 0 ] matrix elements which are determined from data at relatively small p t ∼ < 8 GeV. The theoretical uncertainties include the small-x behaviour of the gluon distribution [51], the evolution of the strong coupling [51], as well as systematic effects inherent in NRQCD, such as the inaccurate treatment of the energy conservation in the hadronisation of the colour-octet cc pairs [55]. Moreover, higher-order QCD corrections are expected to play an important role, as discussed in more detail below.
The inclusion of higher-order QCD corrections is required to reduce the theoretical uncertainty and to allow a more precise prediction of the charmonium cross sections. Next-to-leading order (NLO) calculations for quarkonium production at hadron colliders are presently available only for total cross sections [56,57]. Significant higher-order corrections to differential distributions are expected from the strong renormalisation and factorisation scale dependence of the leading-order results [51]. Moreover, the NLO colour-singlet cross section includes processes like g + g → QQ[1, 3 S 1 ] + g + g which are dominated by t-channel gluon exchange and scale as ∼ α 4 s (2m Q ) 2 /p 6 t . At p t ≫ 2m Q their contribution is enhanced with respect to the the leading-order cross section, Figure 1(a), which scales as ∼ α 3 s (2m Q ) 4 /p 8 t . This is born out by preliminary studies [58] which include part of the NLO hadroproduction cross section and by the complete calculation of NLO corrections to the related process of quarkonium photoproduction [59]. The NLO colour-singlet cross section may be comparable in size to the colour-octet 1 S 0 and 3 P J processes, which scale as ∼ α 3 , and affect the determination of the corresponding NRQCD matrix elements from the Tevatron data. A full NLO analysis is however needed before quantitative conclusions can be drawn.
Another source of potentially large higher-order corrections is the multiple emission of soft or almost collinear gluons from the initial state partons. These corrections, as well as effects related to intrinsic transverse momentum, are expected to modify the shape of the transverse momentum distribution predominantly at relatively low values of p t ∼ < 2m Q . Initial state radiation can be partially summed in perturbation theory [60], but so far only total cross sections have been considered in the literature [61]. An estimate of the effect on the transverse momentum distribution should be provided by phenomenological models where a Gaussian k t -smearing is added to the initial state partons. The result of these calculations, however, not only depends on the average k t , which enters as a free parameter, but also on the details of how the smearing is implemented. Moreover, a lower cut-off has to be provided which regulates the divergences at p t = 0. 10 Qualitatively, the inclusion of a Gaussian k t -smearing leads to an enhancement of the short distance cross section at small p t , which results in smaller values for the fits of the O ψ [8, 1 S 0 ] and O ψ [8, 3 P 0 ] NRQCD matrix elements [62,58]. The actual size of the effect, however, turns out to be different for the two models studied in the literature.
An alternative approach to treat the effect of initial state radiation is by means of Monte Carlo event generators which include multiple gluon emission in the parton shower approximation. Comprehensive phenomenological analyses have been carried out for charmonium production at the Tevatron and at the LHC [63][64][65] using the event generator Pythia [66] supplemented by the leading colour-octet processes [63]. The inclusion of initial state radiation as implemented in Pythia leads to an enhancement of the short-distance cross section. The size of the effect is significantly larger than for the Gaussian k t -smearing mentioned above, and it extents out to large p t . Consequently, the O ψ [8, 1 S 0 ] and O ψ [8, 3 P 0 ] NRQCD matrix elements estimated from the Monte Carlo analysis of the Tevatron cross sections are significantly lower than those listed in Table 1.
Recently, the impact of a non-vanishing initial state transverse momentum has been analysed in the context of the k t -factorisation formalism. This approach has been developed to sum large logarithmic corrections in the high-energy limit s ≫ 4m 2 c . The summation generates initial state gluons with transverse momenta comparable to the hard scale m Q . They are described by so-called unintegrated gluon distribution functions. The k t -factorisation calculations [67,68] indicate that the colour-singlet contribution to the direct J/ψ cross section could be substantially larger than in the conventional collinear factorisation approach. Moreover, the non-zero transverse momentum of the initial state partons allows for the production of high-p t colour-octet cc pairs without emission of hard gluons. The shape of the transverse momentum distribution for the various colour-octet processes in the k t -factorisation formalism is thus different from the curves displayed in Figure 2. The analyses [67,68] suggest that the direct J/ψ cross section might be dominated by the cc [8, 1 S 0 ] and cc [8, 3 P J ] channels up to large transverse momenta of the order p t ∼ < 20 GeV. However, more work is needed to establish the quantitative predictions of the k t -factorisation approach in general. At present there is no reliable estimate of the corresponding theoretical uncertainties due to higher-order corrections and, in particular, due to the poorly known unintegrated gluon distribution. The impact of non-vanishing initial state transverse momentum has to be quantified in various other high-energy QCD processes before the importance of the different quarkonium production mechanisms can be assessed conclusively in this framework.
To illustrate the uncertainty in extracting the NRQCD colour-octet matrix elements from the J/ψ, ψ(2S) and χ c transverse momentum distributions at the Tevatron, the results of various numerical analyses are collected in Tables 2, 3 and 4. It is evident from the comparison that the theoretical error is not very well under control. The discrepancies between the analyses that rely on the leading-order collinear factorisation calculation can be attributed to different choices for the input parameters and different prescriptions for the interpolation between the fusion cross section at low p t and the fragmentation cross section including summation of ln(p 2 t /(4m 2 c ) terms at large p t . More accurate data on the J/ψ and ψ(2S) cross sections could help to reduce some of the ambiguities. Ultimately, the calculation of next-to-leading order corrections to both colour-singlet and colour-octet processes is needed to improve the analysis and reduce the large factorisation and renormalisation scale uncertainty of the leading-order prediction. Effects beyond fixed order in perturbation theory, like the summation of soft gluon radiation, have to be addressed in a more rigorous and systematic way, and a partial summation of the NRQCD expansion might be needed to account for the inaccurate treatment of the kinematics in the hadronisation of the colour-octet cc pairs.
Taking into account the large theoretical uncertainties of present analyses, most results appear to be consistent with the NRQCD power counting rules. The colouroctet matrix elements extracted from the k t -factorisation calculation, however, do not follow the expected pattern. Such extreme values for the NRQCD matrix elements would indicate a serious flaw in our understanding of the power counting.
It is instructive to compare the determination of the colour-octet matrix elements from the Tevatron data with analyses of other production processes. NRQCD matrix elements have been extracted from charmonium hadroproduction at fixed-target experiments and from B-decays. The calculations include O(α s ) corrections to the short-distance coefficients and, in the case of B-decays, ≈ 0.04 ± 0.01 ≈ 6.5 ± 0.5 5 Table 2: NRQCD matrix elements for J/ψ production obtained from various analyses of the J/ψ transverse momentum distribution at the Tevatron. Values are given in units The first error quoted is statistical, the second error, when listed, due to variation of the renormalisation and factorisation scales between µ = 1/2 p 2 t + 4m 2 c and 2 p 2 t + 4m 2 c .

Reference
Potentially large higher-twist effects at small p t may render the application of the NRQCD factorisation formula to quarkonium hadroproduction at fixed-target experiments unreliable. Nevertheless, the values obtained for the non-perturbative matrix elements are in qualitative agreement with the determination at the Tevatron. Bdecays provide a more accurate testing ground for NRQCD factorisation. The linear combination of O ψ [8, 1 S 0 ] and O ψ [8, 3 P 0 ] extracted from the branching ratios BR(B → ψ + X) is consistent with Tevatron analyses that include soft-gluon effects in the parton-shower approximation.
In subsequent numerical analysis the following ranges of values for the non-perturbative matrix elements will be considered: Note that although the theoretical errors seem conservative, they do not accommodate the non-perturbative matrix elements obtained from calculations in the k t factorisation approach. 11

J/ψ and ψ(2S) polarisation
The analysis of J/ψ and ψ(2S) polarisation at large transverse momentum is one of the most decisive tests of the NRQCD factorisation approach. At large p t , the di-rect J/ψ and ψ(2S) cross sections should be dominated by gluon fragmentation into colour-octet 3 S 1 charm quark pairs, Figure 1(c). When p t ≫ 2m c the fragmenting gluons are effectively on-shell and the intermediate heavy quark pairs transversely polarised. According to the standard NRQCD power counting rules, Section 2.3, cc [8, 3 S 1 ] evolves into J/ψ or ψ(2S) predominantly through a double chromoelectric transition which preserves the heavy quark spin. The dominance of gluon fragmentation into colour-octet 3 S 1 charm quark pairs at large p t and NRQCD spin-symmetry thus imply transversely polarised J/ψ and ψ(2S) at large transverse momentum [80].
The polarisation can be measured through the angular distribution in the leptonic decay, where θ denotes the angle between the lepton three-momentum in the ψ rest frame and the ψ three-momentum in the laboratory frame. The polar-angle asymmetry α can take values in the range −1 ≤ α ≤ 1, where α = 1 corresponds to transverse polarisation.
Corrections to the asymptotic NRQCD prediction of pure transverse polarisation can arise from spin-symmetry breaking chromomagnetic interactions in the nonperturbative evolution of the colour-octet 3 S 1 charm quark pair. NRQCD power counting rules imply that spin-symmetry breaking transitions, which lead to random polarisation, are suppressed by a factor v 4 . Hence their effect on α should not exceed ≈ 15% [81]. Additional corrections to transverse polarisation result from higher-order terms in the short-distance part of the gluon fragmentation function. After hard gluon radiation, the intermediate heavy quark pair can be in various spin and angular momentum states and hadronise into unpolarised or longitudinal charmonium. The explicit calculation [81] reveals that the effect is numerically quite small and that radiative O(α s ) corrections can decrease the transverse polarisation by only about 10%. The dominant source of depolarisation comes from the colour-octet 1 S 0 and 3 P J fusion processes, Figure 1(d). They are important for moderate transverse momentum and need to be included to predict quarkonium polarisation in the p t -range probed by present data [51,82].
At O(v 4 ) in the velocity expansion, the polarisation yield from all production channels can be calculated unambiguously in NRQCD in terms of the non-perturbative matrix elements that have been determined from the unpolarised cross section. 12  (8, 1 S 0 , 3 P 0 ) have been considered. No transverse polarisation is expected at p t ∼ 5 GeV, but the angular distribution is predicted to change drastically as p t increases. A first measurement from CDF [83] does not support this prediction, but the experimental errors are too large to draw definite conclusions.
The analysis of J/ψ polarisation is complicated by the fact that the experimental data sample [83] includes J/ψ that have not been produced directly but come from decays of χ c and ψ(2S) mesons. The contribution from radiative χ decays decreases the transverse J/ψ polarisation at large p t , but tends to be cancelled out by the contribution from ψ(2S) decays [69,84]. The overall polarisation yield for J/ψ at large transverse momentum is thus not strongly diluted by decays of higher charmonium states. Again, the first experimental data do not confirm the NRQCD prediction [69,84], see Figure 6.
The absence of transversely polarized J/ψ and ψ(2S) at large transverse momentum, if confirmed with higher-statistics data at the Tevatron Run II, represents a serious challenge for the application of NRQCD to charmonium production. In the following, some uncertainties and potential shortcomings of the present theoretical analysis will be addressed: The uncertainty in the determination of NRQCD matrix elements O ψ [8, translates into an uncertainty in the predicted yield of transversely polarized J/ψ and ψ(2S). In the large-p t limit, the 3 S 1 production channel yields transverse polarisation, while the 1 S 0 and 3 P J channels both yield unpolarised quarkonia in this limit. A decreased value of the O ψ [8, 3 S 1 ] matrix element, relative to O ψ [8, 1 S 0 , 3 P J ] , would thus delay the onset of transverse polarisation and reduce the discrepancy between NRQCD predictions and data. However, extreme values for the non-perturbative matrix elements outside the range specified in Equation 14 would be in conflict with most of the analyses listed in Tables 2 and  3, and moreover be inconsistent with conventional NRQCD power counting.
Higher-order QCD effects in the short-distance cross section should be included to improve the accuracy of the theoretical prediction, in particular in the intermediate p t region. Preliminary studies indicate large NLO effects for the unpolarized p tdistributions, see Section 3.1, but no attempt has been made so far to estimate the impact of higher-order corrections on the spin-dependence of the cross section. First studies within the k t -factorisation approach [85] suggest that the J/ψ and ψ(2S) cross  Figure 5: Polar angle asymmetry α for ψ(2S) production in pp → ψ(2S)(→ µ + µ − ) + X at the Tevatron as a function of p t compared to experimental data from CDF [83]. Parameter specifications as in Figure 2. NLO corrections to the fragmentation contribution [81,51] have been included. The error band is obtained as a combination of the uncertainty (statistical only) in the extraction of the NRQCD matrix elements [ section might not be dominated by gluon fragmentation but rather by the colour-octet 1 S 0 process up to p t ∼ < 20 GeV. Since the 1 S 0 production channel yields unpolarised quarkonia, the overall J/ψ and ψ(2S) polarisation might be greatly reduced. However, as mentioned previously, the theoretical uncertainties in the k t -factorisation prediction need yet to be quantified. 13 Heavy quark spin-symmetry is violated by higher-order terms in the NRQCD Lagrangian, and longitudinal polarisation can arise if the binding of the charm quark pair into J/ψ and ψ(2S) proceeds through two chromomagnetic transitions. The suppression of the spin-symmetry breaking chromomagnetic interactions follows from the conventional NRQCD power counting rules. These rules, however, have not been firmly established for the case of charmonium. Alternative power counting schemes [16,24] imply that the chromomagnetic transitions might not be suppressed with respect to chromoelectric transitions. As a consequence, the O ψ [8, 1 S 0 ] matrix element would ′ Figure 6: Polar angle asymmetry α for prompt J/ψ production (i.e. not from B-decay) in pp → J/ψ(→ µ + µ − ) + X at the Tevatron as a function of p t compared to experimental data from CDF [83]. From [69]. Finally, charmonium production mechanisms beyond NRQCD could be responsible for J/ψ and ψ(2S) depolarisation. For example, the leading-twist formalism of NRQCD factorisation does not include possible rescattering interactions between the intermediate heavy quark pair and a comoving colour field. It has been shown that such rescattering corrections could yield unpolarized quarkonium [87,88]. The analysis [87,88], however, relies on several simplifying assumptions, and further work is needed to establish the importance of comover interactions for charmonium production at large p t .
The distinctive NRQCD prediction of transverse J/ψ and ψ(2S) polarisation at large p t is not supported by first experimental data. However, the measurements still have rather large errors and no conclusions should be drawn until the statistics improves. If, in Run II of the Tevatron, the experimental results continue to dis-agree with NRQCD predictions, one might have to conclude that either higher-order QCD effects could be crucial to describe the the spin-dependence of the charm cross section, or that alternative power counting schemes might be needed for the case of charmonium.

Bottomonium production
The application of NRQCD should be on safer grounds for the bottomonium system. As v 2 ∼ 0.1 for bottomonium, higher-order terms in the velocity expansion (in particular colour-octet contributions) are expected to be less significant than in the case of charmonium. Cross sections for the production of Υ and χ b states have been measured at the Tevatron in the region p t ∼ < 20 GeV [89][90][91]. The leading-order colour-singlet model predictions underestimate the data, but the discrepancy is much less significant than in the case of charmonium. Given the large theoretical uncertainties in the cross section calculation, in particular at small p t ∼ < M Υ , the need for colour-octet contributions is not yet as firmly established as for charmonium production. The inclusion of both next-to-leading order corrections and the summation of soft gluon radiation is required to obtain a realistic description of the Υ cross section in the p t -range probed by present data. Such calculations have not yet been performed.
The inclusive cross section for Υ(1S) production at the Tevatron is compared to the leading-order NRQCD calculation in Figure 7. The theoretical curves are based on the simple choice of non-perturbative parameters listed in Table 5 [92]. The colour-octet matrix elements have been determined from the CDF data for p t > 8 GeV, where the contributions from direct Υ(1S) production and from radiative χ b decays are measured separately. Given the large uncertainty in the leading-order calculation of the short distance cross section, no attempt has been made to quantify the error of the non-perturbative parameters listed in Table 5. The ambiguities in the determination of the non-perturbative matrix elements for bottomonium production are best demonstrated by comparing the results in Table 5 with other numerical analyses [46,93], collected in Table 6 for the case of direct Υ(1S) production.
The impact of initial state gluon radiation on the Υ cross sections at the Tevatron has been estimated by adding a Gaussian k t -smearing as discussed previously in the context of charmonium production. An average k t ∼ 3 GeV and a K-factor ∼ 3 are found to bring the leading-order colour-singlet cross section in line with the inclusive Υ(1S, 2S) data [94]. Similar results have been obtained within a Monte Carlo analysis [95], leading to significantly lower fit values for the colour-octet NRQCD matrix elements than those determined from a leading-order calculation. Moreover, the Monte Carlo results imply that no feeddown from χ states produced through  to experimental data from CDF [36]. Parameters: CTEQ5L parton distribution functions [50]; factorisation and renormalisation scale µ = p 2 t + 4m 2 b ; m b = 4.88 GeV. NRQCD matrix elements as specified in Table 5. colour-octet 3 S 1 bb states is needed to describe the inclusive Υ cross section, in contrast to what is found at leading-order. The calculation of next-to-leading order corrections and a systematic treatment of soft gluon radiation within perturbation theory are required to resolve these issues.
A preliminary measurement shows that inclusive Υ(1S) is produced unpolarised within the experimentally accessible range of transverse momentum. This result is consistent with NRQCD factorisation which predicts the onset of transverse Υ polarisation at p t (Υ) ∼ > 15 GeV [96], beyond what is probed by present data. 14 14 Transverse polarisation of Υ states produced in p-Cu collisions at √ s = 38.8 GeV has been reported recently [97]. The observation of a non-zero polarisation disproves the assumption that gluon emission from the bb pair during hadronisation randomises spin and colour, as implicit in the colour-evaporation model. On the other hand, the degree of transverse polarisation observed experimentally is smaller than predicted by a NRQCD calculation [98,99].  13.6 ± 6.8 +10.8 −7.5 Table 6: NRQCD colour-octet matrix elements for Υ(1S) production obtained from various analyses of the transverse momentum distribution at the Tevatron. Values are given in units 10 −2 GeV 3 . The first error quoted is statistical, the second error, when listed, due to variation of the renormalisation and factorisation scales between µ = 1/2 p 2 t + m 2 b and 2 p 2 t + m 2 b .

Prospects for the Tevatron Run II and the LHC
Run II at the Tevatron will provide a substantial increase in luminosity and will allow to determine the J/ψ, ψ(2S) and χ c cross sections more precisely and up to larger values of p t . A accurate measurement of the J/ψ and ψ(2S) polarisation at large transverse momentum will be the most crucial test of NRQCD factorisation. In addition, improved data on the J/ψ and ψ(2S) cross sections will help to reduce some of the ambiguities in extracting the colour-octet matrix elements.
With increased statistics it might be possible to access other charmonium states like the η c (nS) or the h c (nP ). Heavy-quark spin symmetry provides approximate relations between the non-perturbative matrix elements that describe spin-singlet and spin-triplet states [3]. The matrix elements for η c (nS) are related to those for ψ(nS), while the leading matrix elements for h c (nP ) can be obtained from those for χ c (nP ). Within NRQCD the rates for η(nS) and h(nP ) can thus be predicted unambiguously in terms of the non-perturbative matrix elements that describe the J/ψ, ψ(2S) and χ c cross sections. Comparing the various charmonium production rates would hence provide a stringent test of NRQCD factorisation and heavy-quark spin symmetry. The cross sections for producing η c and h c at Run II of the Tevatron are large [100,101], but the acceptances and efficiencies for observing the decay modes that can be triggered upon are in general small, and detailed experimental studies are needed to quantify the prospects. Other charmonium processes that have been studied in the literature include the production of D wave states [102], J/ψ production in association with photons [103,104] and double gluon fragmentation to J/ψ pairs [105].
The larger statistics expected at Run II of the Tevatron will also allow to improve the measurement of the bottomonium cross sections. Yet undiscovered states like η b (1S) could be detected, for example in the decay η b → J/ψ + J/ψ [93], and the associated production of Υ and electroweak bosons might be accessible [106]. If sufficient statistics can be accumulated, the onset of transverse Υ(nS) polarisation may be visible at p t (Υ) ∼ > 15 GeV.
A comprehensive test of the NRQCD factorisation approach in the bottomonium sector and a measurement of the Υ polarisation at p t ≫ M Υ will have to wait until the LHC starts operating. By that time, new experimental data from the Tevatron and HERA, as well as theoretical progress, e.g. in the calculation of higher-order corrections, will have significantly improved the present picture. Still, in order to provide some benchmark cross sections for quarkonium production at the LHC, the NRQCD predictions for selected charmonium and bottomonium states have been collected in Figure 8 [92]. Given the substantial uncertainty in the determination of the non-perturbative matrix elements from present data, the predictions should be regarded as order-of-magnitude estimates. 15 The polarisation pattern predicted for direct Υ(1S) production at the LHC is presented in Figure 9, based on the NRQCD matrix elements of Table 5. Higher-order corrections to the gluon fragmentation function [81] will lead to a small reduction of the transverse polarisation at large p t and should be included once data become  available. If the charmonium mass is indeed not large enough for a nonrelativistic expansion to be reliable, the onset of transverse Υ polarisation at p t ≫ M Υ may become the most decisive test of the NRQCD factorisation approach.  Table 5, other parameters as in Figure 7. The error band reflects the limiting cases that either O ψ [8, 1 S 0 ] or O ψ [8, 3 P 0 ] is set to zero in the linear combination extracted from the data.

Quarkonium production at ep colliders
The analysis of charmonium and bottomonium cross sections at the high-energy ep collider HERA provides a powerful tool to assess the importance of the different quarkonium production mechanisms and to test the general picture developed in the context of NRQCD factorisation. Section 4.1 focusses on J/ψ and ψ(2S) photoproduction, for which most of the experimental data have been collected. The impact of colour-octet processes and higher-order QCD corrections on cross sections and various differential distributions is reviewed, and the theoretical predictions are confronted with recent experimental results. Quarkonium production in deep-inelastic scattering is discussed in Section 4.2. The high-statistics data expected at the upgraded HERA collider will allow detailed studies of more exclusive observables like J/ψ polarisation, see Section 4.3. Moreover, various other processes may be accessible, including χ c production, associated J/ψ + γ production and the production of bottomonium states. The prospects for quarkonium physics at the upgraded HERA collider will be discussed in some detail in Section 4.4.

J/ψ and ψ(2S) photoproduction
Quarkonium production in high-energy ep collisions at HERA is dominated by photoproduction events where the electron is scattered by a small angle producing photons of almost zero virtuality. Experimental results have been reported for J/ψ and ψ(2S) photoproduction in a wide kinematical region, 30 GeV ∼ < s γp ∼ < 200 GeV. The cross section is dominated by photon-gluon fusion, where the photon interacts as a pointlike particle. Besides the direct production channel, charmonium photoproduction at HERA can also take place through resolved reactions, in which case the photon participates in the hard scattering through its parton content. Resolved processes contribute significantly to the lower endpoint of the charmonium energy spectrum and have been probed at HERA for the first time recently. Higher-twist phenomena like elastic/diffractive quarkonium production [108][109][110] are not included in the leading-twist calculation of NRQCD (see [47] for a detailed discussion) and have to be eliminated by either a cut in the J/ψ transverse momentum p t ∼ > 1 GeV, or by a cut in the J/ψ energy variable z ≡ p P · k ψ /p P · k γ ∼ < 0.9 (in the proton rest frame, z is the ratio of the J/ψ to γ energy, z = E ψ /E γ ).

Photoproduction mechanism
The cross section for quarkonium photoproduction can be written as where f i/γ and f j/p denote the parton distributions in the photon and the proton, respectively. For the direct photon processes i = γ and f γ/γ (x) = δ(1 − x). In the following, it is assumed that the J/ψ transverse momentum is restricted to p t > 1 GeV, in order to suppress the elastic contribution and higher-twist corrections in general. Generic Feynman diagrams which contribute to the parton cross section in the inelastic region away from p t = 0 (or z = 1) are shown in Figure 10. At leading order in the velocity expansion, J/ψ is produced through the colour-singlet channel [1], Figure 10(a). Relativistic corrections at O(v 2 ) modify the large-z and small-p t region [111][112][113][114] but can be neglected for p t ∼ > 1 GeV. The next-to-leading order QCD corrections to the direct photon process [115,59] increase the cross section prediction significantly, in particular at p t ∼ > m c as discussed in detail below.
At O(v 4 ) relative to the colour-singlet contribution, the J/ψ can also be produced through intermediate colour-octet 3 S 1 , 1 S 0 and 3 P J configurations [116][117][118], see Figure 10(b). Next-to-leading order corrections to the colour-octet channels are only known for the total photoproduction cross section (integrated over all z and p t ) [119]. The importance of the colour-octet production processes follows from a kinematical enhancement of the short-distance cross section, similar to the case of J/ψ hadroproduction. The 1 S 0 and 3 P J configurations can be produced through t-channel exchange of a gluon already at lowest order in α s . (For the 3 S 1 octet channel this is true for the resolved process.) This leads to a potentially significant colour-octet contribution, in particular in the large-z region.
The leading-order subchannels which contribute to inelastic J/ψ photoproduction through direct photon processes are γ + g → cc 1, 3 S 1 ; 8, 3 S 1 ; 8, 1 S 0 ; 8, 3 P J + g, (17) γ + q/q → cc 8, 3 S 1 ; 8, 1 S 0 ; 8, where the initial-state parton originates from the target proton. For resolved photon processes, the subchannels are g + q/q → cc 8, 3 S 1 ; 8, 1 S 0 ; 8, q +q → cc 8, 3 S 1 ; 8, 1 S 0 ; 8, 3 P J + g, where one of the initial-state partons originates from the photon and the other originates from the proton. Explicit expressions for the parton differential cross sections are given in [47]. The direct-photon processes dominate in the region z ∼ > 0.2, whereas resolved-photon processes become important in the region z ∼ < 0.2. At HERA energies, photon-quark fusion can contribute about 10%-15% to the cross section at large z. Quark-gluon fusion constitutes about 20%-40% of the resolved cross section at z ∼ < 0.2 and becomes more important than gluon-gluon fusion at larger z. Quark-antiquark fusion is always completely negligible.
Fragmentation contributions in photon-gluon fusion [120][121][122] exist only at the next order in α s and are suppressed for p t ∼ < 10 GeV. The cc[8, 3 S 1 ]-channel is insignificant for the direct-photon contribution, but dominates in resolved-photon interactions at p t ∼ > 5 GeV because it includes a gluon fragmentation component. The resolved photon amplitudes are identical to those relevant to J/ψ production in hadron-hadron collisions, and at HERA energies the relative importance of the various contributions as functions of p t is nearly the same as at Tevatron energies.
The J/ψ energy distribution The J/ψ energy spectrum is a particularly sensitive probe of the different production mechanisms. The colour-singlet and colour-octet contributions to the differential cross section through direct and resolved photon processes are presented in Figure 11 [116,47]. The results are shown as a function of the scaling variable z = p ψ ·p p /p γ ·p p at a typical HERA energy of √ s γp = 100 GeV in the inelastic region p t > 1 GeV. Adopting the NRQCD matrix elements as determined from the Tevatron J/ψ data (see Section 3.2, Equation 14), colour-octet processes are predicted to contribute significantly near the upper and lower endpoint of the energy spectrum. As mentioned previously, the colour-octet enhancement is a consequence of t-channel gluon exchange in the cc [  The theoretical predictions are compared to experimental results obtained by the H1 [124] and ZEUS [125] collaborations at HERA in Figure 11. The colour-octet enhancement in the region z ∼ > 0.6 is not supported by the experimental data. At first sight this could lead to rather stringent constraints on the colour-octet matrix elements O J/ψ [8, 1 S 0 ] and O J/ψ [8, 3 P 0 ] and might point towards an inconsistency with the values obtained for these matrix elements from other processes. However, the peaked shape of the z-distribution is derived neglecting the energy transfer of softparton emission in the non-perturbative transition of the colour-octet charm quark pairs. Such kinematic effects become important near the endpoint of the J/ψ energy spectrum, and a summation of singular higher-order terms in the velocity expansion is necessary to predict the shape of the z-distribution in the region z ∼ > 0.7 [55]. 16 A 16 The situation is analogous to the well-known problem of extracting the CKM matrix element quantitative study [126] reveals that the kinematic impact of soft-gluon emission in the hadronisation of colour-octet cc pairs at large z is not well under control, even after the summation of higher-order NRQCD corrections. A more stringent cut on the J/ψ transverse momentum is needed to remove the sensitivity of the prediction to processes which contribute in the elastic region at z ≈ 1. At present, no reliable constraints on the size of colour-octet contributions in J/ψ photoproduction can be derived from the J/ψ energy spectrum near the upper endpoint.
The low-z region is not expected to be sensitive to higher-order terms in the velocity expansion. Therefore, if the data could be extended to z ∼ < 0.3, an important resolved photon contribution should be visible [127,121,122], if the colour-octet matrix elements are not significantly smaller than extracted from the Tevatron data and expected from power counting rules. A preliminary measurement [128] indeed shows evidence for resolved processes, but the statistical error is still too large to establish the importance of colour-octet channels in the low-z region.
NLO colour-singlet contributions The current experimental results on J/ψ photoproduction are adequately described by the colour-singlet channel alone, once higher-order QCD corrections are included [115,59] and once the normalisation uncertainty due to the choice for the charm quark mass and the QCD coupling is taken into account. The NLO calculation for the J/ψ energy distribution is confronted with the experimental data in Figure 12. The shaded error band illustrates the large normalisation uncertainty of the theoretical prediction: the upper boundary corresponds to a charm quark mass m c = 1.3 GeV and α s (M Z ) = 0.1225, while the lower boundary has been obtained using m c = 1.5 GeV and α s (M Z ) = 0.1175. 17 From Figure 12 one can infer that the shape of the energy distribution in the colour-singlet channel is not very sensitive to the choice of input parameters. The theoretical prediction is in good agreement with the experimental results. A similar conclusion can be drawn from the analysis of the inelastic cross section as a function of the photon-proton energy, Figure 12.
The next-to-leading order QCD corrections are crucial to describe the shape of the J/ψ transverse momentum distribution. The NLO colour-singlet cross section includes processes like γ + g → QQ[1, 3 S 1 ] + g + g which are dominated by t-channel gluon exchange and scale as ∼ α 3 s m 2 Q /p 6 t . At p t ∼ > m Q their contribution is enhanced with respect to the leading-order cross section, which scales as ∼ α 2 s m 4 Q /p 8 t . The comparison with the experimental data, Figure 13, confirms the importance of the higher-order corrections.
|V ub | from the endpoint region of the lepton energy distribution in semileptonic B decay. 17 The variation of the cross section with α s has been studied in a consistent way by adopting the corresponding set of parton distribution functions [129].   Figure 13: Leading order and next-to-leading order colour-singlet prediction for the J/ψ transverse momentum spectrum dσ/dp 2 t at the photon-proton centre-of-mass energy √ s γp = 100 GeV compared to experimental data from H1 and ZEUS [124,125]. Parameters as in Figure 12.
The colour-singlet channel provides a successful description of J/ψ photoproduction as far as present data on energy and transverse momentum distributions is concerned. The substantial normalisation uncertainties of the theoretical prediction, as well as the inaccurate kinematical treatment of soft gluon emission in the hadronisation of colour-octet cc pairs, however, do not allow to place stringent constraints on the role of colour-octet processes. With the expected increase in luminosity at HERA, observables like J/ψ polarisation and various other quarkonium final states will provide more information on the importance of the different production mechanisms, see ψ(2S) photoproduction While a non-relativistic treatment may be appropriate to describe J/ψ photoproduction in many circumstances, significant relativistic corrections are expected for radial excitations like ψ(2S). In the static limit v → 0, the σ(ψ(2S))/σ(J/ψ) cross section ratio is given by the ratio of the NRQCD long-distance factors O ψ(2S) [1, 3 S 1 ] / O J/ψ [1, 3 S 1 ] to all orders in perturbation theory. The colour-singlet matrix elements can be related to the quarkonium decay widths, so that σ(ψ(2S))/σ(J/ψ) = Γ ψ(2S) ee / Γ J/ψ ee = 0.41±0.05 [131] in the static approximation. The non-relativistic estimate is not confirmed by recent HERA measurements which yield σ(ψ(2S))/σ(J/ψ) = 0.210 ± 0.058 [124] and σ(ψ(2S))/σ(J/ψ) = 0.242 ± 0.065 [132], respectively. Higher-order terms in the velocity expansion will modify the theoretical prediction, their impact is, however, hard to quantify. Colour-octet processes might increase or decrease the cross section ratio, depending on the value of the nonperturbative matrix elements. No quantitative prediction is possible given the present theoretical error on O ψ [8, 1 S 0 ] and O ψ [8, 3 P 0 ] . Another source of relativistic corrections originates from NRQCD operators which involve total derivatives on fermion bilinears [55,16]. Although these operators are formally suppressed by powers of v 2 relative to the leading colour-singlet contribution, they are crucial to reproduce the hadronic phase space constraints of the quarkonium production cross section. According to factorisation, all dependence on the specific quarkonium state resides within the NRQCD matrix elements. Since the short-distance coefficients are consequently expressed in terms of the heavy quark mass 2m Q rather than the quarkonium mass m H , the physical phase space restrictions for producing a specific bound state are not manifest at fixed order in the velocity expansion. While 2m Q − m H ∼ m Q v 2 is formally suppressed and numerically small for J/ψ, the difference between parton and hadron kinematics can lead to a large ambiguity in predicting the cross section of radial excitations like ψ(2S) where 2m c − m ψ ≈ 1 GeV. The ambiguity can be resolved by summing certain higher-order terms in the velocity expansion [55]. The summation, however, leads to a non-perturbative distribution function, which can not be calculated except within models. A rough estimate [59], obtained by assuming that the effective charm masses in the short-distance cross sections scale like the corresponding ψ(2S) and J/ψ masses, implies a significant phase space suppression factor of (m J/ψ /m ψ ) 3 ≈ 0.6, consistent with the experimental observation.

J/ψ production in deep-inelastic scattering
Leptoproduction of quarkonium has been extensively studied in the context of the colour-singlet model [133][134][135][136] and in the framework of the NRQCD factorisation approach including colour-octet processes [137,138]. At large photon virtualities Q 2 , the theoretical prediction is under better control than in the case of photoproduction because higher-order corrections and higher-twist contributions should become less important as Q 2 increases. Moreover, diffractive processes [110] are expected to be suppressed at Q 2 ≫ 4m 2 Q , and the total quarkonium leptoproduction cross section can be sensibly compared with the NRQCD prediction, contrary to photoproduction. The leading-order J/ψ production process is O(αα s ) and proceeds purely through intermediate 1 S 0 and 3 P J colour-octet configurations. The leptoproduction cross section in the high-energy limit Q 2 , s ≫ 4m 2 c [137], and O J/ψ [8, 3 P 0 ] would also help to clarify the role of colour-octet processes in J/ψ photoproduction. At O(αα 2 s ) J/ψ is produced through 3 S 1 colour-singlet and 1 S 0 , 3 S 1 and 3 P J colour-octet configurations. Results for the parton cross sections can be found in [138]. The O(αα 2 s ) colour-octet processes are particularly important at large values of the inelasticity z. Note, however, that the upper endpoint region of the z-distribution in J/ψ leptoproduction can not be predicted reliably without summing large higher-order corrections in the NRQCD velocity expansion, similar to photoproduction.
Experimental results on J/ψ leptoproduction at HERA [139] are shown in Figure 14 compared to the leading-order NRQCD predictions. The theoretical cross sections include the O(αα 2 s ) colour-singlet process and the O(αα s ) colour-octet contributions. Taken at face value, the colour-singlet channel underestimates the data by a factor ∼ 2 − 3. There is, however, a considerable amount of uncertainty in the cross section normalisation, in particular from the value of the charm quark mass, the strong coupling, and the choice of renormalisation and factorisation scales. On the other hand, adopting O J/ψ [8, 1 S 0 ] = 0.01 GeV 3 and O J/ψ [8, 3 P 0 ] /m 2 c = 0.005 GeV 3 the sum of colour-singlet and colour-octet contributions exceeds the data by up to a factor of three. The shape of the distributions favour the prediction of the coloursinglet channel, but more experimental information extending to larger Q 2 values is needed before firm conclusions can be drawn.

J/ψ and ψ(2S) polarisation
Polarisation measurements are powerful to assess the relative importance of colouroctet processes in different kinematic regions. First preliminary results [140] indicate a non-zero polarisation signature in J/ψ photoproduction at HERA. With the expected increase in luminosity, the analysis of decay angular distributions could provide information on the relevance of the various quarkonium production processes, complementing the previously discussed polarisation measurement at hadron colliders. Figure 14: Colour-singlet and colour-octet contributions to inelastic J/ψ leptoproduction ep → e J/ψ X [137] compared to experimental data from H1 [139]. The kinematic region is 2 GeV 2 ≤ Q 2 ≤ 80 GeV 2 , 40 GeV ≤ W ≤ 180 GeV. The variable y * denotes the rapidity of the J/ψ in the γ * p centre-of-mass system. The invariant mass of the hadronic final state has been restricted to M X > 10 GeV to suppress diffractive contributions. Parameters: GRV(LO) parton distribution functions [72,123]; factorisation and renormalisation scales µ = Q 2 + 4m 2 c ; m c = 1. The general decay angular distribution can be parametrised as dΓ(J/ψ → l + l − ) dΩ ∝ 1 + λ cos 2 θ + µ sin 2θ cos φ + ν 2 sin 2 θ cos 2φ, where θ and φ refer to the polar and azimuthal angle of the l + three-momentum with respect to a coordinate system defined in the J/ψ rest frame. The parameters λ, µ, ν can be calculated within NRQCD as a function of the kinematic variables like z and p t [47]. Two examples are shown in Figure 15, the polar angle parameter λ and the azimuthal angle parameter ν as a function of the energy fraction z and the transverse momentum p t , respectively. Each plot exhibits the result from the colour-singlet channel alone (dashed line) and the result after including colour-octet contributions. Since the decay angular distribution parameters involve ratios of cross sections, the dependence on parameters that affect the absolute normalisation, such as the charm quark mass, strong coupling, the renormalisation/factorisation scale and parton distribution, cancels to a large extent and does not constitute a significant uncertainty. The z-dependence of the polar angle asymmetry (Figure 15, left plot) is distinctive and can be used to discriminate between NRQCD and the colour-evaporation model, which always predicts unpolarised J/ψ. 18 The transverse momentum distribution, on the other hand, can prove very useful to determine the relative magnitude of coloursinglet and colour-octet contributions: While colour-octet processes tend to predict unpolarised quarkonium in the p t region considered here, large and positive values of the azimuthal parameter ν are expected for p t ∼ > 5 GeV if the cross section is dominated by the colour-singlet channel, see Figure 15, right plot. To observe the unique transverse polarisation signature of gluon fragmentation into colour-octet cc pairs one would have to isolate the resolved photon contribution at small z ∼ < 0.3 and measure the polar angular distribution for large transverse momenta p t ∼ > 10 GeV. A comprehensive analysis of the p t and z dependence of the angular parameters λ, µ, ν for different reference frames can be found in [47].
Polarisation in J/ψ photoproduction through the colour-singlet channel has also been studied in the context of the k t -factorisation approach [141]. The virtuality of the initial state gluon increases the fraction of longitudinally polarised J/ψ, in particular at large transverse momentum. The effect becomes significant at p t ∼ > 3 GeV, and could be probed at the upgraded HERA collider by analysing the polar angle asymmetry. The impact of colour-octet processes on the polarisation prediction in the k t -factorisation approach has not been estimated so far.
Information on the quarkonium production mechanism can also be obtained from the angular distribution in inelastic J/ψ leptoproduction [135,137,138]. The theoretical prediction at leading order only depends on the ratio of the NRQCD colour-octet matrix elements O J/ψ [8, 1 S 0 ] and O J/ψ [8, 3 P 0 ] . Longitudinally polarised J/ψ is expected if the 3 P J channel dominates, while the 1 S 0 channel leads to unpolarised quarkonium. The experimental result, λ = 0.77 +0. 44 −0.38 [139], is consistent with coloursinglet dominance, which implies λ ≈ 0.5 [139]. The experimental uncertainties are, however, still too large to put stringent limits on the size of colour-octet contributions.

Prospects for the upgraded HERA collider
With the HERA luminosity upgrade, a wealth of new quarkonium data will become available. The existing measurements can be extended into new kinematic regions, and observables like J/ψ polarisation and various other quarkonium final states will be accessible. The future analyses of quarkonium production at HERA offer unique possibilities to test the theoretical framework of NRQCD factorisation and to clarify the role of colour-octet processes. Some of the most interesting and promising reactions will be discussed in the following.
the non-perturbative distribution functions that result from a summation of higher-oder corrections in the velocity expansion. Note, however, that the dependence is significant only if the angular parameter varies strongly in the endpoint region and disappears entirely if its distribution is flat [47].
Resolved photon processes As emphasised previously, the low-z region of the J/ψ energy spectrum in photoproduction is dominated by resolved photon processes, which resemble the J/ψ production mechanisms relevant at hadron colliders. Away from the upper endpoint region, the energy spectrum is not expected to be sensitive to higher-order terms in the velocity expansion, and a measurement of the resolved cross section can provide a sensible test of the importance of colour-octet processes.
Colour-singlet and colour-octet contributions to the low-z region of the energy spectrum are displayed in Figure 16, left, at a photon-proton energy of W γp = 170 GeV. Note that, because only a small fraction of the incoming photon energy participates in the hard interaction, the typical W γp values for resolved photon events within the detector acceptance are higher than for the medium-z analyses [128]. The right plot of Figure 16 shows the uncertainty of the leading-order colour-singlet cross section due to the variation of the charm quark mass and the renormalisation scale. Next-to-leading order corrections to inelastic J/ψ photoproduction through resolved photons are not yet available.
A preliminary measurement [128] clearly establishes the existence of resolved processes at small z, but experimental errors and theoretical normalisation uncertainties do not allow to constrain the size of colour-octet contributions. To confirm the significant excess over the colour-singlet process observed at the Tevatron it will be necessary to analyse resolved photon processes at larger values of p t . χ c photoproduction The measurement of inelastic χ c photoproduction is a particularly powerful way to discriminate between NRQCD and the colour-evaporation model. The assumption of a single universal long-distance factor in the colourevaporation model implies a universal σ(χ c )/σ(J/ψ) ratio, and a large χ c cross section is predicted in photon-proton collisions, similar to hadron colliders where σ(χ c )/σ(J/ψ) ≈ 0.5 [52]. In NRQCD, on the other hand, the σ(χ c )/σ(J/ψ) ratio is process-dependent and strongly suppressed in photoproduction. While χ c can be produced copiously in hadron-hadron collisions through gg, gq and qq interactions, the direct photon process γ + g → cc[1, 3 P J ] + g is forbidden at leading order in α s due to the specific colour structure and charge conjugation invariance. Apart from higher-order corrections and resolved photon interactions, χ c can thus only be produced through the colour-octet process γ + g → cc[8, 3 S 1 ] + g [142,127]. According to factorisation, and up to corrections of O(α s , v 2 ), the cross section ratio σ(χ c )/σ(J/ψ) in direct photon interactions can be expressed in terms of the corresponding longdistance factors with NRQCD matrix elements as listed in Table 1. A search for χ c production at HERA resulting in a cross section measurement or upper cross section limit would directly probe the colour-octet matrix element O χ J [8, 3 S 1 ] and test the assumption of a single universal long-distance factor implicit in the colour-evaporation model.
Photoproduction of spin-singlet states The inclusion of colour-octet processes is crucial to describe η c (nS) and h c (nP ) photoproduction. With regard to the P -wave state h c , the colour-octet contribution is required to cancel the infrared divergence present in the colour-singlet cross section [143]. The production of η c , on the other hand, is dominated by colour-octet processes since the colour-singlet cross section vanishes at leading-order due to charge conjugation invariance [144,145], similar to χ c photoproduction. As discussed previously in the context of future quarkonium physics at the Tevatron, Section 3.5, heavy-quark spin-symmetry allows to predict the production of spin-singlet states like η c (nS) and h c (nP ) in terms of the NRQCD matrix elements that describe J/ψ and χ c production. The cross sections for η c (nS) and h c (nP ) photoproduction are sizeable [143,144], but it is not obvious that these particles can be detected experimentally, even with high-statistics data.
Associated J/ψ + γ production The energy spectrum of J/ψ produced in association with a photon, γp → J/ψ + γ X, is a distinctive probe of colour-octet processes [146,127,147,148]. In the colour-singlet channel and at leading-order in α s , J/ψ + γ can only be produced through resolved photon interactions. The corresponding energy distribution is thus peaked at low values of z. The intermediate and large-z region of the energy spectrum is expected to be dominated by the colour-octet process γg → cc[8, 3 S 1 ] γ. Colour-singlet and colour-octet contributions to associated J/ψ + γ photoproduction through direct and resolved photon processes are shown in Figure 17 as a function of the J/ψ energy. Observation of a substantial fraction of J/ψ + γ events at z ∼ > 0.5 would provide clear evidence for the presence of colouroctet processes in quarkonium photoproduction. For a more comprehensive study of associated J/ψ + γ production including various other differential distributions see [148].
Bottomonium production With the significant increase in statistics at the upgraded HERA collider, it might be possible to study inelastic photoproduction of bottomonium states for the first time. The large value of the bottom quark mass makes the perturbative QCD predictions more reliable than for charm production, and the application of NRQCD should be on safer grounds for the bottomonium system where v 2 ≈ 0.1. However, the production rates for Υ states are suppressed, compared with J/ψ, by a more than two orders of magnitude at HERA, a consequence of the smaller bottom electric charge and the phase space reduction by the large b mass. The next-to-leading order colour-singlet cross section in the inelastic region p t > 1 GeV and z < 0.8 is shown in Figure 18 as a function of the photon-proton centre-of-mass energy. Several years of data taking at the upgraded HERA may be required to obtain accurate experimental information on the inelastic Υ photoproduction cross section.

Conclusion and Outlook
The NRQCD factorisation approach is a systematic framework for inclusive quarkonium production and provides a striking and elegant explanation for the large direct J/ψ and ψ(2S) cross sections measured in pp collisions at the Fermilab Tevatron. Still, more theoretical work and experimental information is needed to firmly establish the applicability of NRQCD factorisation to charmonium production. The absence of transversely polarised J/ψ and ψ(2S) at large transverse momentum, if confirmed with higher statistics data, could imply that the conventional NRQCD power counting rules have to be replaced by alternative schemes, or that higher-order QCD effects are essential to describe the spin-dependence of the charm cross section.
The substantial increase in luminosity at the Tevatron Run II will allow to assess the role of colour-octet contributions in bottomonium production more reliably. The calculation of next-to-leading order QCD corrections and a systematic treatment of soft-gluon radiation are, however, required to describe the cross section in the p t region probed by the experiments at the Tevatron. A comprehensive test of the NRQCD factorisation approach in the bottomonium sector will have to wait until the LHC starts operating. If the charmonium mass is indeed not large enough for a non-relativistic approach to be reliable, the onset of transverse Υ polarisation at p t ≫ M Υ may become the most decisive test of the NRQCD factorisation approach.
A global analysis of different production processes is needed to verify the universality of long-distance matrix elements in NRQCD. The colour-octet matrix elements obtained from Tevatron data are in qualitative agreement with analyses of other production processes, like charmonium hadroproduction at fixed-target experiments and charmonium production in B decays. However, the sizeable theoretical uncertainties, mainly related to QCD effects in the short-distance cross sections, prevent a more quantitative test at present.
The experimental data on charmonium photoproduction are adequately described by the colour-singlet channel, including next-to-leading order corrections in α s . But there is a considerable amount of uncertainty in the normalisation of the theoretical prediction, and no constraints on the size of colour-octet processes can be deduced from the photoproduction analysis at present. The upgraded HERA collider has a remarkable potential to test the theoretical framework more thoroughly. The analysis of resolved photon processes and charmonium polarisation, in particular, could be essential in assessing the importance of different quarkonium production mechanisms.
The theoretical issues to be addressed in the future include a critical evaluation of the power counting rules and the calculation of higher-order QCD corrections. The importance of non-vanishing initial state transverse momentum has to be clarified by analysing various high-energy QCD processes. Moreover, systematic effects inherent in NQRCD, such as the inaccurate treatment of energy conservation in the hadronisation of heavy quark pairs and the mismatch between partonic and hadronic phase space need to be addressed quantitatively. Finally, reliable estimates of nonfactorisable corrections and production mechanisms beyond NRQCD are required to improve the theoretical analyses. Theoretical progress and the wealth of new results from the future Tevatron and HERA experiments will lead to exciting and comprehensive phenomenology and result in a more complete and rigorous understanding of quarkonium physics.