Inelastic heavy quarkonium photoproduction in p-p and Pb-Pb collisions at LHC energies

The study of heavy quarkonium has yielded valuable insight into the nature of the strong interaction, $ Q\bar{Q} $ bound states have provided useful laboratories for probing both perturbative and nonperturbative QCD. During the last years, the study of the heavy vector meson produced by photon-induced interactions at hadronic colliders has been strongly motivated by the possibility of constraining the dynamics of the strong interactions at high energies [1−5]. It also sheds light on the low-x physics and helps to constrain the nuclear parton distributions [6−9]. It is well known that this type of mechanism can be theoretically studied using the Weizsäcker-Williams approximation (WWA) [10−12]. The central idea of WWA is that the moving electromagnetic field of charged particles can be treated as a flux of photons. In an ultrarelativistic ion collider, these photons can interact with the target nucleus in the opposite beam (photoproduction) or with the photons of the opposite beam (two-photon reactions). At the CERN Large Hadron Collider (LHC) energies, the intense heavy-ion beams represent a prolific source of quasireal photons, hence enabling extensive studies of photon-induced physics. In the calculations, an important function is the photon flux function, which has different forms for different charged sources.

Although great development has been achieved, the features of WWA in inelastic heavy quarkonium photoproduction in heavy-ion collisions are rarely noticed. WWA is usually employed beyond its applicable scopes, and imprecise statements pertaining to the advantages of WWA were given [13−28]. For instance, the WWA is usually adopted in electroproduction reactions or exclusive processes, where the virtuality $ Q^2 $ of the photon is very small, controlled by $ m_e $ or the coherence condition. However, when the WWA is used in hadronic collisions, $ Q^{2} $ is controlled by the nucleus mass $ m_N $ and it is not obvious that the WWA is still valid. Particularly in the ultrarelativistic heavy-ion collisions at LHC energies, the influence of WWA becomes significant to the accuracy of describing photoproduction processes, since photon flux function scales as $ f_{\gamma}\propto Z^{2}\ln\sqrt{s}/m $, in which the collision energy $ \sqrt{s} $ and the squared nuclear charge $ Z^{2} $ turn into the very large enhancement factors to the cross sections. Thus, heavy-ion collisions have a considerable flux advantage over the proton. For these reasons, it is necessary to present a comprehensive analysis of WWA in inelastic heavy quarkonium photoproduction in heavy-ion collisions, and to estimate the important inaccuracies appeared in its application.

There are different channels which contribute to heavy quarkonium production. From the beam side, the different photon sources need to be considered [29]: coherent-photon emission (coh.), ordinary-incoherent photon emission (OIC), and ultra-incoherent photon emission (UIC). In the first type, the virtual photons are emitted coherently by the whole nucleus which remains intact after photons radiated. In the second and third types, the virtual photons are emitted incoherently by the protons and quarks inside nucleus, respectively, and nucleus will dissociate after photon radiation. To avoid confusion, the terminology "elastic" and "inelastic" describe the case of whether the target nucleus remains intact or is allowed to break up after scattering with photons. When these different photon sources are considered together, we have to weight its relative contributions to avoid double counting. Meanwhile, in the final state, there are two types of inelastic productions need to be distinguished: direct and fragmentation contributions [30−32]. The fragmentation process is described by the fragmentation functions to specify the probability of final partons (gluons or quarks) hadronizing into quarkonia bound states. The fragmentation contribution originates from the large $ p_{T} $ region, where one encounters large logarithms of $ p_{T}^{2}/m_{Q}^{2} $, such large logs are resummed into the fragmentation functions. This essentially means that fragmentation mechanism can be only used in the large $ p_{T} $ region, and one can not naively add the direct and fragmentation contributions together, since this will also cause double counting [32−34]. But in fact, the double counting exists in most works, and the fragmentation formalism is employed beyond its validity range [14−22].

There are numerous studies on these processes, however, the application of UIC, to our knowledge, is insufficient in inelastic heavy quarkonium photoproduction. For instance, Gonçalves et al. systematically studied the exclusive production of vector mesons in hadronic collisions considering different phenomenological models in Refs. [36, 37]. Machado et al. studied the inelastic and exclusive heavy quarkonium photoproductions within the color dipole formalism [38, 39]. In Refs. [40−42], Klein and Nystrand studied the exclusive vector meson production via photon-Pomeron or photon-meson interactions, and discussed the interlay between photoproduction and two-photon interaction [43]. Ducati et al. investigated the exclusive $ J/\psi $ photoproduction and the radially excited $ \psi(2S) $ state off nucleons in p-p collisions according to the light-cone dipole formalism [44]. There are also many other relevant works, however the photon emission types in all of these works are coherent, with the incoherent-photon emission being neglected. Furthermore, the UIC photoproduction, which is best treated as inclusive processes, can provide additional corrections to the central collisions. For instance, the authors in Refs. [18−20] have investigated the inelastic dileptons, photons, and light vector mesons productions at the LHC energies. These works show that the UIC photoproduction enhances the contribution of massless and light final-state particles in the central collisions. However, the correction is not clear for heavy quarkonium due to its large mass.

According to the above purposes, in the present work, we investigate the inelastic photoproductions of charmonium ($ J/\psi $, $ \psi(2S) $) and bottomonium ($ \Upsilon(nS) $) in p-p and $ Pb $-$ Pb $ collisions at LHC energies. An exact treatment is performed that recovers the WWA near the $ Q^{2}\sim0 $ domain and can be considered as the generalization of leptoproduction [45]. The full kinematical relations matched with the exact treatment are also obtained. We present a consistent analysis of the features of WWA in heavy-ion collisions comparing with the exact results, and study the double counting. In addition, we estimate the contribution of ultra-incoherent channel to the inelastic heavy quarkonia photoproductions.

It is necessary to mention that there is a important work recently done by Wangmei Zha et al., where they address similar challenges in reconciling WWA validity [46]. Based on the QED approach, Wangmei Zha et al. developed a spatially-dependent photon flux distribution and established a more precise relationship between the photon transverse momentum distribution and impact parameter of collisions. They justified the inadequacy of the WWA model in describing the photon flux for electron-ion collisions, and pointed out that the QED approach provides a more realistic basis for calculating the impact parameter dependence of photoproduction processes in electron-ion collisions. Within this refined method, they explores the potential of utilizing nuetron tagging from Coulomb excitation of nuclei to effectively determine centrality for exclusive photoproduction in electron-ion collisions. Their study offers a new methodology for exploring the spatial and momentum structure of gluons in nuclei, and offering novel insights for experimental design and data analysis.

We organize the paper as follows. In Sec. II, we presents the formalism of exact treatment for the inelastic heavy quarkonium photoproduction in heavy-ion collisions. Based on the Martin-Ryskin method, we consider the coh., OIC, and UIC processes simultaneously. And according to the BCCKL method, we match the fixed order and fragmentation contributions. In Sec. III, we turn the accurate formula into the WWA one near the region $ Q^{2}\sim0 $, and study the several widely utilized photon fluxes. In Sec. IV, we numerically calculate the $ Q^{2} $, y, z, $ \sqrt{s} $, and $ p_{T} $ dependent differential cross sections, and the total cross sections at LHC energies. Finally, in Sec. V we summarize our paper.

For heavy quarkonium production and decay, an effective field theory known as nonrelativistic QCD (NRQCD), has been proposed to explain the huge discrepancy between the theoretical predictions and experimental measurements of the transverse momentum distribution of $ J/\psi $ production at the Tevatron. This scheme has proven to be highly successful in numerous applications [47]. In this section, we employ this scheme to describe the inelastic heavy quarkonium photoproduction. The NRQCD scheme is based upon a double expansion in $ \alpha_{s} $ and ν (the heavy-quark relative velocity in quarkonium rest frame), and its form for inelastic quarkonium H production is

Here, $ d\sigma_{A+B\rightarrow Q\bar{Q}_{[1,8]}[n]X} $ are the process-dependent short-distance coefficients (SDCs), which can be computed in perturbative QCD by expansion in $ \alpha_{s} $ and correspond to the production of a heavy quark-antiquark pair $ Q\bar{Q}_{[1,8]}[n] $ in a specific color and angular momentum state $ n=^{2S+1}L_{J}^{(c)} $. The term $ \langle{\cal{O}}^{H}_{[1,8]}[n]\rangle $ denotes the NRQCD long-distance matrix elements (LDMEs), which describe the probability of hadronization and have a well-defined scaling with v ($ v\simeq0.08 $ and $ 0.25 $ for bottomonia and charmonia, respectively). Since the color-singlet LDME for quarkonium production $ \langle{\cal{O}}^{\psi}[^{3}S_{1}^{(1)}]\rangle $ is related to the color-singlet LDME for quarkonium decay, it can be determined in lattice QCD, from potential models, or from the ψ decay rates into lepton pairs [33]. On the other hand, it is not known how to compute the color-octet production LDMEs from first principles, and they are usually determined by fitting some experimental data. Notably, if NRQCD factorization holds, the LDMEs should be universal.

For $ 1^{--} $ quarkonia productions, $ J/\psi $, $ \psi\left(2S\right) $ and $ \Upsilon\left(nS\right) $, the sum over n, truncated at order $ v^{4} $, involves four LDMEs [47]: $ \langle{\cal{O}}^{H}[^{3}S_{1}^{(1)}]\rangle $, $ \langle{\cal{O}}^{H}[^{3}S_{1}^{(8)}]\rangle $, $ \langle{\cal{O}}^{H}[^{1}S_{0}^{(8)}]\rangle $, and $ \langle{\cal{O}}^{H}[^{3}P_{J}^{(8)}]\rangle $. Due to their lengthy expressions, we summarized them in Appendix A for completeness and convenience.

The exact treatment for the inelastic heavy quarkonium photoproduction in heavy-ion collisions can be considered as proceeding in two steps. In the first step, the density matrix of virtual photon should be expanded using the polarization operators, according to the fact that the photon radiated from the projectile is off mass shell and no longer transversely polarized. In the second step, the square of the electric form factor $ D\left(Q^{2}\right) $ is adopted as the weighting factor (WF) for different charged sources to avoid double counting.

By comprehensively analyzing the terms neglected in transiting from the exact formula of Fig. 1(a) to the WWA one, we can naturally estimate the features of WWA in heavy-ion collisions. In our case, the details of the result depend essentially on the characteristic behaviour met when moving off mass shell for photoprocess amplitudes. In the first step of exact treatment, we should derive the general form of cross section for the inelastic heavy quarkonium photoproduction in Fig. 1(a):

Figure 1.  (a): The general inelastic heavy quarkonium photoproduction. The virtual photon emitted from α interacts with parton b of nucleus B, α can be the nucleus or its charged parton (protons or quarks). (b): real photoabsorption.

where $ x_{b}=p_{b}/p_{B} $ is the momentum fraction of the massless parton b struck by the virtual photon, and the distribution function of parton b in nucleus B is

where $ R_{B}\left(x,\mu^{2}\right) $ is the nuclear modification factor [48], Z is the proton number, N is the neutron number. $ p\left(x,\mu^{2}\right) $ and $ n\left(x,\mu^{2}\right) $ are the parton distributions of the proton and neutron, respectively. According to the NRQCD scheme, the partonic cross section in Eq. (2) can be written as

Denoting the virtual photo-absorption amplitude by $ M^{\mu} $, we obtain the SDCs in the parton level

where $ e_{\alpha} $ is the charge of α, $ \alpha_{\text{em}} $ is the electromagnetic coupling constant, $ E_{\alpha'} $ is the energy of $ \alpha' $, and Π is a phase space volume of the produced particle system (with total momentum k). Keeping in mind the process in which photons may be emitted by various particles, we present a generalized density matrix of the virtual photon as:

where $ C\left(Q^{2}\right) $ and $ D\left(Q^{2}\right) $ are the general notations of form factors for α. Note that the $ \rho^{\mu\nu} $ is non-diagonal, indicating that the virtual photons are polarized. The expression in Eq. (5) is formulated to naturally introduce the terminology suitable for WWA. Namely, instead of discussing nucleus-nucleus collisions [Fig. 1(a)], one can refer to the collisions of a virtual photon off the nucleus [Fig. 1(b)].

Now we employ the accurate expression Eq. (5) to give the $ Q^{2} $- and y-dependent differential cross section for the inelastic heavy quarkonium photoproduction. It is more convenient to perform the calculations in the rest frame of α, where $ |{\bf{q}}|=|{\bf{p}}_{\alpha'}|=r $, $ Q^{2}=-q^{2}=\left(p_{\alpha}-p_{\alpha'}\right)^{2}= 2m_{\alpha}\left(\sqrt{r^{2}+m_{\alpha}^{2}}-m_{\alpha}\right) $, $ d^{3}p_{\alpha}'=r^{2}drd\cos\theta d\varphi $, and $ y=\left(q\cdot p_{b}\right)/\left(p_{\alpha}\cdot p_{b}\right)=\left(q_{0}-|{p}_{b}|r\cos\theta/E_{b}\right)/m_{\alpha} $ (which measures the relative energy loss of α in the lab-system). By doing the following transformation

the differential cross section of Eq. (5) can be turned into (the details of $ {\cal{J}} $ are given in Appendix B)

and

where $ s_{\alpha b}=(p_{\alpha}+p_{b})^{2} $ and $ \hat{s}=(q+p_{b})^{2} $ are the energy square in the α-b and $ \gamma^{*} $-b partonic processes, respectively. $ p_{\text{CM}} $ and $ \hat{p}_{\text{CM}} $ are the corresponding momenta. The details are summarized in Appendix B.

After integrating over the phase space volume Π, the following quantity will be included in the result Eq. (8):

where $ W^{\mu\nu} $ is the absorptive part of the $ \gamma b $ amplitude [Fig. 1(b)], connected with the cross section in the usual way. The tensors according to which $ W^{\mu\nu} $ is expanded, can be constructed only from the q, $ p_{b} $, and $ g^{\mu\nu} $ tensor. Considering the gauge invariance $ q^{\mu}W^{\mu\nu}=q^{\nu}W^{\mu\nu}=0 $, it is convenient to use the following transverse and longitudinal polarization operators [49]

which satisfy the relations: $ q_{\mu}\epsilon_{T}^{\mu\nu}=q_{\mu}\epsilon_{L}^{\mu\nu}=0 $, $ \epsilon_{T \mu}^{\mu}=-2 $, and $ \epsilon_{L \mu}^{\mu}=-1 $. Furthermore,

is the polarization tensor of an unpolarized spin-one boson with mass $ q^{2} $. Having expanded $ W^{\mu\nu} $ in these tensors, we get

These Lorentz scalar functions $ W_{\text{T}} $ and $ W_{\text{L}} $ are connected with the transverse and longitudinal photon absorption cross sections $ \sigma_{\text{T}} $ and $ \sigma_{\text{L}} $, respectively:

Substituting Eqs. (13), (14) into Eq. (8), we finally obtain

where the relations: $ d\sigma_{T}/d\hat{t}=F_{b}[n]T_{b}[n] $ and $ d\sigma_{L}/d\hat{t}=-2Q^{2}F_{b}[n]L_{b}[n] $, are employed. $ F_{b}[n] $, $ T_{b}[n] $ and $ L_{b}[n] $ are the functions of Mandelstam variables $ \hat{s} $, $ \hat{t} $, $ \hat{u} $, and $ Q^{2} $, which can be found in Ref. [26]. The coefficients $ \rho^{ab} $ are the elements of the density matrix Eq. (6) in the $ \gamma b $-helicity basis:

Here we come to the position to derive the second step of exact treatment, the details of the form factors in Eq. (16) need to be distinguished in each photon emission channel. In the Martin-Ryskin method [50], the probability or weighting factor (WF) of the coherent-photon emission is given by the square of the electric form factor in p-p collisions: $ w_{\text{coh}}=G_{\text{E}}^{2}\left(Q^{2}\right)= 1/ \left(1+Q^{2}/0.71\; \text{GeV}\right)^{4} $, where the effect of the magnetic form factor is neglected. We adopt this central idea to deal with the situation in heavy-ion collisions, where the magnetic form factor is also included. In the case of coherent-photon emission, the photon emitter α is nucleus, and thus the general notations $ C(Q^{2}) $ and $ D(Q^{2}) $ in Eq. (16) are the elastic nucleus form factors. In p-p collisions, α presents proton, $ C(Q^{2}) $ and $ D\left(Q^{2}\right) $ turn into the Sachs combinations [22]

where $ \mu_{p}=2.79 $ is the magnetic dipole moment. In $ Pb $-$ Pb $ collisions, α is lead ion, $ C(Q^{2}) $ and $ D(Q^{2}) $ are

where

is the electromagnetic form factor parameterization from the STARlight MC generator [51], in which $ R_{A}=1.1A^{1/3}\; \text{fm} $, $ a=0.7\; \text{fm} $ and $ Q=\sqrt{Q^{2}} $.

In the case of ordinary- and ultra-incoherent photon emissions, the contributions must be multiplied by the 'remaining' probability, $ 1-w_{\text{coh}} $, to avoid double counting [50]. For ordinary-incoherent photon emission in $ Pb $-$ Pb $ collisions, α is the protons within the lead ion, and the corresponding $ D(Q^{2}) $ and $ C(Q^{2}) $ are

For ultra-incoherent photon emission, α is the individual quarks within the nucleus. The corresponding $ C(Q^{2}) $ and $ D(Q^{2}) $ in p-p collisions have the following form

In $ Pb $-$ Pb $ collisions, since the neutron can not emit photon coherently, the WF for proton and neutron inside lead ion are different

Now we switch the general expression [Eq. (15)] to each specific channel in inelastic photoproduction processes in heavy-ion collisions. In the initial state, the processes may be direct or resolved that are sensitive to the gluon distribution in the nucleus [4]. The photons emitted from the projectile, 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 $ \gamma^{*}b\rightarrow H+X $ receives contributions from both direct and resolved channels. All two contributions are formally of the same order in the perturbative expansion and must be included [52]. Actually, as always with photons, the situation is quite complex. Together with the three different photon emissions mentioned in Section I, the complete description of the heavy quarkonium production requires the calculation of six classes of processes [Figs. 2-3]: coherent-direct (coh.dir.), coherent-resolved (coh.res.), ordinary-incoherent direct (OIC dir.), ordinary-incoherent resolved (OIC res.), ultra-incoherent direct (UIC dir.), and ultra-incoherent resolved (UIC res.) processes. These abbreviations will appear in many places of the remaining content.

Figure 2.  (a): coherent-direct process in which the virtual photon emitted from the whole incident nucleus A interacts with parton b of target nucleus B via the $ \gamma^{*} $-q Compton scattering and $ \gamma^{*} $-g fusion, and A remains intact after photon emitted. (b): incoherent-direct process in which the virtual photon emitted from the quark a within nucleus A interacts with parton b, and A is allowed to break up after photon emitted. Where photon emitter a are proton and quark for ordinary-incoherent photon emission (OIC) and ultra-incoherent photon emission (UIC), respectively.

Figure 3.  (a): coherent-resolved process in which the parton $ a' $ of hadron-like photon emitted from nucleus A, interacts with the parton b of target B via q-g Compton scattering, q-$ \bar{q} $ annihilation and g-g fusion. (b): incoherent-resolved process in which the parton a inside nucleus A emit a resolved virtual photon, then the parton $ a^\prime $ of this resolved photon interacts with parton b inside target B like a hadron, and A is break up after photon emitted. Where photon emitter a are proton and quark for ordinary-incoherent photon emission (OIC) and ultra-incoherent photon emission (UIC), respectively.

In the case of direct photoproduction [Fig. 2], the corresponding differential cross sections are

where the factor of two arises because both nuclei emit photons and thus serve as targets. The partonic cross section can be derived from Eq. (15) with $ m_{\alpha}=m_{q}=0 $ and $ e_{\alpha}=e_{a} $, where $ e_{a} $ is the charge of massless quark a.

In the resolved photoproduction [Fig. 3], the corresponding differential cross sections are

where $ z_{a}'=p_{a'}/q $ and $ f_{\gamma}\left(z_{a'},\mu^{2}\right) $ denote the parton's momentum fraction and the parton distribution function of the resolved photon [53], respectively. The involved partonic cross sections can be found in Ref. [52].

The distributions in $ p_{T} $ and inelastic variable $ z=(p_{H}\cdot p_{b})/(q\cdot p_{b}) $ can be obtained using the Jacobian transformation. In the final state, we need to classify the two types of inelastic photoproductions. In the first type, direct heavy quarkonium produced from the γ-g fusion, annihilation and Compton scattering of partons. In the second type, fragmentation heavy quarkonium produced through the final fragmentation of a parton. In the following, we will take into account all of these aspects.

Before doing the transformation the Mandelstam variables in $ \gamma^{*} $-b parton level should be expressed as:

where $ y_{r} $ is the rapidity, $ M_{T}=\sqrt{M_{H}^{2}+p_{T}^{2}} $ is the transverse mass of heavy quarkonium, $ \hat{E}_{\gamma} $, $ \hat{E}_{b} $, and $ \hat{p}_{\text{CM}} $ are the corresponding energies and momentum. The details are summarized in Appendix B.

In the case of direct-photon processes, the variables $ x_{b} $ and $ \hat{t} $ should be transformed as

and the corresponding differential cross sections are

In the case of resolved processes, we should choose the variables $ \hat{t}^{*} $ and $ z_{a'} $ to do the similar transformation

the corresponding differential cross sections are

where the Mandelstam variables of resolved photoproductions are the same as Eq. (29) but for $ Q^{2}=0 $.

For the z distribution, we should rewrite the Mandelstam variables in Eq. (29) as follows:

thus, we can do the similar Jacobian transformation. The relevant cross sections in z distribution can be obtained in the following manner

The Jacobian factors $ {\cal{J}} $ and corresponding kinematical boundaries are summarized in Appendix B.

The fragmentation heavy quarkonium production is also an important channel which is described by the fragmentation functions (FFs). This factorization theory of quarkonium production that has been proven at present is the collinear factorization method [54], where one can compute rates at leading power (LP) or next to leading power (NLP) in $ m_{Q}^{2}/p_{T}^{2} $. LP contributions can be factorized into partonic cross sections to produce a specific single parton convolved with one particle FFs [30]. NLP contributions can be factorized to produce two specific partons convolved with two-parton FFs [54]. The Feynman diagrams in the cut diagram notation for these two contributions are shown in Fig. 4. This fragmentation picture dominates over all other creation mechanisms at large $ p_{T} $ range. In particular, gluon fragmentation represents the dominant source of high energy prompt quarkonia at hadron colliders. Of course, because the LP and NLP contributions represent the leading and first subleading terms in an expansion in powers of $ m_{Q}^{2}/p_{T}^{2} $, one would not expect them to be valid unless $ p_{T} $ is significantly greater than $ m_{Q} $. Fragmentation predictions for charmonium differential cross sections are therefore unreliable at low $ p_{T} $ domain. In order to suppress NLP contributions and to ensure the validity of fragmentation mechanism, Ref. [35] suggested that a criterion $ p_{T}>3M_{J/\psi} $ be used in comparing data with theory. In Refs. [32−34], a criterion $ p_{T}>10\; \text{GeV} $ was used in their predictions. Similarly, fragmentation results for bottomonium production are untrustworthy throughout the $ p_{T}<15\; \text{GeV} $ region [31].

Figure 4.  pQCD factorization diagrams of heavy quarkonium production. Left: sing-parton (here taking gluon as an example) fragmentation. Right: heavy quark-pair fragmentation.

In present paper, we adopt the LP-factorization formalism to calculate the fragmentation contributions to heavy quarkonium production. Calculations of these fragmentation contributions, at any given order in $ \alpha_{s} $, are much simpler than those of full fixed-order calculations. According to the NRQCD scheme, the inelastic quarkonium H production in two-body collisions can be written as Eq. (1), where the SDCs, $ d\sigma_{A+B\rightarrow Q\bar{Q}_{[1,8]}[n]X} $, describes the production of a $ Q\bar{Q} $ pair in color and angular momentum state n. At large transverse momentum $ p_{T} $, one can apply LP factorization to the SDCs in Eq. (1). The general expression is [30, 55]

where $ d\hat{\sigma}_{A+B\rightarrow c+X} $ are the parton production cross sections (PPCSs) to produce a parton c, and $ D_{c\rightarrow Q\bar{Q}[n]} $ are FFs for a parton c to fragment into a $ Q\bar{Q} $ pair with quantum numbers n. In the LP factorization, the SDCs of Eq. (40) for photoproduction and initial partons hard scattering (had.scat.) can be summarized as following master formula

where $ z_{c} $ is the longitudinal momentum fraction of the hadron relative to the fragmenting parton. $ f_{a/A}\left(x_{a},\mu^{2}\right) $ is the PDF of parton $ a=g,q,\bar{q} $ in nucleus A, or the flux function of photon $ a=\gamma $. $ f_{k/a}\left(z_{a'},\mu^{2}\right) $ is $ \delta_{ak}\delta\left(1-z_{a'}\right) $ in the "direct" case (since the parton a is the photon itself), or the PDF of parton k in the resolved photon a. The parton distributions in the photon behave like $ \alpha/\alpha_{s}(Q^2) $ for large $ Q^2 $. Therefore the additional power of $ \alpha_{s} $ contained in the "resolved" component as compared to the "direct" one is compensated. The hard PPCSs, $ d\sigma_{kb\rightarrow cX} $, and FFs have been calculated to NLO accuracy, they are expansions in powers of $ \alpha_{s} $

In this paper, we use the LP factorization approximation for the SDCs to compute fragmentation contributions that augment the LO fixed order calculations of the SDCs which are presented in the previous sections. Since direct and fragmentation contributions can not be directly added up, we adopt the method developed by Bodwin, Chao, Chung, Kim, and Lee (BCCKL) to avoid double counting, where the matching rule between fixd-order and fragmentation has been performed [32−34]. Therefore, we combine the LO fixed-order calculations with LP fragmentation corrections [Eq. (40)] according to the formula

where we denotes the LO ($ \alpha_{s}^{3} $) contributions to SDCs by $ d\sigma_{ \text{LO}} $, which are discussed in previous sections and start from Eq. (5). We denotes the LP fragmentation corrections by $ d\sigma^{ \text{LP-frag.}} $, and only consider the contributions at order $ \alpha_{s}^{4} $, because the authors in Ref. [33] have shown that the LP fragmentation contributions at order $ \alpha_{s}^{5} $ are small, and only have little effect on the cross section prediction for photoproduction.

For $ d\sigma^{ \text{LP-frag.}} $, we combine the PPCSs and FFs both at order $ \alpha_{s}^{2} $. The LP fragmentation process also include the six classes of sub-processes which are described in Figs. 2, 3. According to Eq. (40), in calculations we just need to replace the LO SDCs [Eqs. (31)-(37)] with the convolution of PPCSs and FFs. The PPCSs of the direct and resolved processes can be found in Refs. [56, 57]. We consider gluon and light-quark fragmentations. The gluon FFs $ D_{g\rightarrow Q\bar{Q}[n]} $ are given for the $ ^{1}S_{0}^{(8)} $ at order $ \alpha_{s}^{2}\; ( \text{LO}) $ in Refs. [58, 59], for the $ ^{3}S_{1}^{(8)} $ at orders $ \alpha_{s}\; ( \text{LO}) $ and $ \alpha_{s}^{2}\; ( \text{NLO}) $ in Refs. [60, 61], and for the $ ^{3}P_{J}^{(8)} $ at order $ \alpha_{s}^{2}\; ( \text{LO}) $ in Refs. [59, 60]. Because the gluon FF for the $ ^{3}S_{1}^{(1)} $ begins at order $ \alpha_{s}^{3}\; ( \text{LO}) $, we do not consider it here. The light quark FF $ D_{q\rightarrow Q\bar{Q}}[n] $ is given for the $ ^{3}S_{1}^{(8)} $ at order $ \alpha_{s}^{2}\; ( \text{LO}) $ in Refs. [61, 62], the other light quark FFs vanishes though order $ \alpha_{s}^{2} $.

For achieving $ p_{T} $ distribution, we should rewrite the Mandelstam variables as follows,

and since the parton c is taken to be lightlike by neglecting the parton mass, $ z_{c}=p'_{c}/p_{c}=2p_{T}\cosh y_{r}/\sqrt{\hat{s}} $. Then the variables $ z_{c} $ and $ \hat{t} $ can do the transformation

The connection between the process in Fig. 1(a) and (b) is evident. By Fourier-transforming the electric and magnetic fields of an ultrarelativistic charged point-like particle, the photoproduction process can be expressed in terms of the real photo-absorption cross section with the photon spectrum. This idea was originally pointed out by Fermi [10], and was independently developed for the process involving relativistic collisions of charged particles by Weizsäcker and Williams, and the method is now known as the Weizsäcker-Williams approximation (WWA) [11]. An essential advantage of WWA consists in the fact that, when using it, it is sufficient to obtain the photo-absorption cross section on the mass shell only. Details of its off mass-shell behavior are not essential. In the present section, we switch the accurate expression Eq. (15) into the WWA form by taking $ Q^{2}\rightarrow0 $, and discuss a number of widely employed photon spectra. There are two simplifications: the scalar photon contribution $ \sigma_{L} $ is neglected; the term of $ \sigma_{T} $ is substituted by its on-shell value. This provides us a powerful approach to study the features of WWA in inelastic heavy quarkonium photoproduction in heavy-ion collisions.

Taking $ Q^{2}\rightarrow0 $, Eq. (15) turns into:

where the contribution of $ \sigma_{L} $ and the terms proportional to $ Q^{2} $ are neglected in the limit $ Q^{2}\rightarrow 0 $. And the general form of the photon spectrum $ f_{\gamma}(y) $, which is associated with various particles, reads

In the case of coherent-photon emission, Ref. [5] presented a modified photon flux function of proton from Eq. (47). Neglecting the magnetic form factor and adopting the dipole form of electric form factor of proton: $ C(Q^{2})=D(Q^{2})=G_{\text{E}}^{2}(Q^{2}) $, and employing the coherent condition $ Q^{2}\leq 1/R^{2}_{A} $ ($ Q^{2}_{\text{max}}=0.027, y_{\text{max}}=0.16 $), one obtains with $ a=2m_{p}^{2}/Q^{2}_{\text{max}} $ and $ b=2m_{p}^{2}/0.71=2.48 $,

where x depends on y,

.

Actually, the origin of various widely employed photon spectra is another plane wave form, which is given in Ref. [49] and can be written as

this form is achieved from the complete form Eq. (47) by assuming that, $ Q^{2}_{\text{min}}=y^{2}m_{\alpha}^{2}/(1-y) $, which is the LO term of the following complete expression

This approximation is only available when $ m_{\alpha}^{2}\ll1\; \text{GeV}^2 $, however $ m_{p}^{2} $ and $ m_{Pb}^{2} $ do not satisfy the condition; this is a error source in various spectra. Especially for lead ion, this approximation cause erroneous results.

Drees and Zeppenfeld (DZ) provided another widely used photon distribution function of proton [15−21], which is the approximate analytic form of Eq. (50). Assuming: $ Q^{2}_{\text{max}}\rightarrow \infty $, $ C\left(Q^{2}\right)=D\left(Q^{2}\right)=G_{\text{E}}^{2}(Q^{2}) $, and $ Q^{2}-Q^{2}_{\text{min}}\approx Q^{2} $, they obtained

where $ A=\left(1+0.71\; \text{GeV}^{2}/Q^{2}_{\text{min}}\right) $. $ f^{\gamma}_{\text{DZ}} $ properly include electric form factor of proton to describe the situation of the proton as photon emitter. Because WWA is usually adopted in electroproduction processes, if one directly obtains the spectrum of proton from that of electron by just replacing the $ m_{e} $ with $ m_{p} $, it would extensively overestimate the cross section. In Ref. [23], Drees, Ellis, and Zeppenfeld (DEZ) also performed a spectrum of lead ion. Assuming $ y\ll1 $, $ Q^{2}_{\text{max}}\sim\infty $, $ C_{Pb}\left(Q^{2}\right)=0 $, and $ D_{Pb}\left(Q^{2}\right)\approx\exp\left(-\dfrac{Q^{2}}{Q^{2}_{0}}\right) $, they achieved

where $ Q^{2}_{\text{min}}=m^{2}_{Pb}y^{2} $, $ \Gamma(a,Q^{2}_{\text{min}}/Q^{2}_{0})=\int_{y}^{\infty}t^{a-1}e^{-t}dt $ is the incomplete Gamma Function. It should be noticed that, $ y\ll1 $ means $ Q^{2}_{\text{max}}\sim0 $, which contradicts with the assumption $ Q^{2}_{\text{max}}\sim\infty $.

Based on Eq. (52), Nystrand derived a modified photon spectrum of proton which includes the $ Q^{2}_{\text{min}} $ term in Eq. (50) and can be presented as [25]

In addition, the effect of including the magnetic form factor of the proton has been estimated by Kniehl, the final expression $ f^{\text{Kn}}(y) $ (Eq. (3.11) of Ref. [22]) is too long to include here, but will be discussed further below.

Another most important approach for coherent photon spectrum is the semiclassical impact parameter description, which excludes the hadronic interaction easily. The calculation of the semiclassical photon spectrum for the case of E1 (electric dipole) excitation is explained in Ref [63], and the result is

where υ is the velocity of the point charge $ Ze $, $ K_{0}(x) $ and $ K_{1}(x) $ are the modified Bessel functions, and $ \xi=\omega b_{\text{min}}/\gamma_{L}v=b_{\text{min}}m_{A}y/v $. Although the semi-classical photon spectrum is $ Q^{2} $-independent, the boundary of y should be constricted by coherence condition [3].

In the case of ultra-incoherent photon emission, Ref. [64] presented a photon spectrum inside a quark from Eq. (50). Neglecting the weighting factor in Eqs. (21) and (22), and setting $ Q_{\text{min}}^{2}=1\; \text{GeV}^{2} $ and $ Q^{2}_{\text{max}}=\hat{s}/4 $, one obtains

Finally, Brodsky, Kinoshita and Terazawa calculated another important incoherent photon spectrum in Ref. [65], which is originally derived for e-p scattering,

In high-energy physics, the WWA is often used in studying hadronic interactions and heavy-ion collisions: Wangmei Zha et al. generalized WWA to calculating the cross section and $ p_{T} $ distribution of the Breit-Wheeler process in relativistic heavy-ion collisions and their dependence on collision impact parameter $ (b) $ [66]. Based on the WWA, theses author also analyzed gluon shadowing in heavy nuclei through Bayesian Reweighting of coherent $ J/\psi $ photoproduction in UPCs [67]; ATLAS collaboration studied the exclusive muon pair production in the two-photon process [68]; d'Enterria and Silveira [69] proposed using equivalent photons from colliding lead ions to observe light-by-light scattering experimentally, later, the direct observations were obtained by the ATLAS [70] and CMS [71] collaborations using lead-lead collisions; Monte-Carlo event generators, such as Pythia [72] and STARlight [41], also utilize the WWA. Generally, particle production in hadronic collisions is often modelled using WWA [41, 43]. Although great development has been achieved, the properties of WWA in inelastic heavy quarkonium photoproduction in heavy-ion collisions are often neglected, and the imprecise statements were given [13−28]. Especially in the ultrarelativistic heavy-ion collisions at LHC energies, WWA becomes deterministic to the accuracy of describing photoproductions, since $ f_{\gamma}\propto Z^{2}\ln\sqrt{s}/m $, $ \sqrt{s} $ and $ Z^{2} $ are the very large factors. Therefore, we will comprehensively analyze the WWA in the mentioned case, and estimate the inaccuracies in above spectra.

Now, we are ready to show the numerical results. In the calculations, we adopt an approximation $ m_{q}=M_{H}/2 $ for the quark mass, where $ M_{H} $ is the mass of heavy quarkonium. All the masses are taken from PDG [73]. We use $ m_{p}=0.938\; \text{GeV} $, $ \alpha_{\text{em}}=1/137.036 $, $ p_{T\text{min}}=M_{H} $ and the two-loop QCD coupling constants $ \alpha_{s} $ with $ n_{f}=3 $ and $ \Lambda=0.2\; \text{GeV} $ [74]. We adopt MSHT20 LO (NLO) set for PDFs with $ n_{f}=3 $ [75], and the factorization scale $ \mu=\sqrt{M_{H}^{2}+Q^{2}} $. The LDMEs and the full kinematical relations are summarized in Appendices.

First of all, we choose the $ J/\psi $ as an example, to comprehensively study the properties of WWA in inelastic heavy quarkonium photoproduction in heavy-ion collisions. Therefore, the distributions in $ Q^2 $, y, z, and $ \sqrt{s} $ are plotted in Figs. 5-8, where all the results are the sum of direct and resolved contributions, and do not include the fragmentation contributions. The left panels of Figs. 5-8 show the relative errors with respect to the exact results, the central and right panels show the exact results of cross sections in p-p and $ Pb $-$ Pb $ collisions, respectively.

Figure 5.  (color online) $ Q^{2} $ distribution of the $ J/\psi $ photoproduction at LHC energies. (a): relative error of the WWA result with respect to the exact one. (b), (c): exact results of the $ Q^{2} $-dependent differential cross sections in p-p and $ Pb $-$ Pb $ collisions, respectively. Black solid line---coherent-photon emission [coh.(dir.+res.)]. Red dashed line---ultra-incoherent photon emission [UIC (dir.+res.)]. Blue dotted line---ordinary-incoherent photon emission [OIC (dir.+res.)].

Figure 6.  (color online) Same as Fig. 5 but for y distribution.

Figure 7.  (color online) (a): relative error of the WWA result with respect to the exact one. (b), (c): exact results of $ d\sigma/dz $ in p-p and $ Pb $-$ Pb $ collisions. Black solid and dark cyan dot-dashed lines denote the coherent-direct and resolved photon emissions, respectively. Red dashed and magenta dot-dot dashed lines denote the ultra-incoherent direct and resolved photon emissions, respectively. Blue dotted and dark yellow short dashed lines denote the ordinary-incoherent direct and resolved photon emissions, respectively.

Figure 8.  (color online) Same as Fig. 5 but for the total cross sections as a function of $ \sqrt{s} $.

In Fig. 5, there is only one curve in panel $ (a) $, since the curves of coh., OIC and UIC are consistent with each other. We observe that the relative error can be neglected in small $ Q^{2} $ region, but becomes non-negligible at $ Q^{2}\gtrsim1\; \text{GeV}^{2} $ and shows the rapid growing at high $ Q^{2} $ region. At $ Q^{2}=1\; \text{GeV}^{2} $, the relative error is about 7%; at $ Q^{2}=10\; \text{GeV}^{2} $, the error reaches up to 86%. Therefore, WWA has a good accuracy only in very small $ Q^{2} $ domain where exact treatment can nicely recover to WWA. But in the large $ Q^{2} $ region, WWA becomes inapplicable where the deviation is prominent, especially when $ Q^{2}>1\; \text{GeV}^{2} $.

In panels $ (b) $ and $ (c) $, the coherent and ultra-incoherent photon emissions dominate the small and large $ Q^{2} $ regions, respectively. They become comparable at $ Q^{2}=10^{-1} $ and $ 10^{-2}\; \text{GeV}^{2} $ in p-p and $ Pb $-$ Pb $ collisions, respectively. We notice that, except for the region $ 10^{-2}<Q^{2}<10^{-1}\; \text{GeV}^{2} $, the contribution of ordinary-incoherent photon emission can be neglected safely comparing with other two channels. According to the views derived from panel $ (a) $, one can deduce that WWA is valuable in coherent and ordinary incoherent-photon emissions. Especially in $ Pb $-$ Pb $ collisions, the WWA has a much higher accuracy in coherent process which quickly approaches to zero before $ Q^{2}=10^{-2}\; \text{GeV}^{2} $, this effectively avoid WWA error. However, WWA is not valid for ultra-incoherent photon emission, which concentrates on the large $ Q^{2} $ domain where the WWA error is prominent.

In Fig. 6, the results are expressed as a function of y. In panel (a), the curves of coh. and OIC share the same trend, and are consistent with each other in the large y domain. We find that $ y=0.1 $ is a key point: when $ y<0.1 $ the WWA results nicely agree with the exact ones for both channels, and WWA has a better accuracy in coherent process rather than OIC one (at $ y=10^{-4} $, the deviations are about 1.6‰ and 6.7‰ in each case, respectively); but when $ y>0.1 $, the relative errors become prominent, especially the curves show a pronounced peak near $ y=1 $. Therefore, WWA is only effective when $ y<0.1 $ and it has evident error at large values of y. Additionally, the relative error of UIC is prominent (about $ 0.3 $-$ 0.6 $) in the whole y regions.

In panels (b) and (c), the coherent processes dominate the very small y regions, this behaviour guarantees the accuracy of WWA. Conversely, the UIC contributions are predominantly from the regions $ y>10^{-1} $ and $ y>10^{-3} $ in p-p and $ Pb $-$ Pb $ collisions, respectively. Comparing with the views derived from panel (a), we can verify that the validity range of WWA is compatible with coh. and OIC processes. Particularly in $ Pb $-$ Pb $ collisions, the WWA has a better accuracy in coherent process, since its contribution fall off far more rapidly. However, in UIC process, WWA is inapplicable in the entire y regions.

In Fig. 7, we provide the cross section $ d\sigma/dz $ vs z. The panel (a) explicitly shows us that, WWA is applicable in coh. and OIC processes in the entire z range, and has the highest accuracy in coh. process (the deviation is about 1.1‰ − 4.7‰). Whereas the deviation of UIC process reaches up to 153%-268% in the whole z regions. In panels (b) and (c), we find that the z behaviour of each channel is quite different. Firstly, the resolved contributions dominate the lower z region and become smaller than those of direct contributions when $ z>0.3 $; this feature agrees with the traditional perspective that the resolved process contributes appreciably only at $ z\leq0.3 $ [76]. Secondly, the results are divergent near the endpoint $ z=1 $, which are mainly from the color-octet $ ^{1}S_{0} $ and $ ^{3}P_{J} $ processes. The reason is that the NRQCD prediction breaks down and the color-octet channels exhibit collinear singularities in the region of $ z\lesssim1 $, where diffractive production takes place. In order to screen the collinear singularities and suppress the elastic production, the traditional way is to impose the cuts: $ z<0.9 $, $ M_{X'}>10\; \text{GeV} $ and $ p_{T}>1\; \text{GeV} $ or $ M_{J/\psi} $ (actually, if $ p_{T\; \text{min}}\neq0 $, $ z_{\text{max}} $ will naturally less than one). Another possibility to suppress the elastic production at $ z\lesssim1 $ would be to require that $ Q^{2} $ be sufficiently large [26]. However, then also the bulk of the inelastic contribution would be sacrificed. Finally, the coherent contributions are larger than those of UIC processes in the whole z regions. Especially in $ Pb $-$ Pb $ collisions, coherent contribution dominates the entire z regions about an order of magnitude than other two channels. Since the equivalent photon flux scales as $ Z^{2} $, which is a large enhancement factor for the cross section.

In Fig. 8, we show the total cross sections distributed in $ \sqrt{s} $. In panel (a), the curves of coh. and OIC are consistent with each other. We do not show the curve of UIC process, where the deviation is $ 0.38 $-$ 0.49 $ in the whole $ \sqrt{s} $ ranges. It can be seen that the curves show a pronounced rising when $ \sqrt{s}<400 $ and $ 1400\; \text{GeV} $ in p-p and $ Pb $-$ Pb $ collisions, respectively. And they slowly decreased with increasing $ \sqrt{s} $. The trend of curves is similar to the results of Kniehl [22]. Therefore, WWA has the significant errors in small $ \sqrt{s} $ domain (RHIC energies), and has a good accuracy at high energies (LHC energies). In panels (b) and (c), the coherent process is slightly smaller than UIC one in p-p collisions; while the situation is opposite in $ Pb $-$ Pb $ collisions, where the coherent process starts to play a fundamental role (it is about one and two orders of magnitude larger than OIC and UIC processes, respectively). The reason is that the coherent process scale with $ Z^{2} $, whereas the OIC and UIC processes scale approximately with Z and $ N_{A} $, respectively. As for $ Z\gg1 $, the coherent part dominates the production processes. Based on the views derived from the left panel and previous discussions, we can deduce the further conclusion that WWA can reach the high accuracy in high energy range and in $ Pb $-$ Pb $ collisions, where the UIC contribution is greatly suppressed; however it is not a good approximation in p-p collisions, where the UIC contribution is comparable with the coherent one.

In the comprehensive discussions of the results displayed in above four figures [Figs. 5-8], we obtained the validity scopes of WWA. Therefore, in deriving of the equivalent photon spectra based on WWA, the kinematical boundaries are crucial to the accuracy. Now, we display the total cross sections in Tables 1-3, to discuss the accuracies and its sources of the widely employed photon spectra mentioned in Section III. In p-p collisions [Table 1], the relative error of $ f^{\gamma}_{\text{DZ}} $ is the largest, since the integration is performed in the entire kinematical allowed regions: $ Q^{2}_{\text{max}}=\infty $ and $ y_{\text{max}}=1 $, we know that from Figs. 5 and 6, the $ Q^{2}>1\; \text{GeV}^{2} $ and $ y>0.1 $ domains are able to give a large fictitious contribution. For the spectrum $ f^{\gamma}_{\text{Ny}} $, the deviation has a obvious reduction compared to $ f^{\gamma}_{\text{DZ}} $, since it includes the $ Q^{2}_{\text{min}} $ term in Eq. (50) which is omitted in $ f^{\gamma}_{\text{DZ}} $. Actually this term is inversely proportional to $ Q^{4} $ and thus has the noticeable effect in small $ Q^{2} $ region, which can not be neglected when performing the photon spectra; the perspective agrees with Kniehl [22]. The relative error of $ f^{\gamma}_{\text{Kn}} $ is higher than that of $ f^{\gamma}_{\text{Ny}} $ but lower than that of $ f^{\gamma}_{\text{DZ}} $, since the magnetic form factor of proton is also included in this form, which effects only the large $ Q^{2} $ range and should be essentially excluded [56]. Finally, the modified proton spectra $ f^{\gamma}_{\text{MD}} $ nicely agrees with the exact ones ($ \delta\sim1 $%). There are two reasons: $ f^{\gamma}_{\text{MD}} $ is derived from the complete form Eq. (47) which includes the $ Q^{2}_{\text{min}} $ term and properly excludes the effects of magnetic form factor; it adopts the coherence condition, which means that the wavelength of the photon is larger than the size of the nucleus, and the charged constituents inside the nucleus should act coherently. This condition effectively cut the WWA errors. There are also other limitations which can reach the such high accuracy, the key point is that, in most of the physically interesting cases a dynamical cut off exists, such that, the photo-absorption cross sections differ slightly from their values on the mass shell. The details can be found in Ref. [49].

coh.(dir.+res.)	Exact	$ f_{\text{DZ}}^{\gamma} $	$ f_{\text{Ny}}^{\gamma} $	$ f_{\text{Kn}}^{\gamma} $	$ f_{\text{SC}}^{\gamma} $	$ f_{\text{MD}}^{\gamma} $
$ \sigma [\text{nb}] $	22.4241	43.3561	35.1059	39.9058	20.3327	22.6512
δ [%]	0.00	93.35	56.55	77.96	9.33	1.01
Relative error with respect to the exact result: $ \delta=|\sigma/\sigma_{\text{Exact}}-1| $.

Table 1.  $ J/\psi $ photoproduction in the channel of coherent-photon emission in p-p collisions [$ 7\; \text{TeV} $].

coh.(dir.+res.)	Exact	$ f_{\text{DEZ}}^{\gamma} $	$ f_{\text{SC}}^{\gamma} $
$ \sigma [\text{mb}] $	2.4854	18.9783	1.9274
$ \delta [ $%$ ] $	0.00	663.59	22.45

Table 2.  Same as Table 1 but in $ Pb $-$ Pb $ collisions [$ 2.76\; \text{TeV} $].

$ \sigma_{\text{UIC (dir.+res.)}.} $	Exact	$ f_{q}^{\gamma} $	$ f_{\text{BKT}}^{\gamma} $
$ \sigma_{pp}\ [\text{nb}] $	24.2773	39.2290	101.3291
$ \delta_{pp}\ [ $%$ ] $	0.00	61.59	317.38
$ \sigma_{PbPb}\ [\text{mb}] $	0.2494	0.6023	0.7254
$ \delta_{PbPb}\ [ $%$ ] $	0.00	141.51	190.87

Table 3.  Same as Table 1 but for ultra-incoherent photon emission in p-p [$ 7\; \text{TeV} $] and $ Pb $-$ Pb $ [$ 2.76\; \text{TeV} $] collisions.

In $ Pb $-$ Pb $ collisions [Table 2], the deviation of $ f^{\gamma}_{\text{DEZ}} $ reaches up to 680.19%, since $ f^{\gamma}_{\text{DEZ}} $ is constructed on a contradictory assumption: $ Q^{2}_{\text{max}}\sim\infty $ and $ y\ll1 $ (actually, $ Q^{2}_{\text{max}}\sim\infty $ is equivalent to $ y_{\text{max}}\sim1 $). For $ f^{\gamma}_{\text{SC}} $, it roughly agrees with the exact ones in both p-p and $ Pb $-$ Pb $ collisions, but the deviations still can not be neglected. Since $ f^{\gamma}_{\text{SC}} $ is calculated from the semiclassical impact parameter description, where the coherence condition is used. One should be careful that when using $ f^{\gamma}_{\text{SC}} $, setting $ y_{\text{max}}=1 $ will cause the erroneous results.

In the case of ultra-incoherent photon emission [Table 3], the derivations are prominent and turn to much serious in $ Pb $-$ Pb $ collisions. This quantitatively verifies the inapplicability of WWA in UIC process. In particular, the errors of $ f^{\gamma}_{\text{BKT}} $ are the largest, since it is originally derived from e-p scattering, but is directly expanded to describe the probability of finding a photon in any relativistic fermion [20−22], this will overestimate the cross sections. Therefore, the UIC process should be treated in exact treatment, and the results in Refs. [13−28] are not accurate enough, where the mentioned spectra are adopted and the serious double counting exists.

Next, in Table 4, we estimate the contributions of the UIC channel to the heavy quarkonium photoproduction, and discuss the double counting encountered in most of the works when the different photon emissions are considered at the same time. In p-p collisions, UIC result is comparable with coh. one for $ J/\psi $, and becomes about two times larger than coh. one for $ \Upsilon(1S) $. In $ Pb $-$ Pb $ collisions, the ratio of UIC to the coh. one is about 10% for $ J/\psi $, and about 29% for $ \Upsilon(1S) $. Therefore, the UIC channel plays the very important role in p-p collisions, and can also provide the meaningful contribution in $ Pb $-$ Pb $ collisions. In addition, the relative contribution of UIC channel to the inelastic heavy quarkonium photoproduction becomes much more obvious along with the increasing quarkonium mass.

Exact	Quarkonium	coh.	UIC	Total	Total (NWF)
$ \sigma_{pp}\ [\text{nb}] $	$ J/\psi $	22.4241	24.2773	46.7014	103.9845
$ \sigma_{PbPb}\ [\text{mb}] $	$ J/\psi $	2.4854	0.2494	2.7348	3.1055
$ \sigma_{pp}\ [\text{pb}] $	$ \Upsilon(1S) $	13.6794	26.9652	40.6446	82.6189
$ \sigma_{PbPb}\ [\mu\text{b}] $	$ \Upsilon(1S) $	0.7345	0.2113	0.9458	1.1262

Table 4.  Exact results of the $ J/\psi $ and $ \Upsilon(1S) $ photoproductions in p-p [$ 7\; \text{TeV} $] and $ Pb $-$ Pb $ [$ 2.76\; \text{TeV} $] collisions. All the results are the sum of the direct and resolved contributions, and do not include fragmentation contributions.

On the other hand, Table 4 also shows us that the double counting is serious. The total results with no weighting factor (NWF) are much larger than those of exact ones. These large fictitious contributions are mainly from the lower $ Q^{2} $ domain. The traditional way of excluding these unphysical results is to utilize a cutoff $ Q^{2}>1\; \rm{GeV^{2}} $ when performing the equivalent photon spectra. Such as the calculation of $ f_{q}^{\gamma} $ [Eq. 56], the cutoff is adopted to replace the remaining probability factor "$ 1-w_{\text{coh}} $". However, we can see that the corresponding errors in Table 3 are still evident. Therefore, when the different photon emissions are considered simultaneously, each channel should be multiplied by the probability or weighting factor to avoid double counting.

Now, we turn to study the contribution of photoproduction processes to the heavy quarkonium productions in p-p and $ Pb $-$ Pb $ collisions at LHC energies. Before the discussion, there is one critical issue must be addressed: heavy quarkonium productions in hadronic collisions (especially $ Pb $-$ Pb $) are spatially limited, with a collision parameter range within $ 2R_{A} $ ($ R_{A}=A^{1/3}1.2\; \text{fm} $ is the size of the nucleus). When considering the contribution of photonproduction to the hadronic process, the equivalent photon flux should be integrated to $ b_{ \text{max}}=2R_{A} $. Integrating to infinity artificially enhances the cross section by including unphysical contributions from non-overlapping nuclei at large, violating the spatial constraints of the collision geometry. However, in present paper, we performed the calculations and developed our formalism based on the parameterized form factor method, which do not include impact parameter b explicitly (where the form factor is derived from the Fourier transform of the charge density ρ of the nucleus and the approximate plane waves are used). For the consistency of our formalism, we adopt the empirical approximation to include this truncation affects, which is also used in the literature [3, 4, 63]. Where the effect of $ b_{ \text{max}}=2R_{A} $ is appear in the cut $ Q^{2}_{ \text{min}} $, the details see Appendix B. The strict method to incorporate b is the semiclassical impact parameter description, which is more appropriate to calculate the case of heavy ions collisions and often used in the study of UPCs (recently, Wangmei Zhang et al. developed a spatially-dependent photon flux distribution and established a more precise relationship between the photon transverse momentum distribution and collision impact parameter [46]).

In Figs. 9 and 10, we adopt the exact treatment to plot the $ p_{T} $ distributions of charmonium ($ J/\psi $, $ \psi(2S) $) and bottomonium ($ \Upsilon(nS) $) photoproductions, respectively. Where all the results are the sum of the direct and resolved photon contributions. For completely showing the photoproduction contribution and studying the $ p_{T} $ behaviour of each photon source, we plot the Figs. 9 and 10 in different $ p_{T} $ ranges which compensate with each other. In Fig. 9, we plot the $ p_{T} $ distribution of charmonium in the range of $ 2\; \text{GeV}<p_{T}<20\; \text{GeV} $. Fragmentation mechanism is unreliable at this $ p_{T} $ range, because the LP and NLP contributions is valid only for $ p_{T}\gg m_{H} $. (Authors in Ref. [32, 34] have shown that the NLO LP SDC are very different to the fixed-order SDC accurate through NLO at small $ p_{T} $. They suggested that a criterion $ p_{T}>3M_{J/\psi}\simeq10\; \text{GeV} $ to ensure the validity of fragmentation mechanism and suppress NLP contributions.) In order to consider the NLO effects, the LO results are multiplied by the K factor in Ref. [79], we choose $ K=1.8 $ and $ 1.3 $ for direct and resolved contributions, respectively. The spectra of photoproductions are compared to the hard scattering of initial partons. In panels (a) and (c), we compare the results with data of CMS Collaboration [77]. At very small $ p_{T} $, the color-singlet and color-octet distributions are corrupted by collinear divergences and soft gluon effects [31]. Since the intrinsic motion of incident partons renders the differential cross section uncertain for $ p_{T}\leq2\; \text{GeV} $, our predictions only concern the range $ p_{T}\geq2\; \text{GeV} $.

Figure 9.  (color online) The $ p_{T} $ distribution of $ J/\psi $ and $ \psi(2S) $ photoproductions in p-p and $ Pb $-$ Pb $ collisions. Red dashed line denotes the initial partons hard scattering (had.scat.). Magenta dot-dashed line is for coherent-photon emissions [coh.(dir.+res.)]. Blue dot-dot dashed line is for ultra-incoherent photon emissions [UIC (dir.+res.)]. Dark yellow dotted line is for ordinary-incoherent photon emissions [OIC (dir.+res.)]. Black solid line is for the sum of had.scat. and photoproduction processes. The $ J/\psi $ and $ \psi(2S) $ data are from the CMS collaboration [77]. The rapidity $ y_{r} $ is integrated over the experimental range $ |y_{r}|<1.2 $.

Figure 10.  (color online) Same as Fig 9, but for $ \Upsilon(nS) $ productions, all the results are the sum of direct and fragmentation $ \Upsilon(nS) $ productions. The branching fractions are $ B_{r}(\Upsilon(1S)\rightarrow \mu^{+}\mu^{-})=2.48 $%, $ B_{r}(\Upsilon(2S)\rightarrow \mu^{+}\mu^{-})=1.93 $%, and $ B_{r}(\Upsilon(3S)\rightarrow \mu^{+}\mu^{-})=2.18 $%. The $ \Upsilon(nS) $ data are from the CMS collaboration [78]. The rapidity $ y_{r} $ is integrated over the experimental range $ |y_{r}|<1.2 $.

In p-p collisions [panels (a), (c)], the UIC processes are larger than the coh. ones, this verifies again the significance of UIC channel in p-p collisions. However, the coh. contributions are still comparable with UIC ones ($ d\sigma_{\text{coh.}}/d\sigma_{\text{UIC}}\sim0.18-0.13 $ in the $ p_{T} $ ranges considered). This is different to the results in Ref. [18], where the coh. contributions are neglected. In $ Pb $-$ Pb $ collisions [panels (b), (d)], the coh. processes are larger than the UIC ones, since they scale with $ Z^{2} $, whereas the UIC ones scale with $ N_{A} $. We disagree with the results in Refs. [20, 21] where the situation is opposite to us, the UIC contributions are about an order of magnitude larger than coh. ones. Finally, we notice that the charmonium photoproduction processes give the valuable corrections to the had.scat. when $ p_{T}>4\; \text{GeV} $, and the corrections become more obvious along with the increasing $ p_{T} $. And compared to p-p collisions, the corrections are more obvious in $ Pb $-$ Pb $ collisions. This is also different to the the results in Refs. [20, 21] where the contribution of photoproduction even larger than the hadronic process. The reason is that they do not consider the cut $ b_{ \text{max}}=2R_{A} $, artificially enhancing the cross section by including unphysical contributions from non-overlapping nuclei at large. This cut will exclude the large contribution from small $ Q^{2} $ region [see Fig. 5], and causing the coherent contribution becomes much smaller (for incoherent channel, this cut only cause the small effect).

In Fig. 10, we plot the exact results of bottomonium ($ \Upsilon(nS) $) photoproduction and fragmentation processes. For avoiding double counting, we adopt the BCCKL method where the matching rule between fixed-order and fragmentation has been performed [32−34]. Since the fragmentation bottomonium productions are untrustworthy throughout the $ p_{T}<15\; \text{GeV} $ region [31], and in order to suppress the next to leading power (NLP) contributions. We present our theoretical predictions for the values of $ p_{T}>20\; \text{GeV} $. In upper panels, the exact results are compared to the CMS data [78]. In the case of p-p collisions [panels (a)-(c)], the differences between UIC processes and coh. ones become much larger than those of $ J/\psi $ productions in Fig. 9. In the case of $ Pb $-$ Pb $ collisions [panels (d)-(f)], the UIC channel becomes larger than the coh. one at sufficiently large $ p_{T} $ ranges. This feature is opposite to the situation in panels (b), (d) of Fig. 9. We disagree with the results in Ref. [21] where the UIC contributions are much larger than the coh. ones in the whole $ p_{T} $ ranges, and the fragmentation formalism is used beyond its validity scope and the double counting exists. Thus, the relative contribution of UIC channel to the heavy quarkonium photoproduction increases with the growing quarkonium mass. And UIC process begins to dominate the photoproduction processes at large $ p_{T} $ domain. Finally, we find that the contributions of photoproduction and fragmentation processes are still non-negligible in bottomonium production.

In the framework of the NRQCD, we studied the production of heavy quarkonium originating from inelastic photoproduction and fragmentation processes in p-p and $ Pb $-$ Pb $ collisions at LHC energies. We derived the exact treatment by performing a consistent analysis of the terms neglected in transiting from the exact formula to the WWA one, where the Martin-Ryskin method is adopted to weight the different photon sources, and the BCCKL method is adopted to match the fixed-order and fragmentation contributions. The full partonic kinematics matched with the exact treatment were also given. By comprehensively discussing the $ Q^{2} $-, y-, z- and $ \sqrt{s} $-dependence behaviours of cross sections, we obtained the features of WWA in inelastic heavy quarkonium photoproductions in heavy-ion collisions. In order to estimate the errors existing in the widely employed photon spectra, we calculated the total cross sections. Meanwhile, the double counting problems and the relative contributions of the ultra-incoherent channel were also discussed. Finally, we also calculated the $ p_{T} $-dependent cross sections to discuss the contributions of photoproduction and fragmentation processes to the inelastic heavy quarkonium productions.

We found that the inelastic photoproduction and fragmentation processes provide valuable contributions to the heavy quarkonium production, especially in the large $ p_{T} $ regions. While the ultra-incoherent photon emission plays very important role in p-p collisions, and can also provide the meaningful contribution in $ Pb $-$ Pb $ collisions; its relative contribution to inelastic heavy quarkonium photoproductions rapidly increases along with the growing quarkonium mass, and begins to dominate the photoproduction processes at large $ p_{T} $ ranges.

The validity scopes of the WWA in the processes discussed in this work are highly restricted, which are compatible with the characteristics of coh. and OIC processes, and possess high accuracy only within the ranges of $ Q^{2}<1\; \text{GeV}^{2} $, $ y<0.1 $, $ Z\gg1 $, and $ \sqrt{s}>400 $ and $ 1400\; \text{GeV} $ in p-p and $ Pb $-$ Pb $ collisions respectively. In particular, the WWA has a much higher accuracy in $ Pb $-$ Pb $ collisions, where the coh. process is enhanced by $ Z_{Pb}^{2} $ and the UIC one is greatly suppressed. Conversely, the kinematical behaviours of the UIC process are contradict to those of WWA, which concentrate on the regions where the WWA errors are prominent. Therefore, WWA is not a suitable approximation in p-p collisions, where the UIC contribution becomes dominant. Furthermore, we found that the aforementioned equivalent photon spectra are generally derived beyond the applicable scopes of WWA, and the serious double counting exists when different channels are considered simultaneously. Indeed, the exact treatment is necessary to accurately address the inelastic heavy quarkonium photoproduction in p-p and $ Pb $-$ Pb $ collisions at LHC energies.

It should be noted that, when considering the contribution of photonproduction to the hadronic process, the equivalent photon flux should be integrated to $ b_{ \text{max}}=2R_{A} $. Integrating to infinity artificially enhances the cross section by including unphysical contributions from non-overlapping nuclei at large. In present paper, we adopt the empirical approximation to include this truncation affects. The proper method to incorporate b is the semiclassical impact parameter description (Glauber model), which is more appropriate to calculate the case of heavy ions collisions (recent work, see Ref. [46]). In our future work, we will perform our work based on the b-dependent scenario.

For the reader's convenience and for completeness, we list here the involved LDMEs. The mass of heavy quarkonia are $ M_{J/\psi}=3.097\; \text{GeV} $, $ M_{\psi(2S)}=3.686\; \text{GeV} $, $ M_{\Upsilon(1S)}=9.460\; \text{GeV} $, $ M_{\Upsilon(2S)}=10.023\; \text{GeV} $, and $ M_{\Upsilon(3S)}=10.355\; \text{GeV} $ [73]. The LDMEs of the charmonium are given by [79−81],

The LDMEs for the bottomonium [81] are,

And the multiplicity relations

are used, where $ m_{c} $ and $ m_{b} $ are the charm quark and bottom quark mass, respectively.

We give here a detailed treatment of the partonic kinematics which is matched with the exact treatment.

The energy and momentum in α-b parton level read

where the $ s_{\alpha b} $ for each photon emission are

where $ s=\left(p_{A}+p_{B}\right)^2=\left(N_{A}+N_{B}\right)^{2}s_{NN}/4 $ is the energy square in A-B process. While the energy and momentum in $ \gamma^{*} $-b parton level are

The Mandelstam variables involved in the case of direct photoproduction are

while those for resolved photoproduction are

The kinematical boundaries of the mentioned distributions are given in Tables B1, Tables B2. When we calculate the contribution of photoproduction to hadronic process in Figs. 9 and 10, the cut $ b_{ \text{max}}=2R_{A} $ should be considered for avoiding the unphysical contributions from non-overlapping nuclei at large. We adopt the central ideal in Refs. [3, 63], the squared transverse momentum of the photon can be written as

Variables	Coherent direct	UIC direct	Coherent resolved	UIC resolved
$ z_{\text{min}} $	$ \left[(M_{H}^{2}+\hat{s})-\sqrt{(\hat{s}-M_{H}^{2})^{2}-4p_{T \text{min}}^{2}\hat{s}}\right]/(2\hat{s} $)
$ z_{\text{max}} $	$ \left[(M_{H}^{2}+\hat{s})+\sqrt{(\hat{s}-M_{H}^{2})^{2}-4p_{T \text{min}}^{2}\hat{s}}\right]/(2\hat{s} $)
$ \hat{t}_{\text{min}} $	$ -(1-z_{\text{min}})(\hat{s}+Q^{2}) $		$ -(1-z_{\text{min}})\hat{s}^{*} $
$ \hat{t}_{\text{max}} $	$ -(1-z_{\text{max}})(\hat{s}+Q^{2}) $		$ -(1-z_{\text{max}})\hat{s}^{*} $
$ z_{a' \text{min}} $	$ \backslash $	$ \backslash $	$ \hat{s}^{*}_{\text{min}}/[y(s_{\alpha b}-m_{\alpha}^{2})] $	$ \hat{s}^{*}_{\text{min}}/(yx_{a}x_{b}s_{NN}) $
$ z_{a' \text{max}} $	$ \backslash $	$ \backslash $	1	1
$ x_{b \text{min}} $	$ (\hat{s}_{\text{min}}+Q^{2})/[y(s_{\alpha b}|_{x_{b\text{max}}}-m_{\alpha}^{2})] $	$ (\hat{s}_{\text{min}}+Q^{2})/(yx_{a}s_{NN}) $	$ \hat{s}^{*}_{\text{min}}/[y(s_{\alpha b}|_{x_{b\text{max}}}-m_{\alpha}^{2})] $	$ \hat{s}^{*}_{\text{min}}/(yx_{a}s_{NN}) $
$ x_{b \text{max}} $	1
$ x_{a \text{min}} $	$ \backslash $	$ (\hat{s}_{\text{min}}+Q^{2})/(ys_{NN}) $	$ \backslash $	$ \hat{s}^{*}_{\text{min}}/(ys_{NN}) $
$ x_{a \text{max}} $	$ \backslash $	1	$ \backslash $	1
$ y_{\text{min}} $	$ (\hat{s}_{\text{min}}+Q^{2})/(s_{\alpha b}|_{x_{b\text{max}}}-m_{\alpha}^{2}) $		$ \hat{s}^{*}_{\text{min}}/(s_{\alpha b}|_{x_{b\text{max}}}-m_{\alpha}^{2}) $
$ y_{\text{max}} $	$ \left[\sqrt{Q^{2}(4m_{\alpha}^{2}+Q^{2})}-Q^{2}\right]/(2m_{\alpha}^{2}) $
$ Q^{2}_{\text{min}} $	$ -2m_{\alpha}^{2}+\left[(s_{\alpha b}+m_{\alpha}^{2})(s_{\alpha b}-\hat{s}+m_{\alpha}^{2})-(s_{\alpha b}-m_{\alpha}^{2})\sqrt{(s_{\alpha b}-\hat{s}+m_{\alpha}^{2})^{2}-4s_{\alpha b}m_{\alpha}^{2}}\right]/(2s_{\alpha b}) $
$ Q^{2}_{\text{max}} $	$ y(s_{\alpha b}|_{x_{b\text{max}}}-m_{\alpha}^{2})-\hat{s}_{\text{min}} $

Table B1.  The kinematical boundaries of $ Q^{2} $, y distributions. $ \hat{s}_{\text{min}}=\hat{s}^{*}_{\text{min}}=(m_{T}+p_{T\text{min}})^{2} $, $ p_{T}^{2}=\hat{t}(\hat{s}\hat{u}+Q^{2}M_{H}^{2})/(\hat{s}+Q^{2})^{2} $.

Variables	Coherent direct	UIC direct	Coherent resolved	UIC resolved
$ Q^{2}_{\text{min}} $	$ x^{2}_{1}m^{2}_{\alpha}/(1-x_{1}) $
$ Q^{2}_{\text{max}} $	$ 1/R_{\alpha}^{2} $	$ (1-x_{1})s_{NN} $	$ 1/R_{\alpha}^{2} $	$ (1-x_{1})s_{NN} $
$ |y_{r \text{max}}| $	$ \ln[(\hat{s}_{\text{max}}+M_{H}^{2}+\sqrt{(\hat{s}_{\text{max}}-M_{H}^{2})^{2} -4p^{2}_{T}\hat{s}_{\text{max}}})/(\hat{s}_{\text{max}}+M_{H}^{2}-\sqrt{(\hat{s}_{\text{max}}-M_{H}^{2})^{2} -4p^{2}_{T}\hat{s}_{\text{max}}})]^{1/2} $
$ p^{2}_{T\text{min}} $	$ M^{2}_{H} $
$ p^{2}_{T\text{max}} $	$ (1-z)\left[z(s_{\alpha b}|_{x_{b\text{max}}}-m_{\alpha}^{2})-M_{H}^{2}\right] $

Table B2.  The kinematical boundaries of $ p_{T} $, z distributions. Where $ x_{1}=\hat{s}/s_{\alpha b} $, other boundaries are the same as Table B1.

where $ b_{ \text{max}}=2R_{A} $ ($ R_{A}=A^{1/3}1.2\; \text{fm}) $.

Finally, we give here the Jacobian determinant $ {\cal{J}} $ for each distribution. In the case of the $ Q^{2} $ and y distributions, we have

while those for z distribution are

and that for $ p_{T} $ distribution is

for coherent-direct process. And the relations between Eq. (72) and the rest cases are: $ {\cal{J}}_{\text{incoh.dir.}}={\cal{J}}/x_{a} $, $ {\cal{J}}_{\text{coh.res.}}={\cal{J}}/x_{b} $, and $ {\cal{J}}_{\text{incoh.res.}}={\cal{J}}/x_{a}x_{b} $. In the case of fragmentation processes, $ {\cal{J}}=\left(\hat{s}+Q^{2}\right)/\left(\cosh y_{r}\sqrt{\hat{s}}\right) $.