Charmonium pair production in ultraperipheral collision

As the simplest hadron consisting of two heavy quarks, heavy quarkonium provides an ideal laboratory to probe the perturbative and nonperturbative properties of quantum chromodynamics (QCD). Due to the large mass of the constituent quarks, quarkonium production process is assumed to be factorizable into two stages. First, a heavy quark-antiquark pair whose invariant mass near the bound-state mass is produced perturbatively, then the pair binds into a quarkonium state nonperturbatively. The non-relativistic QCD (NRQCD) factorization formalism, which is proposed by Bodwin, Braaten, and Lepage Bodwin:1994jh , provides a systematic framework for the theoretical study of quarkonium production. In NRQCD factorization approach, the nonperturbative hadronization effects are encoded into the long distance matrix elements (LDMEs), which are sorted by means of velocity scaling rules. As a result, the theoretical prediction takes the form of a double expansion in the strong coupling constant $\alpha_{s}$ and the heavy quark relative velocity $v$.

Quarkonium production in photon-photon collisions is a compelling area of study, offering a sensitive probe for testing NRQCD factorization formalism. The cross section for inclusive $J/\psi$ production via photon-photon scattering was measured by the LEP II experiment DELPHI:2003hen , yet the underlying production mechanisms remain incompletely understood. Theoretical researches indicate that the color-singlet (CS) sectors are inadequate to explain the experimental data, while the color-octet (CO) contributions based on different LDME sets diverse quite a lot, which even lead to opposite conclusions in some cases Ma:1997bi ; Japaridze:1998ss ; Klasen:2001mi ; Klasen:2001cu ; Qiao:2003ba ; Li:2009zzu ; Chen:2016hju ; Butenschoen:2011yh . On the other hand, with the advent of high-luminosity colliders, some exclusive processes, like the $J/\psi$-pair production in photon-photon collision, can be measured with high accuracy. In comparison with $J/\psi$ inclusive production, the production mechanism of $J/\psi$-pair is much simpler, due to the fact that the CO contributions, which are suppressed by $v^{8}$, are expected to be insignificant. Therefore, more clear conclusions can be expected from the study of this exclusive process.

In our previous work Yang:2020xkl , we calculated the next-to-leading order (NLO) QCD corrections to the $\gamma+\gamma\to J/\psi+J/\psi$ process, and investigated the corresponding cross section at the SuperKEKB $e^{+}e^{-}$ collider. It should be remarked that the ultraperipheral collision (UPC) of heavy ions is an important alternative test window to study the photon-photon physics. In the UPC, the ion impact parameter is larger than twice the ion radius, hence the ions are kept unbroken during the collision, and the interaction is pure electromagnetic. Due to the coherent action of all the protons in the nucleus, the electromagnetic field is very strong, and the resulting flux of equivalent photons is large. Explicitly, the photon flux from the nuclei scaling as the square of the charge number $Z^{2}$, giving a two-photon luminosity scaling as $Z^{4}$, which eventually leads to a large cross section. In addition, compared with photon-photon collision from electron bremsstrahlung, the UPC event topologies are featured as lacking pileup and large rapidity gap, hence the event signals can be reconstructed with high efficiency. For these reasons, in this work, we investigate the $\gamma+\gamma\to J/\psi+J/\psi$ process at the heavy ion UPC.

The rest of this paper is organized as follows. In Sect. II, we present the primary formulas employed in the calculation. In Sect. III, the total cross sections, their uncertainties, and the differential distributions are discussed. The last section is reserved for summary and conclusions.

The process of quarkonium pair production via photon-photon fusion in an UPC of ions A and B is displayed in FIG. 1. In the framework of equivalent photon approximation (EPA) vonWeizsacker:1934nji ; Williams:1934ad , the total cross section can be expressed as the convolution of the $\gamma+\gamma\to\mathcal{H}+\mathcal{H}$ cross section with the equivalent photon spectral functions

$\displaystyle\sigma(A+B\to A+B+\mathcal{H}+\mathcal{H})$

$\displaystyle=\int\dfrac{dx_{1}}{x_{1}}\dfrac{dx_{2}}{x_{2}}n_{1}(x_{1})n_{2}(x_{2})\hat{\sigma}(\gamma+\gamma\to\mathcal{H}+\mathcal{H}),$

(1)

where $x_{i}=E_{\gamma,i}/E_{\text{beam}}$ is the energy fraction of each photon, $n_{i}(x)$ is the photon spectral function. For ion with charge $Z$, the photon spectral function takes the form Cahn:1990jk

$\displaystyle n(x)=\dfrac{2Z^{2}\alpha}{\pi}\left\{\chi K_{0}(\chi)K_{1}(\chi)-(1-\gamma_{\text{L}}^{-2})\dfrac{\chi^{2}}{2}\left[K_{1}^{2}(\chi)-K_{0}^{2}(\chi)\right]\right\},$

(2)

where $K_{i}$ is the modified Bessel function, $\gamma_{\text{L}}=E_{\text{beam}}/m_{\text{p}}$ is the Lorentz boost factor with $m_{\text{p}}$ denotes the proton mass. The variable $x$ is absorbed in to $\chi=xm_{\text{p}}b_{\text{min}}$, in which the minimum of the impact parameter $b_{\text{min}}$ will be set as the nucleus radius.

In our previous work Yang:2020xkl , we calculated the NLO QCD corrections to the $\gamma+\gamma\to\mathcal{H}+\mathcal{H}$ process. In the following, we briefly review the ingredients of the calculation.

Up to NLO, the $\hat{\sigma}(\gamma+\gamma\rightarrow\mathcal{H}+\mathcal{H})$ of Eq. (1) incorporate the LO contribution and the NLO correction term, that is

$$\hat{\sigma}(\gamma+\gamma\rightarrow\mathcal{H}+\mathcal{H})=\hat{\sigma}_{\rm born}+\hat{\sigma}_{\rm NLO}+\mathcal{O}(\alpha^{2}\alpha^{4}_{s})\ ,$$

(3)

where

$$\begin{split}&\hat{\sigma}_{\rm born}=\frac{1}{2}\frac{1}{2\hat{s}}\int\overline{\sum}|\mathcal{M}_{\rm tree}|^{2}d{\rm PS}_{2}\ ,\\ &\hat{\sigma}_{\rm NLO}=\frac{1}{2}\frac{1}{2\hat{s}}\int\overline{\sum}2{\rm Re}(\mathcal{M}^{*}_{\rm tree}\mathcal{M}_{\rm oneloop})d{\rm PS}_{2}\ .\end{split}$$

(4)

Here, $\frac{1}{2}$ is the statistical factor for identical particles in the final state, $\hat{s}$ is the center-of-mass energy square of the two-photon system, $\overline{\sum}$ means sum (average) over the polarizations and colors of final (initial) state, $d{\rm PS}_{2}$ denotes two-body phase space.

In the calculation of the amplitude, the covariant projection method is applied. The standard form of the spin and color projector is Bodwin:2002cfe :

$$v(p_{\bar{c}})\bar{u}(p_{c})\to\frac{1}{4\sqrt{2}E(E+m_{c})}(\not\!p_{\bar{c}}-m_{c})\not\!\epsilon^{*}_{S}(\not\!P+2E)(\not\!p_{c}+m_{c})\otimes(\frac{\bf{1}_{c}}{\sqrt{N_{c}}})\,$$

(5)

where $\epsilon^{*}_{S}$ denotes the polarization vector of the spin-triplet state, $p_{c}$ and $p_{\bar{c}}$ are the momenta of heavy quark and antiquark respectively, $P=p_{c}+p_{\bar{c}}$ is the momentum of quarkonium, $E=\sqrt{P^{2}/4}$ is the energy of heavy (anti) quark in quarkonium rest frame, $\bf{1}_{c}$ represents the unit color matrix, and $N_{c}=3$ is the colors number in QCD. At the leading order of the relative velocity expansion, it is legitimate to take the approximation of $p_{c}=p_{\bar{c}}=P/2$ and $E=m_{c}$. As for the spin-singlet state, the projection operator can be obtained by replacing the $\not\!\epsilon^{*}_{S}$ with $\gamma_{5}$.

In the calculation of the one-loop amplitude, we encounter ultraviolet (UV) and infrared (IR) singularities. Both of them are regularized by the dimensional regularization with $D=4-2\epsilon$. According to Ref.Bodwin:1994jh , the soft IR singularities are canceled with each other, and the collinear IR singularity in the hard-collinear region will not emerge due to the final charm quarks are massive. The Coulomb IR singularities, where a potential gluon is exchanged between the constituent quarks of a quarkonium, vanish in the dimensional regularization since we set $p_{Q}=p_{\bar{Q}}$ before calculating the Feynman integrals Beneke:1997zp . As for the UV singularities, they are removed by the standard renormalization procedure. The relevant renormalization constants include $Z_{2}$, $Z_{m}$, and $Z_{g}$, corresponding to heavy quark field, heavy quark mass, and strong coupling constant, respectively. We define $Z_{2}$, $Z_{m}$ in the on-mass-shell (OS) scheme, $Z_{g}$ in the modified minimal-subtraction ($\overline{\rm MS}$) scheme. The counter terms are

$$\begin{split}&\delta Z^{\rm OS}_{2}=-C_{F}\frac{\alpha_{s}}{4\pi}\left[\frac{1}{\epsilon_{\rm UV}}+\frac{2}{\epsilon_{\rm IR}}-3\gamma_{E}+3\ln\frac{4\pi\mu^{2}}{m_{Q}^{2}}+4\right],\\ &\delta Z^{\rm OS}_{m}=-3C_{F}\frac{\alpha_{s}}{4\pi}\left[\frac{1}{\epsilon_{\rm UV}}-\gamma_{E}+\ln\frac{4\pi\mu^{2}}{m_{Q}^{2}}+\frac{4}{3}\right],\\ &\delta Z^{\overline{\rm MS}}_{g}=-\frac{\beta_{0}}{2}\frac{\alpha_{s}}{4\pi}\left[\frac{1}{\epsilon_{\rm UV}}-\gamma_{E}+\ln(4\pi)\right].\end{split}$$

(6)

Here, $\mu$ is the renormalization scale, $\gamma_{E}$ is the Euler’s constant, $\beta_{0}=\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{f}$ is the one-loop coefficient of the QCD $\beta$-function, and $n_{f}$ denotes the active quark flavor numbers. The color factors $C_{A}=3,C_{F}=\frac{4}{3}$ and $T_{F}=\frac{1}{2}$.

The input parameters taken in the numerical calculation go as follows:

	$\displaystyle\alpha=1/137.065,\ m_{c}=1.5\pm 0.1\ {\rm GeV},\ m_{\text{p}}=0.9315\ {\rm GeV},$
	$\displaystyle|R_{J/\psi}^{\rm LO}(0)|^{2}=0.528\ {\rm GeV}^{3},\ |R_{J/\psi}^{\rm NLO}(0)|^{2}=0.907\ {\rm GeV}^{3},\ R_{\eta_{c}}(0)=R_{J/\psi}(0).$		(7)

Here, the charm quark mass is set to be roughly one half of the $J/\psi$ or $\eta_{c}$ meson mass. The $J/\psi$ radial wave function at the origin is extracted from its leptonic width

$$\Gamma(J/\psi\to e^{+}e^{-})=\frac{4\alpha^{2}}{9m_{c}^{2}}|R_{H}(0)|^{2}\left(1-4C_{F}\frac{\alpha_{s}(\mu_{0})}{\pi}\right),$$

(8)

with $\mu_{0}=2m_{c}$, $\Gamma(J/\psi\to e^{+}e^{-})=5.55$ keV ParticleDataGroup:2024cfk . According to the heavy quark spin symmetry of NRQCD at the leading order in relative velocity expansion Bodwin:1994jh , here we take $R_{\eta_{c}}(0)=R_{J/\psi}(0)$. Note that, in the LO calculation, LO extraction of the wave function and one-loop formula of the strong coupling are applied; while in the NLO, NLO extraction and two-loop formula of the strong coupling are adopted. The two-loop formula for the running coupling constant takes the form

$$\frac{\alpha_{s}(\mu)}{4\pi}=\frac{1}{\beta_{0}L}-\frac{\beta_{1}\ln L}{\beta^{3}_{0}L^{2}},$$

(9)

where, $L=\ln(\mu^{2}/\Lambda^{2}_{\rm QCD})$, $\beta_{0}=\tfrac{11}{3}C_{A}-\tfrac{4}{3}T_{F}n_{f}$, $\beta_{1}=\frac{34}{3}C^{2}_{A}-4C_{F}T_{F}n_{f}-\frac{20}{3}C_{A}T_{F}n_{f}$, with $n_{f}=4$, $\Lambda_{\rm QCD}=297$ MeV. The one-loop formula can be obtained by dropping the $\beta_{1}$ term of Eq. (9).

We first consider the Pb-Pb collision at the HL-LHC, for which we set $\sqrt{s_{\text{NN}}}=5.52$ TeV, $R_{\text{A}}=7.1$ fm. The corresponding LO and NLO cross sections at different renormalization scale are presented in Table 1, wherein, the central values refer to the results at $m_{c}=1.5$ GeV, the superscripts and subscripts corresponding to the results at $m_{c}=1.4$ GeV and $m_{c}=1.6$ GeV, respectively. It can be seen that, As indicated in the previous work, the NLO corrections are significant, and display different behaviors for $J/\psi$-pair and $\eta_{c}$-pair productions. For $J/\psi$-pair production, the NLO corrections are negative at $\mu=2m_{c}$ and $\mu=\sqrt{4m_{c}^{2}+p_{t}^{2}}$, while positive at $\mu=\sqrt{\hat{s}_{\gamma\gamma}}$;

The interaction of heavy ions at large impact parameters is purely electromagnetical, such an interaction can be treated as real photon-photon fusion where the flux is enhanced by $Z^{2}$, Z is the charge of ion. Compared with photon-photon collision from electron bremsstrahlung, the UPC event topologies are featured as lacking pileup and large rapidity gap, in this way the event signatures can be reconstructed with high efficiency. The generic characteristics of photon-photon fusion in ultraperipheral ion-ion collisions at HL-LHC Bruce:2018yzs ; dEnterria:2022sut and FCC Dainese:2016gch ; FCC:2018vvp energies are collected into TABLE. 2.

The total LO and NLO cross sections for double $J/\psi$ and $\eta_{c}$ via various UPCs are given in TABLE.2, here no any dynamical cut is imposed. As indicated in the previous work Qiao:2001wv ; Yang:2020xkl ; He:2024lrb , the NLO QCD radiative and relativistic corrections both turn out to be significant and negative for $J/\psi$, whereas positive for $\eta_{c}$. The NLO K-factor is around 0.25 (2.35) for $J/\psi$ ($\eta_{c}$) through various ion-ion UPCs.

Supposing the integrated luminosities in TABLE. 2, the produced $J/\psi\mbox{-}J/\psi$ numbers via various heavy ion UPCs at the HL-LHC are 140-194, while the produced numbers can reach 835 for p-p collision due to its relative high luminosity. For $\eta_{c}\mbox{-}\eta_{c}$, the produced numbers increase to 325-456 for heavy ion case, and reach 1980 for p-p collision. In experiment, $J/\psi$ can be reconstructed through its leptonic decays, and $\eta_{c}$ can be reconstructed through its $K\bar{K}\pi$ decay channel. According to Ref. ParticleDataGroup:2024cfk , the branching ratios ${\rm Br}(J/\psi\to l^{+}l^{-}(l=e,\mu))=12\%$ and ${\rm Br}(\eta_{c}\to K\bar{K}\pi)=7.3\%$, and hence the numbers of reconstructed candidates per year would reach 2-12 for $J/\psi$-pair production and 1-10 for $\eta_{c}$-pair. Due to the event topologies for ultraperipheral collision are vary clear, the background from various QCD interactions can be suppressed, hence the experimental investigation are feasible. As the collision energies and luminosities are highly improved at the FCC, the yield events for $J/\psi$ pair will reach 180-200, leaving open opportunity for detailed study of differential distributions, and further provides new test window for NRQCD.

As the number of events corresponding to $\gamma+\gamma\to J/\psi\mbox{-}J/\psi\ (\eta_{c}\mbox{-}\eta_{c})$ are considerable, it is worthy to perform a more elaborate phenomenological analysis. The differential transverse momentum distribution $d\sigma/dp_{T}$ and invariant mass of charmonium pair $d\sigma/dm_{\gamma\gamma}$ for $J/\psi$-pair and $\eta_{c}$-pair productions at the HL-LHC and FCC are shown in Fig. 2, 3 respectively. As the effective photon-photon luminosities decrease dramatically with $W_{\gamma\gamma}$, the main contribution can be only related to small $W_{\gamma\gamma}$ values, therefore we restrict the $p_{T}$ of charmonium to be 0-7 GeV in the transverse momentum distribution, and the $m_{\gamma\gamma}$ of the charmonium pair is set to be 6-10 GeV. The cross sections decease rapidly versus high $p_{T}$, showing a logarithmic dependence of $p_{T}$. The invariant mass distribution reaches the maximum around threshold energy. In FIG.4, we show the rapidity difference $\Delta y$, the absolute value of the difference of charmonium-pair, distributions via Pb-Pb and p-p ultraperipheral collisions. The produced charmonium-pair mainly shares similar rapidity and decreases rapidly with the rapidity difference. While for the region $\Delta y>2$, the $J/\psi-J/\psi$ cross sections at NLO are negative, this may be explained by a large negative loop corrections. The unphysical behavior of the NLO charmonium production cross section in UPC could be cured via resetting the renormalization scale, see TABLE 3. We note that similar negative differential cross section (small $p_{T}$ or large rapidity) can be found in Refs.Lansberg:2010cn ; Feng:2015cba ; ColpaniSerri:2021bla ; Chen:2024dkx , resetting the scale or resummation the logarithms may be useful to alleviate the behavior. A high order correction could also be necessary which deserves further investigation.

The fully charm tetraquark state X(6900), which is recently observed by LHCb Collaboration LHCb:2020bwg (another two similar resonances around 6552 and 7287 MeV are established by CMS Collaboration CMS:2023owd ), can be also produced via photon-photon fusion in the UPC and then decays into double charmonium, hence contributes in the phenomenological analysis of double charmonium final state. Using the Low formula Low:1960wv ; Budnev:1975poe ; Shao:2022cly , the production cross section of X(6900) via UPC can be written in terms of the two-photon decay width $\Gamma_{X(6900)\to\gamma\gamma}$ and the effective luminosity $\frac{d\mathcal{L}_{\gamma\gamma}^{AB}}{dW_{\gamma\gamma}}$ at the resonance mass,

$\displaystyle\sigma(AB\xrightarrow{\gamma\gamma}AB+X(6900))=4\pi^{2}(2J+1)\frac{\Gamma_{X\to\gamma\gamma}}{m^{2}_{X}}\left.\frac{d\mathcal{L}_{\gamma\gamma}^{AB}}{dW_{\gamma\gamma}}\right|_{W_{\gamma\gamma}=m_{X}},$

(10)

where J is the spin of X(6900), we take J = 0 in the following discussion. The effective luminosity is calculated by

$\displaystyle\frac{d\mathcal{L}_{\gamma\gamma}^{AB}}{dW_{\gamma\gamma}}=\frac{2}{W_{\gamma\gamma}}\int_{\frac{W^{2}_{\gamma\gamma}}{s_{NN}}}^{1}\frac{dx}{x}n_{\gamma/A}(x)n_{\gamma/B}(\frac{W^{2}_{\gamma\gamma}}{xs_{NN}}).$

(11)

For Pb-Pb collision with nucleon-nucleon c.m. energy $\sqrt{s_{NN}}=5.52\ \rm TeV$, the effective two photon luminosity $\left.\frac{d\mathcal{L}_{\gamma\gamma}^{PbPb}}{dW_{\gamma\gamma}}\right|_{W_{\gamma\gamma}=6.9\ \rm GeV}=9.303\times 10^{3}\ \rm GeV^{-1}$, the estimated cross sections are given in TABLE 4. Supposing that X(6900) mainly decays into di-charmonium, the produced di-charmonium is much great than the direct $J/\psi J/\psi+\eta_{c}\eta_{c}$ yields (93.1 nb). This abnormal behavior may be related to the overestimation of di-photon and di-charmonium decay fractions, which deserves further investigation.

In this work, we investigate the double charmonium production via photon-photon fusion with various ultraperipheral ion-ion collisions in the framework of NRQCD factorization formalism at NLO accuracy. Both $J/\psi$ and $\eta_{c}$ are considered in color-singlet configuration. The total cross sections and differential distributions versus transverse momentum, invariant mass and rapidity difference at the HL-LHC and FCC are given. The theoretical uncertainties caused by charm quark mass and renormalization scale are also estimated.

Numerical results show that the NLO QCD corrections are large but negative for $J/\psi$ pair, while positive for $\eta_{c}$ pair. Based on designed luminosities of each ion-ion collision at the HL-LHC and FCC, a considerable number of double charmonium events can be expected. Due to the event topologies for ultraperipheral collision are vary clear, the background from various QCD interactions can be suppressed, hence the detailed experimental investigation are feasible. The production rate of exotic fully-charmed tetraquark X(6900) through Pb-Pb UPC is also discussed via Low formula with different diphoton widths approximation.

Acknowledgments We thank Long-Bin Chen for useful discussion. This work is supported by the National Natural Science Foundation of China (NSFC) under the Grants Nos. 12205061, 12275185 and 12335002.