{'fraction_non_alphanumeric': 0.042754668671638456, 'fraction_numerical': 0.01469183719087196, 'mean_word_length': 4.227533460803059, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
J/psi plus jet associated production in two-photon collisions at next-to-leading order

We calculate the cross section of J/psi plus jet inclusive production in gamma gamma collisions at next-to-leading order within the factorization formalism of nonrelativistic quantum chromodynamics (NRQCD) focusing on direct photoproduction. Apart from direct J/psi production, we also include the feed-down from directly-produced chi_{cJ} and psi' mesons. We discuss the analytical calculation, in particular the treatment of the various types of singularities and the NRQCD operator renormalization, in some detail. We present theoretical predictions for the future e^+e^- linear collider TESLA, taking into account both brems- and beamstrahlung.


Introduction
Since the discovery of the J/ψ meson in 1974, charmonium has provided a useful laboratory for quantitative tests of quantum chromodynamics (QCD) and, in particular, of the interplay of perturbative and nonperturbative phenomena. The factorization formalism of nonrelativistic QCD (NRQCD) [1,2] provides a rigorous theoretical framework for the description of heavy-quarkonium production and decay. This formalism implies a separation of short-distance coefficients, which can be calculated perturbatively as expansions in the strong-coupling constant α s , from long-distance matrix elements (MEs), which must be extracted from experiment. The relative importance of the latter can be estimated by means of velocity scaling rules; i.e., the MEs are predicted to scale with a definite power of the heavy-quark (Q) velocity v in the limit v ≪ 1. In this way, the theoretical predictions are organized as double expansions in α s and v. A crucial feature of this formalism is that it takes into account the complete structure of the QQ Fock space, which is spanned by the states n = 2S+1 L (a) J with definite spin S, orbital angular momentum L, total angular momentum J, and colour multiplicity a = 1, 8. In particular, this formalism predicts the existence of colour-octet (CO) processes in nature. This means that QQ pairs are produced at short distances in CO states and subsequently evolve into physical, coloursinglet (CS) quarkonia by the nonperturbative emission of soft gluons. In the limit v → 0, the traditional CS model (CSM) [3] is recovered. The greatest triumph of this formalism was that it was able to correctly describe [4] the cross section of inclusive charmonium hadroproduction measured in pp collisions at the Fermilab Tevatron [5], which had turned out to be more than one order of magnitude in excess of the theoretical prediction based on the CSM.
Apart from this phenomenological drawback, the CSM also suffers from severe conceptual problems indicating that it is incomplete. These include the presence of logarithmic infrared singularities in the O(α s ) corrections to P -wave decays to light hadrons and in the relativistic corrections to S-wave annihilation [6], and the lack of a general argument for its validity in higher orders of perturbation theory. While the k T -factorization [7] and hard-comover-scattering [8] approaches manage to bring the CSM prediction much closer to the Tevatron data, they do not cure the conceptual defects of the CSM. The colour evaporation model [9], which is intuitive and useful for qualitative studies, also significantly improves the description of the Tevatron data as compared to the CSM [10]. However, it does not account for the process-specific weights of the CS and CO contributions, but rather assumes a fixed ratio of 1 : 7. In this sense, a coequal alternative to the NRQCD factorization formalism is presently not available.
In order to convincingly establish the phenomenological significance of the CO processes, it is indispensable to identify them in other kinds of high-energy experiments as well. Studies of charmonium production in ep photoproduction [11], ep [12,13] and νN [14,15] deep-inelastic scattering (DIS), e + e − annihilation in the continuum [16], Z-boson decays [17], γγ collisions [18,19,20,21], and b-hadron decays [22] may be found in the literature; for reviews, see Ref. [23]. Furthermore, the polarization of ψ ′ mesons produced directly [24] and of J/ψ mesons produced promptly [25,26], i.e., either directly or via the feed-down from heavier charmonia, which also provides a sensitive probe of CO processes, was investigated. Until recently, none of these studies was able to prove or disprove the NRQCD factorization hypothesis. However, H1 data of ep → e + J/ψ + X in DIS at the DESY Hadron Electron Ring Accelerator (HERA) [27] and DELPHI data of γγ → J/ψ + X at the CERN Large Electron Positron Collider (LEP2) [28] provide first independent evidence for it [13,21].
The verification of the NRQCD factorization hypothesis is presently hampered both from the theoretical and experimental sides. On the one hand, the theoretical predictions to be compared with existing experimental data are, apart from very few exceptions [29,30,31], of lowest order (LO) and thus suffer from considerable uncertainties, mostly from the dependences on the renormalization and factorization scales and from the lack of information on the nonperturbative MEs. On the other hand, the experimental errors are still rather sizeable. The latter will be dramatically reduced with the upgrades of HERA (HERA II) and the Tevatron (Run II) and with the advent of CERN LHC and hopefully a future e + e − linear collider (LC) such as the TeV-Energy Superconducting Linear Accelerator (TESLA), which is presently being designed and planned at DESY. On the theoretical side, it is necessary to calculate the next-to-leading-order (NLO) corrections to the hard-scattering cross sections and to include the effective operators which are suppressed by higher powers in v.
In this paper, we take a first step in this comprehensive research programme, by studying the inclusive production of J/ψ mesons in high-energy γγ collisions. As mentioned above, this process was studied at LEP2 [28], where the photons originated from hard initial-state bremsstrahlung. At high-energy e + e − LCs, an additional source of hard photons is provided by beamstrahlung, the synchrotron radiation emitted by one of the colliding bunches in the field of the opposite bunch. The highest possible photon energies with large enough luminosity may be achieved by converting the e + e − LC into a γγ collider via back-scattering of high-energetic laser light off the electron and positron beams.
In order for a J/ψ meson to acquire finite transverse momentum (p T ), it must be produced together with another particle or a hadron jet (j). From coupling and phase-space considerations it is evident that J/ψ plus jet associated production yields the dominant contribution. In the following, we thus consider the process γγ → J/ψ + j + X, where X denotes the hadronic remnant possibly including a second jet. Here, we take j and X to be free of charm assuming that charmed hadrons or charmonia besides the J/ψ meson would be detectable. The process γγ → J/ψ + γ + X, where γ represents a prompt photon, will be considered in a forthcoming publication [32].
The incoming photons can interact either directly with the quarks participating in the hard-scattering process (direct photoproduction) or via their quark and gluon content (resolved photoproduction). Thus, the process γγ → J/ψ + j + X receives contributions from the direct, single-resolved, and double-resolved channels. All three contributions are formally of the same order in the perturbative expansion. This may be understood by observing that the parton density functions (PDFs) of the photon have a leading behaviour proportional to α ln(M 2 /Λ 2 QCD ) ∝ α/α s , where α is the fine-structure constant, M is the factorization scale, and Λ QCD is the asymptotic scale parameter of QCD. In the following, we focus our attention on the direct channel. The other channels are left for future work. The J/ψ mesons can be produced directly; or via radiative or hadronic decays of heavier charmonia, such as χ cJ and ψ ′ mesons; or via weak decays of b hadrons. The respective decay branching fractions are B(χ c0 → J/ψ + γ) = (1.02 ± 0.17)%, B(χ c1 → J/ψ + γ) = (31.6 ± 3.2)%, B(χ c2 → J/ψ + γ) = (18.7 ± 2.0)%, B(ψ ′ → J/ψ + X) = (55.7 ± 2.6)%, and B(B → J/ψ + X) = (1.15 ± 0.06)% [33]. The b hadrons can be detected by looking for displaced decay vertices with dedicated vertex detectors, and the J/ψ mesons originating from their decays can thus be treated separately. Therefore and because of the smallness of B(B → J/ψ + X), J/ψ production through b-hadron decay is not considered here. The cross sections of the four residual indirect production channels may be approximated by multiplying the direct-production cross sections of the respective intermediate charmonia with their decay branching fractions to J/ψ mesons.
To summarize, the goal of the present analysis is to calculate the inclusive cross section of γγ → J/ψ +X in direct photoproduction for finite values of p T at NLO within NRQCD allowing for the J/ψ meson to be promptly produced. The LO result, also including the single-and double-resolved contributions, may be found in Refs. [20,21] and the references cited therein. The leading relativistic correction, which originates from the cc Fock state n = 1 P and is of O(v 4 ) relative to the LO result, was evaluated in Ref. [19]. In that paper, also J/ψ plus dijet associated production in two-photon collisions was studied for direct photoproduction and direct J/ψ production. The 2 → 3 partonic cross sections obtained therein constitute a starting point for the evaluation of the real radiative corrections in our present study. They have to be complemented with the corresponding results for direct χ cJ production. Furthermore, the kinematic cuts imposed in Ref. [19] to separate the final-state objects have to be removed. This yields soft and collinear singularities, which are collectively denoted as infrared (IR) singularities. On top of this, the virtual radiative corrections, which arise from 2 → 2 Feynman diagrams 1 involving closed loops, have to be added. They involve IR, ultraviolet (UV), and Coulomb singularities. The cancellation of all these singularities is nontrivial and requires the UV renormalization of masses, couplings, and wave-functions; the redefinition of NRQCD MEs so as to absorb IR and Coulomb singularities; the factorization of initial-state collinear singularities in photon PDFs; and the operation of the Kinoshita-Lee-Nauenberg theorem [34] on the cancellation of final-state IR singularities. Apart from being of general phenomenological relevance, our analysis should thus also be of conceptual interest for the theoretical heavyquarkonium community. After all, this is the first time that the full NLO corrections are evaluated for an inclusive 2 → 2 process within the NRQCD framework.
This paper is organized as follows. In Section 2, we describe our analytical calculation in some detail. Specifically, we discuss the structure of the various types of singularities and the mechanisms by which they are removed. Lengthy expressions are relegated to the Appendix. In Section 3, we present our numerical results appropriate for the e + e − mode of TESLA, and discuss their phenomenological implications. Our conclusions are 1 Here, the cc bound state is considered as one particle.

Analytic results
We start this section with a few general remarks. In our analytic calculation, we take the colour gauge group to be SU(N c ) with a generic value of N c , which is put equal to 3 in our numerical analysis. Colour factors appearing in our formulas include T F = 1/2, We work in the fixedflavour-number scheme, with n f = 3 active quark flavours q = u, d, s, which we treat as massless. The charm quark c and antiquark c, with mass m, only appear in the final state. We denote the fractional electric charge of quark q by e q .
The cc Fock states contributing at LO in v are specified for H = J/ψ, χ cJ , ψ ′ in Table 1. Their MEs satisfy the multiplicity relations which follow to LO in v from heavy-quark spin symmetry. We employ dimensional regularization with d = 4−2ǫ space-time dimensions to handle the UV and IR singularities, and we introduce a 't Hooft mass µ and a factorization mass M as unphysical scales. We formally distinguish between UV and IR poles, which we denote as 1/ǫ UV and 1/ǫ IR , respectively. We apply the projection method of Refs. [31,35], which is equivalent to the d-dimensional matching procedure of Ref. [36], in order to extract the short-distance coefficients that multiply the MEs. However, in order to conform with common standards, we adopt the normalizations of the MEs from Ref. [2] rather than from Refs. [30,31]; i.e., the MEs include spin and colour average factors.
There is only one partonic subprocess at LO, namely which is mediated by the diagrams depicted in Fig. 1. The analogous process with n = 1 P (8) 1 yields a relativistic correction of O(v 4 ) and was studied in Ref. [19]. With the fourmomentum assignments indicated within the parentheses in Eq. (2), we have k 1 + k 2 = p + k 3 , k 2 1 = k 2 2 = k 2 3 = 0, and p 2 = 4m 2 . We define the Mandelstam variables as s = (k 1 + k 2 ) 2 , t = (k 1 − p) 2 , and u = (k 2 − p) 2 , so that s + t + u = 4m 2 . For convenience, we also introduce In the NLO analysis, we need to evaluate the cross section of process (2) in d dimensions and retain terms of O(ǫ) because UV counterterms appear in multiplicative renormalization. We have where the first factor on the right-hand side stems from the flux, T 0 is the LO transitionmatrix (T ) element of process (2), and it is averaged (summed) over the spin and colour states of the incoming (outgoing) particles. Here and in the following, we denote the Lorentz-invariant N-particle phase-space element in d dimensions as where p i and m i are the four-momenta and masses of the final-state particles and P is the total four-momentum of the incoming particles. In the case of process (2), we have where Γ(x) is Euler's Γ function, so that Note that the spin average factor of an incoming photon or gluon in d dimensions is 1/(d − 2). The final result is listed in Eq. (1) of Ref. [19]. We now move on to the process e + e − → e + e − H + X, where H denotes a generic charmonium state. For later use, we now also include the resolved contributions. We select the e + e − centre-of-momentum (CM) frame and denote the nominal e + e − energy by √ S, and the transverse momentum and rapidity of the H meson by p T and y, respectively. Invoking the Weizsäcker-Williams approximation (WWA) [37] and the factorization theorems of the QCD parton model [38] and NRQCD [2], the differential cross section of e + e − → e + e − H + X can be written as , where m T = p 2 T + 4m 2 is the transverse mass of the H meson. For a given value of √ S, the accessible phase space is defined by

Virtual corrections
Representative examples of diagrams that generate the virtual corrections to the cross section of process (2) are shown in Fig. 2. They fall into two classes. The diagrams of the first class are obtained by attaching one virtual gluon line in all possible ways to the treelevel seed diagrams of Fig. 1. Their contribution is proportional to e 2 c . They include selfenergy, triangle, box, and pentagon diagrams. Loop insertions in external gluon or charmquark lines are accommodated in the respective wave-function renormalization constants and are not displayed in Fig. 2. The self-energy and triangle diagrams are in general UV divergent; the triangle, box, and pentagon diagrams are in general IR divergent. The pentagon diagrams without three-gluon vertex also contain Coulomb singularities, which are cancelled after taking into account the corresponding corrections to the operator . In the practical calculation, the Coulomb singularities are first regularized by a small gluon mass. This regularization prescription is then transformed into one implemented with a small relative velocity v between the c and c quarks [29].
The diagrams of the second class contain a light-quark triangle or box. The triangle diagrams produce contributions proportional to e 2 q or e q e c , which vanish individually by Furry's theorem [39]. The box diagrams are proportional to e 2 q . They contain UV and IR singularities, but their sum is finite.
Calling the T -matrix element comprising all these loop diagrams T v , the virtual corrections to the cross section of process (2) may be evaluated as We apply two independent approaches to calculate the one-loop diagrams. The first one uses FeynArts [40] to generate the diagrams and self-written Mathematica codes to apply the projectors and provide expressions, which are afterwards treated with a FORM program to perform the tensor reduction and the extraction of the UV and IR singularities. The result is then transformed into a Fortran code to be used for the numerical evaluation. The second approach utilizes QGRAF [41] for the generation of the diagrams, FeynCalc [42] for the tensor reduction, and LoopTools [43] for the numerical evaluation of the IR-safe integrals. These packages are surrounded by self-written interface programs, which allow for a completely automated computation [44]. As in the first approach, the relevant IRdivergent n-point functions are implemented in a Fortran code. They can be extracted from Refs. [29,45,46]. For the analytical treatment of the abelian five-point functions, we refer to Ref. [46]. The pentagon diagrams involving a nonabelian coupling can all be reduced to integrals with a lower number of external legs by partial fractioning.
Our analytic result is too lengthy to be presented here. However, the interference of the light-quark box amplitude T q with T 0 is sufficiently compact to be listed, which is done in Eq. (A.1).

Parameter and wave-function renormalization
The self-energy and triangle diagrams of the type indicated in Fig. 2 contain UV singularities, which are cancelled upon the renormalizations of the QCD gauge coupling g s = √ 4πα s , the charm-quark mass m and field ψ, and the gluon field A µ . Specifically, the renormalization transformations read where the superscript 0 labels bare quantities and Z i = 1 + δZ i , with i = g, m, 2, 3, are renormalization constants. The quantities δZ i are of O(α s ) and contain UV singularities and, in general, also finite pieces. The UV singularities are unique, while the finite pieces depend on the choice of renormalization scheme. In order to comply with the Lehmann-Symanzik-Zimmermann [47] reduction formula without any additional finite adjustments, we define Z 2 and Z 3 in the on-mass-shell (OS) scheme, the results being where γ E is Euler's constant and β 0 = (11/3)C A − (4/3)T F n f is the one-loop coefficient of the QCD beta function. We also define m in the OS scheme, by setting However, we adopt the modified minimal-subtraction (MS) scheme to define g s , by putting

Operator renormalization
A crucial feature of effective field theories, such as NRQCD, is that the composite operators are generally subject to renormalization. In the case of NRQCD, this is essential in order to ensure the complete cancellation of IR and Coulomb singularities at NLO and so to overcome the conceptual problems of the CSM mentioned in Section 1. To be consistent with the rest of our calculation, we also employ dimensional regularization here. We adopt the technique described in Refs. [2,31] to directly evaluate the NLO corrections to the NRQCD operators. In this way, we avoid having to match partonic cross sections evaluated in NRQCD with their counterparts in full QCD.
In the case under consideration, we have to renormalize the CO ME O H 3 S , which appears at LO, in |T 0 | 2 . In d space-time dimensions, this ME has mass dimension d − 1. We thus introduce the 't Hooft mass scale of NRQCD, λ, to keep its renormalized version, which we wish to extract from experimental data, at mass dimension 3.
The four-quark operator is related to the amplitude for the elastic scattering of a cc pair. The corresponding tree-level diagram is depicted in Fig. 3(a). The one-loop corrections to this amplitude are obtained by attaching a virtual gluon line in all possible ways to the external heavy-quark legs, and they involve self-energy and vertex corrections [see Fig. 3 Using the NRQCD Feynman rules in the quarkonium rest frame, expanding the one-loop integrands as Taylor series in 1/m, and performing the integration over the loop momentum, we obtain the unrenormalized one-loop result where the subscript 0 labels the tree-level quantity and µ is the 't Hooft mass scale of QCD that enters through the d-dimensional loop integration. The term proportional to 1/v represents the Coulomb singularity, which arises from the exchange of a longitudinal gluon between the outgoing c and c quarks [see Fig. 3(c)]. Obviously, NRQCD operators of different cc Fock states n start to mix at one loop. Furthermore, the presence of UV singularities indicates that they need renormalization. In the following, we choose the MS scheme for that. We thus write where the subscript r labels the renormalized quantity and we identify λ with the NRQCD renormalization scale. Inserting Eq. (16) into Eq. (15), we obtain Combining Eq. (17) with Eq. (7) generates an IR counterterm at O(α s ) that is indispensible to render the total NLO result finite. The Coulomb singularity present in Eq. (17) is necessary to cancel similar terms in the virtual corrections.

Renormalization group equation
The renormalization group equation that determines the λ dependence of O H 3 S (8) 1 r (λ) may be derived by differentiating Eq. (17) with respect to λ and then taking the physical limit ǫ → 0. It reads There is no obvious physical reason to distinguish between the scales µ and λ, which refer to the ccg and cccc vertices in the one-loop diagrams of Figs. 3(b)-(e), respectively. Both scales should essentially be of O(m). In the following, we thus identify µ = λ. Integration of Eq. (18) then yields where λ 0 is a reference scale for which the value of O H 3 S (8) 1 r (λ 0 ) is assumed to be known and (m) with its λ-independent LO value, which is known from the literature [25]. Since we consider the NLO corrections to the LO cross section (7), which is already of O(α s ), we have to employ in Eq. (19) the two-loop formula for α s (λ), which reads where L = ln λ 2 /Λ 2 QCD .

Real corrections
The real corrections to the cross section of process (2) arise from the partonic subprocesses where n = 3 P (1) where n = 1 S J . The respective diagrams are presented in Figs. 4(a) and 5. From the technical point of view, it is convenient to evaluate the gluon polarization sum as pol ε a µ (q)ε b * ν (q) = −δ ab g µν , at the expense of allowing for the Faddeev-Popov ghosts of the gluon, u g and u g , to appear in the final state. The additional partonic subprocesses are where n = 3 S 1 . The corresponding diagrams are depicted in Fig. 4 1 is forbidden by Furry's theorem [39], as may be understood by observing that the 3 S (1) 1 projector effectively closes the charm-quark line and acts like a vector coupling and that the two gluons are then in a CS state, so that we are dealing with a closed fermion loop containing five vector couplings. This was also verified by explicit calculation. On the other hand, process (22) is prohibited by colour conservation, since the two quark lines are connected by a single gluon, which ensures that the cc and qq pairs are both in a CO state.
The diagrams in Fig. 5 can be divided into three classes. The diagrams of the first (second) class are proportional to e 2 q (e 2 c ) and contain a timelike virtual gluon that splits into a cc (qq) pair. The diagrams of the third class are proportional to e q e c and contain a spacelike virtual gluon.
As for the kinematics of processes (21)-(23), we now have k 1 + k 2 = p + k 3 + k 4 , k 2 1 = k 2 2 = k 2 3 = k 2 4 = 0, and p 2 = 4m 2 . We can define 5 2 = 10 Mandelstam variables, 5 of which are linearly independent. In addition to s 1 , t 1 , and u 1 defined in Eq. (3), we introduce If we take s 1 , t 1 , u 1 , s 14 , and s 24 to be the set of independent Mandelstam variables, then the others are given by Generically denoting the T -matrix elements of processes (21)-(23) by T r , their cross sections may be evaluated as In the case of process (21), a factor of 1/2 has to be included on the right-hand side of Eq. (26) to account for the fact that there are two identical particles in the final state. Integrating Eq. (26) over the three-particle phase space while keeping the value of p T finite, we encounter IR singularities, which can be of the soft and/or collinear type. In order to systematically extract these singularities in a Lorentz-invariant way, it is useful to slice the phase space by introducing infinitesimal dimensionless cut-off parameters δ i and δ f , which are connected with the initial and final states, respectively [48]. In the case of processes (21) and (23), we are thus led to distinguish between the following regions of phase space. In the soft region, where either s 3 + s 34 < δ f s or s 4 + s 34 < δ f s, one of the outgoing gluons is soft and, in addition, they may be collinear. In the final-state collinear region, where s 34 < δ f s < s 3 +s 34 , s 4 +s 34 , the outgoing gluons are both hard and collinear. In the hard region, which comprises the residual phase space, none of the outgoing gluons is soft and they are not collinear. In the case of process (22), we differentiate between the following regions of phase space. In the initial-state collinear region, where one of the inequalities s 13 < δ i s, s 14 < δ i s, s 23 < δ i s, and s 24 < δ i s is satisfied, one of the outgoing q and q quarks is collinear to one of the incoming photons. In the final-state collinear region, where s 34 < δ f s < s 3 + s 34 , s 4 + s 34 , the outgoing q and q quarks are both hard and collinear. In the hard region, which includes the left-over phase space, none of the outgoing q and q quarks is soft and they are not collinear. While the various contributions depend on the cut-off parameters δ i and/or δ f , their sum must be independent of them. The verification of the cut-off independence provides a nontrivial and powerful check for the correctness of our calculation.

Hard region
The integration of Eq. (26) over the hard region can be facilitated by decomposing the three-particle phase space as where k 34 = k 3 + k 4 and On the other hand, writing we have where the boundaries of integration are fixed by the kinematical constraints specified above. Since the contribution from the hard region is UV and IR finite, we may perform the phase-space integration in d = 4 space-time dimensions.

Soft region
As an example, let us consider the kinematical situation where k 3 in processes (21) or (23) becomes sufficiently soft, so that it does not influence the kinematics of the hard process involving the residual four-momenta. In practice, this is arranged by imposing the condition s 3 + s 34 < δ f s with a sufficiently small value of δ f . The three-particle phase space can then be decomposed as where At the same time, |T r | 2 factorizes as where F 3 is an appropriate Eikonal factor, the specific form of which depends on the cc Fock state n. Specifically, soft singularities occur for processes (21) and (23) with J . The appropriate expressions for F 3 are listed in Table 2. Notice that, in Eq. (33), T 0 always refers to process (2) with n = 3 S 1 . The integration of the Eikonal factors over the soft parts of the phase space can be performed analytically, and the results for the soft integrals, and similarly for k 4 , may be found in Table 3, where Process Limit Eikonal factor γγ → cc 3 P J .

Final-state collinear region
In the final-state collinear region, where s 34 < δ f s < s 3 + s 34 , s 4 + s 34 , it is useful to introduce the parameterization where z, with 0 < z < 1, determines the fractions of longitudinal momentum that the finalstate partons receive from the splitting gluon, k 2 T = z(1 − z)s 34 , and terms of O (k 4 T /p 4 ) are neglected. We can then decompose the three-particle phase space as where Using the relation

Process
Soft and collinear integrals where the terms of O(ǫ 2 ) are inconsequential in the physical limit d → 4 because, in the collinear limit, |T r | 2 only develops simple poles in ǫ, we recover the familiar combination of Euler's Gamma functions. At the same time, |T r | 2 factorizes as where F 34 is an appropriate eikonal factor. In the case under consideration, final-state collinear singularities occur for processes (21) and (22) with n = 3 S 1 , where they are related to g → gg and g → qq splitting, respectively. The corresponding expressions for F 34 may be found in Table 2. They are related to the timelike g → g and g → q splitting functions, as respectively. The integration of F 34 over the collinear parts of the phase space can be performed analytically, and the resulting collinear integrals, may be found in Table 3. In the case of process (21), the range of integration in Eq. (44) has to be further constrained by the condition δ f s/s 1 < z < 1 − δ f s/s 1 to ensure that the final-state gluons are both hard.

Initial-state collinear region
As an example, let us consider the kinematical situation where the incoming photon with four-momentum k 2 in process (22) splits into a hard, collinear, quasireal antiquark q, with four-momentum k ′ 2 , and the outgoing quark q, with four-momentum k 3 . The collinearity is enforced by the constraint s 23 < δ i s with a sufficiently small value of δ i . It is then helpful to introduce the parameterization where k 2 T = −(1 − z)s 23 and terms of O (k 4 T /k 4 2 ) are neglected. We can then decompose the three-particle phase space as where At the same time, |T r | 2 factorizes as where the cc Fock state n is the same as appears in T r and F 23 is an appropriate Eikonal factor. In the case under consideration, initial-state collinear singularities occur for process (22) with n = 1 S J , where they are always related to γ → qq splitting. The corresponding expressions for the Eikonal factors F ij , where i = 1, 2 and j = 3, 4, may be found in Table 2. They are related to the spacelike γ → q splitting function, where P g→q (z) is defined in Eq. (42), as The integration of Eq. (48) over the collinear region of the phase space, with s 23 < δ i s, can be performed analytically. In terms of differential cross sections, the result reads where s is the two-photon invariant mass square defined above Eq. (3). According to the mass factorization theorem [38], the form of this collinear singularity connected with an incoming-photon leg is universal and can be absorbed into the bare PDF of the antiquark q inside the resolved photon, f q/γ (z). To this end, one defines the renormalized photon PDF as where M is the factorization scale and the finite function C γ (z) is introduced to fix the factorization scheme. In the modified minimal-subtraction (MS) scheme, we have C γ (z) = 0, while in the so-called DIS γ scheme, which is used, e.g., in Ref. [49], we have Incorporating the second term on the right-hand side of Eq. (52) in Eq. (51), we are left with the finite contribution, which depends on the choice of M. In turn, the resolved-photon contribution is evaluated with the renormalized photon PDF f q/γ (z, M 2 ), which also depends on M. In the sum of these two contributions, this M dependence cancels up to terms beyond NLO. In this sense, the notion of direct-photon contribution ceases to be separately well defined at NLO.

Assembly of the NLO cross section
The NLO result for the cross section of process (2) is obtained by adding the virtual and real corrections to the LO result. Collecting the various contributions discussed above, arising from the parameter and wave-function renormalization (ct), the operator redefinition (op), the initial-state (is) and final-state (fs) collinear configurations, the softgluon radiation (so), and the hard-parton emission (ha), we can schematically write the resulting differential cross section of γγ → H + X, including the MEs, as where the dependences on the unphysical mass scales µ, λ, and M and the regulators ǫ IR , ǫ UV , v, δ i , and δ f are indicated in parentheses for each term. Everywhere in Eq. (55), α s (µ) is evaluated in the MS scheme using the two-loop formula (20), m is defined in the OS scheme, and the MEs are understood as their MS values O H [n] r (λ). This is necessary to ensure the exact cancellation of the singularities. The right-hand side of Eq. (55) is manifestly finite. The UV divergences cancel between δ vi and δ ct ; the IR singularities among δ vi , δ ct , δ op , δ fs , and dσ so ; and the Coulomb singularities between δ vi and δ op . Note that δ op is UV finite upon operator renormalization and that dσ is is IR finite upon photon PDF renormalization; therefore, δ op and dσ is in Eq. (55) do not depend any more on ǫ UV and ǫ IR , respectively. The full details about the cancellation of the IR singularities in the various cc Fock states n are presented in Table 4, where the following short-hand notation is used: Subprocess Source IR-singular term γγ → cc 3 S The right-hand side of Eq. (55) is independent of the cut-off parameters δ i and δ f . Specifically, the δ i dependence cancels between dσ is and dσ ha and the δ f dependence among δ fs , dσ so , and dσ ha . While δ fs , dσ is , and dσ so are known in analytic form, the phase-space integrals occurring in the evaluation of dσ ha are rather cumbersome and are thus solved numerically. Consequently, the cancellation of δ i and δ f has to be established numerically, too. The goodness of this cancellation is assessed in Section 3.
The right-hand side of Eq. (55) is independent of the renormalization scales µ and λ up to terms that are formally beyond NLO; this cancellation is not exact because the running of α s (µ) and O H [n] r (λ) is determined from the respective renormalization group equations, which resum logarithmic corrections to all orders. However, the righthand side of Eq. (55) depends on the factorization scale M already in NLO; this M dependence is cancelled up to terms beyond NLO once the cross section of single-resolved photoproduction, which depends on M through the photon PDFs f a/γ (x a , M), is added. If M is to be varied, this must thus be done simultaneously in the cross sections of direct and single-resolved photoproduction.

Numerical results
We are now in a position to present our numerical analysis of the inclusive production of prompt J/ψ mesons in two-photon collisions with direct photon interactions at NLO in the NRQCD factorization framework. We consider TESLA in its e + e − mode with √ s = 500 GeV, where the photons are produced via bremsstrahlung and beamstrahlung.
The discussion of the numerical analysis proceeds in three steps. We first specify our input parameters. We then verify that our numerical evaluation is independent of the technical cut-off parameters δ i and δ f and assess its dependence on the renormalization and factorization scales µ, λ, and M. Finally, we explore the phenomenological consequences of our analysis by studying the size and impact of the NLO corrections on the p T and y distributions of the cross section.

Input parameters
We use m = 1.5 GeV and α = 1/137.036. For direct photoproduction at NLO (LO), we employ the two-loop (one-loop) formula for α (n f ) s (µ) [33] with n f = 3 active quark flavours and Λ (3) QCD = 299 MeV (204 MeV) [49]. We use the photon PDFs from Glück, Reya, and Schienbein (GRS) [49], which are the only available ones that are implemented in the fixed-flavour-number scheme, with n f = 3. When we combine the contribution due to single-resolved photoproduction, which is so far only known at LO, with the NLO one due to direct photoproduction, then we nevertheless evaluate it using the NLO formula for α (n f ) s (µ) and the NLO set of the GRS photon PDFs so as to maximize the compensation of the M dependence. In order to render this choice fully consistent, we will have to include the NLO correction to the single-resolved contribution and also the double-resolved contribution at NLO, once they become available. On the other hand, we consistently evaluate the single-resolved contribution using the LO formula for α (n f ) s (µ) and the LO set of the GRS photon PDFs when we consider it separately or in combination with the direct contribution at LO. Our default choice of renormalization and factorization scales is µ = M = m T and λ = m. In want of NLO sets of J/ψ, χ cJ , and ψ ′ MEs, we adopt the LO sets determined in Ref. [25] using the LO set of proton PDFs from Martin, Roberts, Stirling, and Thorne (MRST98LO) [50]. Specifically, O ψ(nS) 3 S were extracted from the measured partial decay widths of ψ(nS) → l + l − and χ c2 → γγ [33], respectively, while O ψ(nS) 1 S were fitted to the transverse-momentum distributions of ψ(nS) and χ cJ inclusive hadroproduction [5] and the cross-section ratio σ χ c2 /σ χ c1 [51] measured at the Tevatron. The fit results for O ψ(nS) 1 S are strongly correlated, so that the linear combination with a suitable value of r, is quoted. Unfortunately, Eq. (55) is sensitive to linear combination of O ψ(nS) 1 S that is different from the one appearing in Eq. (57). In want of more specific information, we thus make the democratic choice O ψ(nS) 1 S We now discuss the photon flux functions that enter our predictions for photoproduction in the e + e − mode of TESLA. The energy spectrum of the bremsstrahlung photons is well described in the WWA by Eq. (27) of Ref. [52]. We assume that the scattered electrons and positrons will be antitagged, as was usually the case at LEP2, and take the maximum scattering angle to be θ max = 25 mrad [53]. The energy spectrum of the beamstrahlung photons is approximately described by Eq. (2.14) of Ref. [54]. It is controlled by the effective beamstrahlung parameter Υ, which is given by Eq. (2.10) of that reference. Inserting the relevant TESLA parameters for the √ S = 500 GeV baseline design specified in Table 1.3.1 of Ref. [55] in that formula, we obtain Υ = 0.053. We coherently superimpose the WWA and beamstrahlung spectra.
With the NLO corrections to single-and double-resolved photoproduction yet to be evaluated, we are not in a position to present a complete phenomenological prediction that could be confronted with experimental data as it stands. Therefore, we refrain from presenting a full-fledged quantitative estimate of the theoretical uncertainties. However, in the next section, we do investigate the dependences on the renormalization and factorization scales µ, λ, and M.

Academic study
For the purpose of the following technical study, it is sufficient to consider a typical kinematic situation. We thus choose as our reference quantity the differential cross section d 2 σ/dp T dy at p T = 5 GeV and y = 0.
In Figs. 6(a) and (b), the NLO result of direct photoproduction (solid lines), its hard, noncollinear component (dot-dashed lines), which corresponds to dσ ha in Eq. (55), and the remainder (dotted lines) are plotted as functions of δ i for δ f = 5 × 10 −3 and as functions of δ f for δ i = 5 × 10 −4 , respectively. For comparison, also the LO result (dashed lines), which is, of course, independent of δ i and δ f , is shown. Notice that δ i and δ f are varied over several orders of magnitude. They have to be chosen judiciously: if they are chosen too large, then the collinear approximation underlying δ fs and dσ is and the soft approximation underlying dσ so break down; if they are chosen too small, then the numerical evaluation of dσ ha becomes inaccurate. Our default values for the rest of the numerical analysis are δ i = 5 × 10 −4 and δ f = 5 × 10 −3 . Since the remainder depends logarithmically on δ i and δ f , it is represented by straight lines in Figs. 6(a) and (b), the abscissae of which are logarithmic. The same is true, to good approximation, for the hard, noncollinear component, and the combined NLO result is practically independent of δ i and δ f , as it should. This suggests that the numerical integration is rather precise and stable. We conclude that this error source can be safely neglected in the determination of the theoretical uncertainty of the NLO result.
In Fig. 7, the NLO (solid line) and LO (dashed line) contributions due to the cc Fock state n = 3 S are shown as functions of µ, while λ and M are kept fixed at their reference values. We focus on this particular channel because it is the only one already open at LO, so that a compensation of the µ dependence can occur. Notice that µ is varied by more than one order of magnitude, from m T /4 to 4m T . Passing from LO to NLO, the µ dependence is appreciably reduced, reflecting the partial compensation of the µ dependence of dσ 0 by the one of δ ct in Eq. (55). In fact, the related theoretical uncertainty amounts to +70 −30 % at LO and to ±18% at NLO. In Fig. 8, the NLO result is shown as a function of λ (solid line), while µ and M are kept fixed at their reference values. Notice that λ is varied by almost one order of magnitude, from m/2 to 4m. The related theoretical uncertainty amounts to ±50%. The λ-dependent terms in Eq. (55) are dσ 0 and δ op . In order to exhibit the partial compensation in λ dependence between these two terms, we also include in Fig. 8

Phenomenological study
In Fig. 10, we study d 2 σ/dp T dy (a) for y = 0 as a function of p T and (b) for p T = 5 GeV as a function of y. In each case, the LO (dashed line) and NLO (solid line) results of direct photoproduction as well as the LO result of single-resolved photoproduction evaluated with the LO versions of α (n f ) s (µ) and the photon PDFs (dotted line) are shown. Notice that our analysis is only valid for finite values of p T ; in the limit p T → 0, additional IR singularities occur, which require a more sophisticated scheme of phase space slicing. Therefore, we do not consider p T values below 2 GeV in Fig. 10(a). For a more detailed discussion of this point, we refer to Ref. [21]. From Fig. 10(a), we observe that, with increasing value of p T , the NLO result of direct photoproduction falls of considerably more slowly than the LO one. This is also evident from Fig. 11(a), where the NLO to LO ratio is shown for y = 0 as a function of p T (solid line). This feature may be understood by observing that so-called fragmentationprone partonic subprocesses [56] start to contribute to direct photoproduction at NLO, while they are absent at LO. Such subprocesses contain a gluon with small virtuality, q 2 = 4m 2 , that splits into a cc pair in the Fock state n = 3 S qq that are proportional to e 2 q ; the corresponding diagrams are the last two ones in Fig. 2 and the last six ones in Fig. 5, respectively. In single-resolved photoproduction, a fragmentation-prone partonic subprocess already contributes at LO; the relevant T -matrix elements is the one of γq → cc 3 S q that is proportional to e q . This explains why the solid and dotted curves in Fig. 10(a) run parallel in the upper p T range. At low values of p T , the fragmentation-prone partonic subprocesses do not matter, and the relative suppression of direct photoproduction is due to the fact that, at LO, this is a pure CO process.
At p T = 5 GeV, single-resolved photoproduction is still overwhelming, as is evident from Fig. 10(b). In Fig. 10(b), the two pronounced maxima in the dotted line may be traced to the T -matrix element of γq → cc[n]q that is proportional to e c , which contains a virtual gluon in the t channel that can become almost collinear with the incoming q quark, and to the one that is proportional to e q , which contains a virtual q quark in the u channel that can become almost collinear with the incoming photon.
Furthermore, we observe from Fig. 10(b) that the NLO correction to direct photoproduction dramatically increases towards the forward and backward directions. This feature is also nicely exhibited in Fig. 11(b), where the NLO to LO ratio is shown for p T = 5 GeV as a function of y (solid line). This is due to the finite remainders of the initial-state collinear singularities that were absorbed into the photon PDFs.
In Fig. 11, two kinds of QCD correction (K) factor are studied (a) for y = 0 as a function of p T and (b) for p T = 5 GeV as a function of y. The one that is defined as the ratio of the NLO and LO results of direct photoproduction (solid lines) was already discussed above. An alternative definition of K factor is obtained by adding the LO result of single-resolved photoproduction to both numerator and denominator, where it is evaluated with the NLO and LO versions of α (n f ) s (µ) and the photon PDFs, respectively (dashed lines). The rationale for this definition is the circumstance that the NLO result of direct photoproduction and the LO result of single-resolved photoproduction are interconnected by mass factorization as explained in Section 2.3.4. For the reasons given in Section 3.1, this definition of K factor is slightly inconsistent and the resulting values should, therefore, be taken with a grain of salt. Since the LO result of single-resolved photoproduction dominates both numerator and denominator of this K factor, the latter takes moderate values. The fact that these values are even below unity in the considered regions of phase space is due to the different choices of α (n f ) s (µ) and the photon PDFs for the single-resolved contribution in the numerator and denominator of this K factor.

Conclusions
In the case of charmonium, the experimental verification of the NRQCD factorization hypothesis in two-photon collisions at LEP2, photoproduction at HERA, and hadroproduction at the Tevatron has now practically come to a halt, one reason being that the available theoretical predictions for the production of prompt J/ψ mesons with finite values of p T in high-energetic photon-photon, photon-hadron, and hadron-hadron collisions are only of LO and thus suffer from considerable uncertainties, mostly from the dependences on the renormalization and factorization scales and from the lack of information on the nonperturbative MEs.
Having recently considered the cross sections of processes (2) [21] and (21)-(23) [19] in connection with the associated production of prompt J/ψ mesons with one or two identified energetic hadron jets, respectively, in this paper, we took the next step by studying at NLO the inclusive production of prompt J/ψ mesons with finite values of p T . This is the first time that an inclusive 2 → 2 process was treated at NLO in the NRQCD factorization framework. The technical difficulties that needed to be tackled include the treatment of UV, collinear, soft, and Coulomb singularities, the NRQCD operator renormalization, the mass factorization, and the analytic evaluation of five-point one-loop integrals with UV, soft, and Coulomb singularities.
As for the real corrections, we employed the phase-space slicing method to demarcate the regions of phase space containing soft and collinear singularities from the hard regions, where the phase-space integrations were carried out numerically. We verified that the combined result is, to very good approximation, independent of the choices of the cutoff parameters δ i and δ f , over an extended range of values. We worked in dimensional regularization in connection with the MS renormalization and factorization schemes, so that our NLO result depends on the QCD and NRQCD renormalization scales µ and λ, respectively, and on the factorization scale M connected with the collinear splitting of the incoming photons into massless qq pairs. While the µ and λ dependences are formally cancelled up to terms beyond NLO within direct photoproduction, the M dependence is only compensated by the LO cross section of single-resolved photoproduction. By the same token, the strong M dependence of the latter is considerably reduced by the inclusion of our new result. This is a crucial phenomenological merit of our work, as may be appreciated by observing that, at LO, the overwhelming bulk of the cross section of prompt J/ψ production in two-photon collisions is due to single-resolved photoproduction [20,21].
The K factor of direct photoproduction turned out to be very substantial at large values of p T because fragmentation-prone channels start to open up at NLO [see Fig. 11(a)]. In fact, at p T ≫ 2m, the p T distribution of direct photoproduction at NLO is rather similar to the LO result of single-resolved photoproduction both in shape and normalization [see Fig. 10(a)].
In order to complete the NLO treatment of prompt J/ψ production in two-photon collisions, we still need to evaluate the NLO corrections to single-and double-resolved photoproduction. Then, also prompt J/ψ production in photoproduction at HERA and hadroproduction at the Tevatron can be described at NLO. This will provide a solid basis for an ultimate test of the NRQCD factorization framework.