Double quarkonium hadroproduction as a probe of gluon Sivers function

TL;DR

This paper predicts single spin asymmetry in double J/psi production in proton-proton collisions using non-relativistic QCD and gluon Sivers function, suggesting measurable asymmetries for probing gluon spin effects.

Contribution

It introduces a prediction of SSA in double J/psi production based on the gluon Sivers function fitted from PHENIX data, incorporating color singlet and octet states.

Findings

Sizable asymmetries estimated as functions of kinematic variables.

Results suggest SSA can probe gluon Sivers function.

Uses non-relativistic QCD framework with recent gluon Sivers fit.

Abstract

The prediction of single spin asymmetry (SSA) in double $\rm J/\psi$ production in proton-proton collision is given out within the framework of non-relativistic QCD, using the recently obtained best fit parameters for the gluon Sivers function extracted from PHENIX data in $\rm p+p^\uparrow\to\pi^0+X$. The color singlet state $\rm ^3 S_1^{(1)}$ and color octet state $\rm ^3 S_1^{(8)}$ are considered to the SSA contribution in the double $\rm J/\psi$ hadroproduction. Our result shows that a sizable asymmetry can be estimated as functions of different kinematic variables.

Tables

Table 1. Table 1: Best fit parameters of Sivers function.

	$N_a$	$\alpha$	$\beta$	$M_1^2$ GeV2	$\rho$	$⟨k_{⟂}^2⟩$ GeV2	Notation
g	0.65	2.8	2.8		0.687	0.25	SIDIS1
g	0.05	0.8	1.3		0.576	0.25	SIDIS2
u	0.18	1.0	6.6	0.8		0.57	BV-a
d	-0.52	1.9	10.0	0.8		0.57	BV-b

Equations

$${\mathrm{A}}_{\mathrm{N}}=\frac{\mathrm{d}{\sigma }^{\uparrow }-\mathrm{d}{\sigma }^{\downarrow }}{\mathrm{d}{\sigma }^{\uparrow }+\mathrm{d}{\sigma }^{\downarrow }}$$

$$\mathrm{p}({\mathrm{P}}_{\mathrm{a}})+\mathrm{p}({\mathrm{P}}_{\mathrm{b}})\to {\mathcal{Q}}_1({\mathrm{p}}_1)+{\mathcal{Q}}_2({\mathrm{p}}_2)+\mathrm{X}$$

$$\frac{\mathrm{d}\sigma }{\mathrm{d\Omega dy}{\mathrm{d}}^2{\mathbf{q}}_{\mathrm{T}}}$$

$${\delta }^4({\mathrm{k}}_{\mathrm{a}}+{\mathrm{k}}_{\mathrm{b}}-{\mathrm{p}}_1-{\mathrm{p}}_2)=\frac{2}{\mathrm{s}}\delta ({\mathrm{x}}_{\mathrm{a}}-\frac{\mathrm{M}{\mathrm{e}}^{\mathrm{y}}}{\sqrt{\mathrm{s}}})\delta ({\mathrm{x}}_{\mathrm{b}}-\frac{\mathrm{M}{\mathrm{e}}^{-\mathrm{y}}}{\sqrt{\mathrm{s}}}){\delta }^2({\mathbf{k}}_{\perp \mathrm{a}}+{\mathbf{k}}_{\perp \mathrm{b}}-{\mathbf{q}}_{\mathrm{T}}).$$

$${\mathrm{x}}_{\mathrm{a}}=\frac{\mathrm{M}{\mathrm{e}}^{\mathrm{y}}}{\sqrt{\mathrm{s}}}{\mathrm{x}}_{\mathrm{b}}=\frac{\mathrm{M}{\mathrm{e}}^{-\mathrm{y}}}{\sqrt{\mathrm{s}}}$$

$${\mathrm{P}}_{\mathrm{a}}=\frac{\sqrt{\mathrm{s}}}{2}(1,0,0,1),{\mathrm{P}}_{\mathrm{b}}=\frac{\sqrt{\mathrm{s}}}{2}(1,0,0,-1).$$

$${\mathrm{P}}_{\mathrm{a}}^{\mu }=\sqrt{\frac{\mathrm{s}}{2}}{\mathrm{n}}_{+}^{\mu },{\mathrm{P}}_{\mathrm{b}}^{\mu }=\sqrt{\frac{\mathrm{s}}{2}}{\mathrm{n}}_{-}^{\mu },$$

$${\mathrm{n}}_{+}^{\mu }=\frac{1}{\sqrt{2}}(1,0,0,1),{\mathrm{n}}_{-}^{\mu }=\frac{1}{\sqrt{2}}(1,0,0,-1),$$

$${\mathrm{p}}_1^{\mu }=\left(\sqrt{{\mathrm{m}}_1^2+{\mathrm{p}}_1^2},{\mathrm{p}}_1\sin {\theta }_1\cos {\varphi }_1,{\mathrm{p}}_1\sin {\theta }_1\sin {\varphi }_1,{\mathrm{p}}_1\cos {\theta }_1\right)$$

$${\mathrm{q}}^{\mu }=\left(\sqrt{{\mathrm{M}}^2+{\mathrm{q}}_{\mathrm{T}}^2}\cosh \mathrm{y},{\mathrm{q}}_{\mathrm{T}}\cos \varphi ,{\mathrm{q}}_{\mathrm{T}}\sin \varphi ,\sqrt{{\mathrm{M}}^2+{\mathrm{q}}_{\mathrm{T}}^2}\sinh \mathrm{y}\right).$$

$${\mathrm{k}}_{\mathrm{a}}^{\mu }$$

$${\mathrm{k}}_{\perp \mathrm{b}}^{\mu }$$

$$\int \delta ({\mathrm{q}}^2-2q{\mathrm{p}}_1)\mathrm{f}({\mathrm{M}}^2)\mathrm{d}{\mathrm{M}}^2$$

$$\displaystyle\rm\bigg{(}t_{0}^{2}-q_{T}^{2}-2\left[t_{0}\Big{(}\cosh y\sqrt{m_{1}^{2}+p_{1}^{2}}-\sinh y\,p_{1}\cos\theta_{1}\Big{)}-q_{T}p_{1}\sin\theta_{1}\cos(\phi-\phi_{1})\right]\bigg{)}=0$$

$$\hat{\mathrm{s}}$$

$$\hat{\mathrm{t}}$$

$$\hat{\mathrm{u}}$$

$$\mathrm{d}{\sigma }^{\uparrow }-\mathrm{d}{\sigma }^{\downarrow }$$

$$\mathrm{d}{\sigma }^{\uparrow }+\mathrm{d}{\sigma }^{\downarrow }$$

$${\mathrm{\Delta }}^{\mathrm{N}}{\mathrm{f}}_{g/{\mathrm{p}}^{\uparrow }}({\mathrm{x}}_{\mathrm{b}},{\mathbf{k}}_{\perp \mathrm{b}})={\mathrm{f}}_{g/{\mathrm{p}}^{\uparrow }}({\mathrm{x}}_{\mathrm{b}},{\mathbf{k}}_{\perp \mathrm{b}})-{\mathrm{f}}_{g/{\mathrm{p}}^{\downarrow }}({\mathrm{x}}_{\mathrm{b}},{\mathbf{k}}_{\perp \mathrm{b}})={\mathrm{\Delta }}^{\mathrm{N}}{\mathrm{f}}_{g/{\mathrm{p}}^{\uparrow }}({\mathrm{x}}_{\mathrm{b}},{\mathrm{k}}_{\perp \mathrm{b}})\hat{\mathbf{S}}\cdot (\hat{\mathbf{P}}\times {\hat{\mathbf{k}}}_{\perp \mathrm{b}}).$$

$${\mathrm{f}}_{g/p}({\mathrm{x}}_{\mathrm{a}},{\mathbf{k}}_{\perp \mathrm{a}})={\mathrm{f}}_{g/p}({\mathrm{x}}_{\mathrm{a}})\frac{1}{\pi ⟨{\mathrm{k}}_{\perp \mathrm{a}}^2⟩}{\mathrm{e}}^{-{\mathbf{k}}_{\perp \mathrm{a}}^2/⟨{\mathrm{k}}_{\perp \mathrm{a}}^2⟩}.$$

$$\hat{\mathbf{S}}\cdot (\hat{\mathbf{P}}\times {\hat{\mathbf{k}}}_{\perp \mathrm{g}})=\sin ({\varphi }_{\mathrm{b}}-{\varphi }_{\mathrm{s}}).$$

$${\mathrm{\Delta }}^{\mathrm{N}}{\mathrm{f}}_{g/{\mathrm{p}}^{\uparrow }}({\mathrm{x}}_{\mathrm{b}},{\mathrm{k}}_{\perp \mathrm{b}})=2{\mathcal{N}}_{\mathrm{g}}({\mathrm{x}}_{\mathrm{b}}){\mathrm{f}}_{g/p}({\mathrm{x}}_{\mathrm{b}})\mathrm{h}({\mathrm{k}}_{\perp \mathrm{b}})\frac{{\mathrm{e}}^{-{\mathbf{k}}_{\perp \mathrm{b}}^2/⟨{\mathrm{k}}_{\perp \mathrm{b}}^2⟩}}{\pi ⟨{\mathrm{k}}_{\perp \mathrm{b}}^2⟩},$$

$${\mathcal{N}}_{\mathrm{g}}({\mathrm{x}}_{\mathrm{b}})={\mathrm{N}}_{\mathrm{g}}{\mathrm{x}}_{\mathrm{b}}^{\alpha }(1-{\mathrm{x}}_{\mathrm{b}}{)}^{\beta }\frac{(\alpha +\beta {)}^{(\alpha +\beta )}}{{\alpha }^{\alpha }{\beta }^{\beta }}.$$

$$\mathrm{h}({\mathrm{k}}_{\perp \mathrm{b}})=\sqrt{2e}\frac{{\mathrm{k}}_{\perp \mathrm{b}}}{{\mathrm{M}}_1}{\mathrm{e}}^{-{\mathbf{k}}_{\perp \mathrm{b}}^2/{\mathrm{M}}_1^2}.$$

$$\mathrm{h}({\mathrm{k}}_{\perp \mathrm{b}})\frac{{\mathrm{e}}^{-{\mathbf{k}}_{\perp \mathrm{b}}^2/⟨{\mathrm{k}}_{\perp \mathrm{b}}^2⟩}}{\pi ⟨{\mathrm{k}}_{\perp \mathrm{b}}^2⟩}=\frac{\sqrt{2e}}{\pi }\sqrt{\frac{1-\rho }{\rho }}{\mathrm{k}}_{\perp \mathrm{b}}\frac{{\mathrm{e}}^{-{\mathbf{k}}_{\perp \mathrm{b}}^2/\rho ⟨{\mathrm{k}}_{\perp \mathrm{b}}^2⟩}}{⟨{\mathrm{k}}_{\perp \mathrm{b}}^2{⟩}^{3/2}},$$

$$\rho =\frac{{\mathrm{M}}_1^2}{⟨{\mathrm{k}}_{\perp \mathrm{b}}^2⟩+{\mathrm{M}}_1^2}.$$

$$(\mathrm{a}){\mathcal{N}}_{\mathrm{g}}({\mathrm{x}}_{\mathrm{b}})$$

$$(\mathrm{b}){\mathcal{N}}_{\mathrm{g}}({\mathrm{x}}_{\mathrm{b}})$$

$$\mathrm{d}{\sigma }^{\uparrow }-\mathrm{d}{\sigma }^{\downarrow }$$

$$\mathrm{d}{\sigma }^{\uparrow }+\mathrm{d}{\sigma }^{\downarrow }$$

$${\mathrm{A}}_{\mathrm{N}}^{\sin (\varphi -{\varphi }_{\mathrm{s}})}=2\frac{\int \mathrm{d}\varphi [\mathrm{d}{\sigma }^{\uparrow }-\mathrm{d}{\sigma }^{\downarrow }]\sin (\varphi -{\varphi }_{\mathrm{s}})}{\int \mathrm{d}\varphi [\mathrm{d}{\sigma }^{\uparrow }+\mathrm{d}{\sigma }^{\downarrow }]}$$

Full text

I INTRODUCTION

Since it was first observed, single spin asymmetry (SSA) is a topic in spin physics of significant theoretical and experimental interest Adams et al. (1991a, b); Arsene et al. (2008). SSA appears in scattering process when one of the colliding proton or the target is transversely polarized with respect to the scattering plane. It can provide information on the three-dimensional structure of the nucleons. From the theoretical point of view, to explain the SSA requires the nonperturbative correlators of quark or gluon and there are two methods for it. The first one is the generalized parton model (GPM) Ji et al. (2004, 2005), where the inclusive cross section is written as a convolution of Transverse Momentum Dependent Partonic Distribution Functions (TMD-PDFs), Transverse Momentum Dependent Fragmentation Functions (TMD-FFs) and QCD partonic cross sections. This method is phenomenologically well studied in refs Ji et al. (2004, 2005); Echevarria et al. (2012); Bacchetta et al. (2007); Anselmino et al. (2003); Boer (1999); Arnold et al. (2009); Boer et al. (1997); Anselmino et al. (2007). The second approach describes the SSA as a twist-3 effect in the collinear factorization, and is suited for describing SSA in the large $\rm p_{T}$ region. This formalism was originally proposed and further developed by Efremov and Teryaev (1985); Qiu and Sterman (1999); Kanazawa and Koike (2000); Kouvaris et al. (2006); Eguchi et al. (2007); Kanazawa et al. (2014). TMD factorization extends collinear factorization by accounting for the parton transerve momentum, generally denoted by $\rm{\bf k}_{T}$. It applies to processes in which a momentum transfer is much larger than $\rm{\bf k}_{T}$, for instance at the LHC when a pair of particles (e.g. two quarkonium states $\rm\mathcal{Q}$) is produced with a large invariant mass ($\rm M_{\mathcal{Q}\mathcal{Q}}$) or a large individual transverse momenta as compared to the sum transverse momenta ($\rm P_{\mathcal{Q}\mathcal{Q}T}$) Lansberg et al. (2018).

Among the single spin asymmetries, the Sivers asymmetry plays a vital important role and it is well studied both theoretically and experimentally. The Sivers function Sivers (1990) represents an azimuthal dependence of the number density of unpolarized quarks inside a transversely polarized proton, and has been measured at HERMES Airapetian et al. (2005, 2009, 2014), COMPASS Adolph et al. (2012, 2014, 2017); Aghasyan et al. (2017), JLAB Qian et al. (2011); Zhao et al. (2014) and RHIC Adamczyk et al. (2016) respectively. It has been found that the initial and final state interactions (gauge links) contribute to the Sivers asymmetries significantly, therefore, they are process dependent Boer et al. (2003). For example, the Sivers function probed in semi-inclusive deep inelastic scattering (SIDIS) is expected to be the same in magnitude but opposite in sign compared to the one probed in the Drell-Yan (DY) process. Among the Sivers functions, the quark Sivers functions have been widely studied over the years Anselmino et al. (2017a), however, the gluon Sivers function (GSF) still remains poorly measured. An indirect estimation of the GSF was obtained, within the GPM framework in Ref.D’Alesio et al. (2015), by fitting the midrapidity data on SSA in $\rm\pi^{0}$ production at RHIC.

Quarkonium processes, through both single and double productions, can be used to probe gluons inside hadron Brambilla et al. (2011). For example, through single quarkonium production in electron-proton ($\rm ep$) Godbole et al. (2012); Mukherjee and Rajesh (2017); Anselmino et al. (2017b); Boer et al. (2016) and proton-proton ($\rm pp$) Anselmino et al. (2004); Godbole et al. (2017); D’Alesio et al. (2017) collision, the GSFs are comprehensively studied. The formation of quarkonium out of the two heavy quarks is a nonperturbative process and is treated in terms of different models where the non-relativistic QCD (NRQCD) Bodwin et al. (1995) factorization is one of them. It has been successfully used to explain the $\rm J/\psi$ hadroproduction at Tevatron Abe et al. (1997); Acosta et al. (2005), also data from $\rm J/\psi$ photoproducton at HERA Aaron et al. (2010); Chekanov et al. (2003); Abramowicz et al. (2013). In the NRQCD, the production and decay of heavy quarkonium are factorized into two stages. First, a heavy quark-antiquark pair is perturbatively built at short distances, which is worked out by expansion in the strong coupling constant $\rm\alpha_{s}$. Then, the pair nonperturbatively evolves into quarkonium at a long distance. The short distance coefficients are calculated perturbatively by the projection technique and the long distance matrix elements (LDMEs) are extracted from the experimental data. The LDMEs behave as powers of v, which is the typical heavy-quark (or antiquark) velocity in the quarkonium rest frame Lepage et al. (1992). Hence, the NRQCD factorization expression can be thought of as double expansions in terms of $\rm v$ as well as $\rm\alpha_{s}$. In fact, the asymmetry is very sensitive to the production mechanism.

Sivers effect has been studied theoretically in several processes, i.e., in $\rm ep^{\uparrow}$ collision, the heavy quark pair and dijet production Boer et al. (2016), the inelastic $\rm J/\psi$ photoproduction Rajesh et al. (2018), $\rm e+p^{\uparrow}\to e+J/\psi+X$ Godbole et al. (2012, 2013, 2015), and also in $\rm pp^{\uparrow}$ scattering, i.e., $\rm pp^{\uparrow}\to J/\psi+X$ Godbole et al. (2017), D-meson production Anselmino et al. (2004); Godbole et al. (2016) and back to back jet correlations Boer and Vogelsang (2004), etc. It is found that some processes are not safe in measuring gluon Sivers function due to the problem of TMD factorization breaking contributions Rogers and Mulders (2010), though in some aspects they do probe TMDs. On the other hand, among the hadronic collisions the processes having one or two color singlets in the final state would, in any cases, be safe D’Alesio et al. (2015) in measuring the GSF. Therefore, both single and double heavy quarkonium productions are considered to be clean as probes of the GSF.

In this work, we investigate feasibility of using double charmonium production to obtain information on the Sivers function and present predictions for SSA through the process $\rm pp^{\uparrow}\to J/\psi+J/\psi+X$. Then we estimate the asymmetry by using NRQCD framework with both color singlet and color octet contributions in pp collision. We provide predictions on asymmetry for future proposed experiments at AFTER@LHC which belong to fixed target experiments with $\rm\sqrt{s}=115$ GeV, and also for $\rm\sqrt{s}=200,\ 500$ GeV which will be explored at the RHIC. Two recent extractions D’Alesio et al. (2015); Anselmino et al. (2017a) were used for the gluon Sivers function from the SSA data in the pp collision at the RHIC. The paper is organized as follows. SSA of double quarkonium production is presented in Sec.II. In Sec.III we present the numerical results along with the conclusion in Sec.IV.

II theoretical framework

Single spin asymmetry for the inclusive process $\rm A^{\uparrow}+B\to C+X$ is described as

[TABLE]

where $\rm d\sigma^{\uparrow(\downarrow)}$ denotes the differential cross section for scattering of a transversely polarized hadron A off an unpolarized hadron B, with A upwards (downwards) transversely polarized with respect to the scattering plane.

II.1 Model for $J/\psi$ pair production

We consider first the inclusive production of a quarkonium pair in unpolarized proton-proton scattering

[TABLE]

where $\rm P_{i}$ and $\rm p_{i}$ are the four momenta of the particles given in parentheses. We assume that the heavy quark-antiquark pairs produced in the final state are in a bound state described by a nonrelativistic wave function with spin S=1, orbital angular momentum L=0 and total angular momentum J=1. In the following, we adopt the spectroscopic notation $\rm\mathcal{Q}=Q\bar{Q}[^{2S+1}L_{J}^{(1,8)}]$ with $\rm Q=c$, where the color assignments for the quark pair are generally specified by the singlet or octet superscripts, (1) or (8). The leading contribution must be the color-singlet channel $\rm\mathcal{Q}[^{3}S_{1}^{(1)}]+\mathcal{Q}[^{3}S_{1}^{(1)}]$. The color octet channel $\rm\mathcal{Q}[^{3}S_{1}^{(8)}]+\mathcal{Q}[^{3}S_{1}^{(8)}]$ is suppressed by a velocity factor compared to the color singlet channel. Nevertheless, the color octet channel can also contribute with a large enhancement factor in some kinematic condition Ko et al. (2011). While the remaining color octet channels suppressed by velocity factors do not have the enhancement factor. Therefore, the contributions to the quarkonium pair production that we consider in this work are $\rm\mathcal{Q}[^{3}S_{1}^{(1)}]+\mathcal{Q}[^{3}S_{1}^{(1)}]$ and $\rm\mathcal{Q}[^{3}S_{1}^{(8)}]+\mathcal{Q}[^{3}S_{1}^{(8)}]$. This consideration is also suitable for the polarized case.

TMD factorization requires restriction to the so-called correlation limit where the sum transverse momentum of the $\rm J/\psi$’s is small compared to the individual transverse momenta. The individual $\rm p_{T}$ of one of the $\rm J/\psi$ (for example, $\rm p_{1T}$), needs to be large and considered fixed by 5 GeV ($\rm\gg\Lambda_{QCD}$) to avoid evolution as this is related to the hard scale, and the transverse momentum of the other $\rm J/\psi$ ($\rm p_{2T}$) is restricted by $\rm p_{2T}\gg\Lambda_{QCD}$. For the sum transverse momentum of the $\rm J/\psi$’s we have $\rm q_{T}\sim\Lambda_{QCD}$, where $\rm\Lambda_{QCD}$ is the QCD scale.

In order to guarantee the TMD factorization, the cross section of double $J/\psi$ production should be differential in the sum transverse momentum($\rm P_{\mathcal{Q}\mathcal{Q}T}$) den Dunnen et al. (2014). Within the generalized parton model formalism, the differential cross section for the process is given by

[TABLE]

where $\rm{\bf k}_{\perp a}$ and $\rm{\bf k}_{\perp b}$ are the transverse momentum of two initial gluons, q and y are the transverse momentum and rapidity of the final system respectively. We have written the four momentum conservation delta function by

[TABLE]

The fractions $\rm x_{a}$ and $\rm x_{b}$ of the initial protons’ longitudinal momenta are given by

[TABLE]

where M is the invariant mass of the final mesons. In addition, we define the four-momentum and mass of $\rm\mathcal{Q}_{1}$ as $\rm p_{1}$ and $\rm m_{1}$, the corresponding solid angle is denoted by $\Omega$. The unpolarized gluon TMD, $\rm f_{g/p}$, represents the density of gluons inside an unpolarized proton. $\rm\mathcal{\overline{M}}_{g+g\rightarrow J/\psi+J/\psi}[^{3}S_{1}^{(1)},^{3}S_{1}^{(1)}]$ ($\rm\mathcal{\overline{M}}_{g+g\rightarrow J/\psi+J/\psi}[^{3}S_{1}^{(8)},^{3}S_{1}^{(8)}]$) is the amplitude of gluon-gluon fusion process with two final color singlet(octet) and its square is given in Appendix.A. The delta function $\rm\delta(q^{2}-2q\cdot p_{1})$ is solved in Sec.II.2. In order to focus on the contribution from initial gluons we impose a cut: $\rm sx_{a}x_{b}=\hat{s}\gg\Lambda_{QCD}$, since the gluon distribution function is much larger than the quark distribution function in the case with relative large x Qiu et al. (2011). As a result, the parton process with $\rm q\bar{q}$ initial states is neglected and we only consider gg initial state in the present paper.

II.2 Kinematics

We consider a frame in which the unpolarized and polarized protons are moving along $\rm z$ and -$\rm z$-axes respectively. The four momenta of the protons are given by

[TABLE]

The center of mass (CM) energy of proton-proton system is $\rm\hat{s}=(P_{a}+P_{b})^{2}=M^{2}$. The above four momenta in light-cone coordinate system can be written as

[TABLE]

where $\rm n_{+}$ and $\rm n_{-}$ are two light-like vectors with $\rm n_{+}.n_{-}=1$ and $\rm n_{+}^{2}=n^{2}_{-}=0$, and given by

[TABLE]

where $\rm x_{a}=\frac{k_{a}^{+}}{P_{a}^{+}}$ is the light-cone momentum fraction. The four momentum of one $\rm\rm J/\psi$ is given by

[TABLE]

and the sum 4-momentum of double $\rm J/\psi$ is parameterized as

[TABLE]

Then the gluons four momenta are given by

[TABLE]

By using the above relations, we can work out the delta function $\rm\delta(q^{2}-2q\cdot p_{1})$ as

[TABLE]

where $\rm t_{0}$ satisfies

[TABLE]

and $\rm f(M^{2})$ is a function with respect to $\rm M^{2}$. Then we can write down the expressions of Mandelstam variables as below

[TABLE]

II.3 The Sivers effect

Now, we are in the position to write down the expression of numerator and denominator terms of Eq.(1) when the target proton is polarized. The expression reads

[TABLE]

and

[TABLE]

where $\rm\Delta^{N}f_{g/p^{\uparrow}}(x_{g},{\bm{k}}_{\perp g})$, GSF, represents the parton density function of gluon inside the transversely polarized proton and is defined as below

[TABLE]

In estimating the SSA numerically, we have to discuss about the parameterization of TMDs. In this paper, we assume that the unpolarized gluon TMD follows the general Gaussian distribution. The Gaussian parameterization of unpolarized TMD is

[TABLE]

Here, $\rm x_{a}$ and $\rm k_{\perp a}$ dependencies of the TMD are factorized. The collinear PDF is denoted by $\rm f_{g/p}(x_{a})$ which is measured at the fixed scale $\rm\mu=\sqrt{m_{1}^{2}+p_{1T}^{2}}$.

We choose a frame as discussed in the last subsection where the polarized proton is moving along -z axis with momentum $\rm{\bf P}$, and transversely polarized $\rm\hat{\bf S}=(\cos\phi_{s},\sin\phi_{s},0)$. The transverse momentum of the initial polarized gluon is $\rm{\bf k}_{\perp b}=k_{\perp b}(\cos\phi_{b},\sin\phi_{b},0)$, then

[TABLE]

We have taken $\rm\phi_{s}=\pi/2$ for numerical estimation. The parameterization of GSF is given by D’Alesio et al. (2015); Anselmino et al. (2017a)

[TABLE]

here $\rm f_{g/p}(x_{b})$ is the usual collinear gluon PDF and

[TABLE]

The definition of function $\rm h(k_{\perp b})$ is given by

[TABLE]

The $\rm k_{\perp b}$ dependent part of Sivers function can be rewritten as

[TABLE]

where we defined

[TABLE]

At RHIC, the GSF from pion production data has been extracted by D’Alesio et al. D’Alesio et al. (2015), and we define two sets of two best fit parameters obtained as SIDIS1 and SIDIS2. In addition, the quark and anti-quark Sivers functions have been extracted by Anselmino et al. Anselmino et al. (2017a) using the latest SIDIS data. However, GSF has not been extracted yet from SIDIS data. Therefore, in order to estimate the asymmetry, we have to use some parameterizations Boer (1999) listed below to reach the best fit parameters of gluon Sivers function

[TABLE]

We call the parameterization (a) and (b) as BV-a and BV-b respectively. The best fit parameters are tabulated in Table 1.

Generally, we take $\rm\langle k_{\perp a}^{2}\rangle=\langle k_{\perp b}^{2}\rangle=\langle k_{\perp}^{2}\rangle$, then the final expressions of numerator and denominator terms of Eq.(1) are given by

[TABLE]

and

[TABLE]

Then we reach the Sivers asymmetry

[TABLE]

where $\phi$ is the azimuthal angle of double $\rm J/\psi$ system.

III Numerical results

In this section, we present predictions of transverse single spin asymmetry in $\rm p+p^{\uparrow}\to J/\psi+J/\psi+X$, obtained using the recent direct fits D’Alesio et al. (2015) and the BV models Boer and Vogelsang (2004) of the GSF with the corresponding best fit parameters of quark Sivers functions. The GSF parameterizations denoted by ”SIDIS1” and ”SIDIS2” are given in refs Anselmino et al. (2005, 2009). MSTW2008 Martin et al. (2009) is used for collinear PDFs probed at the fixed scale $\rm\mu=\sqrt{m_{1}^{2}+p_{1T}^{2}}$, and $\rm m_{1}=3.096$ GeV is the mass of $\rm J/\psi$. We use the wave function at origin and the LDME $\rm\langle O_{8}(^{3}S_{1})\rangle$ whose values are listed as Li et al. (2009); Ko et al. (2011)

[TABLE]

The transverse momentum of the two initial gluons $\rm k_{\perp 1}$, $\rm k_{\perp 2}$ in Eq.(3) are integrated within the limits $\rm 0<k_{\perp 1},k_{\perp 2}<1$ GeV. We have noticed that the higher values of $\rm k_{\perp 1max}$ and $\rm k_{\perp 2max}$ do not make a difference with the SSA and the unpolarized differential cross section.

Our predictions of SSA are given for three different CM energies $\rm\sqrt{s}$ = 115 GeV (AFTER@LHC), $\rm\sqrt{s}$ = 200 GeV (RHIC1) and $\rm\sqrt{s}$ = 500 GeV (RHIC2). As mentioned above, we impose the following overall cuts for three different CM energies so as to ignore $\rm q\bar{q}$ channel contributions to $\rm J/\psi$ pair, and apply TMD factorization:

[TABLE]

Moreover, we impose cuts on both individual $\rm J/\psi$’s rapidity $\rm-2.8<y_{1},y_{2}<0.2$ for $\sqrt{s}=115$ GeV, $\rm 2<y_{1},y_{2}<3$ and $\rm 3<y_{1},y_{2}<3.8$ for $\rm\sqrt{s}=200$ GeV, and $\rm 3<y_{1},y_{2}<4$ for $\rm\sqrt{s}=500$ GeV. The given rapidity ranges were chosen keeping in mind the proposed forward sPHENIX (fsPHENIX) upgrade Barish (2012); Aschenauer et al. (2015). In practice, we abandon the cut $\rm 3<y_{1},y_{2}<3.8$ for $\sqrt{s}=200$ GeV, since it cuts off almost all events. We present the asymmetry predictions as a function of the sum transverse momentum ($\rm q_{T}$) and the rapidity ($\rm y$) of two quarkonium system as well as the rapidity difference of two $\rm J/\psi$ ($\rm\Delta y$).

The results are shown in FIG.1-4 and the conventions in these figures are the following. The ”SIDIS1” and ”SIDIS2” curves are obtained by using D’Alesio at el. D’Alesio et al. (2015) fit parameters of GSF. The obtained asymmetry using Anselmino et al. Anselmino et al. (2017a) fit parameters is represented by ”BV-a” and ”BV-b”. As aforementioned, we have considered the final heavy quarks produced to be in the both CS and CO states for calculating the numerator and denominator part of Eq.(1). We find that the SSA with also color octet state contribution is very close to the one without color octet state contribution.

From FIG.1-4, the asymmetries are estimated to be positive and negative with SIDIS and BV parameterization respectively, as functions of $\rm q_{T}$, $\rm y$ and $\rm\Delta y$. The estimated asymmetry using ”SIDIS2” fit is close to zero for all $\rm\sqrt{s}$ and the physical quantities $\rm q_{T}$, $\rm y$ and $\rm\Delta y$, considering our adopted cuts. While the estimated $\rm|A_{N}|$ using ”BV-b” fit is larger than all the other results in all CM energies and all considered physical quantities. In addition, the estimated asymmetries using all four parameters vary a little in the considered ranges of $\rm q_{T}$, $\rm y$ and $\rm\Delta y$. The obtained asymmetry using ”BV-b” parameters is maximum about 12% as functions of $\rm q_{T}$ and $\rm y$ at AFTER@LHC $\rm\sqrt{s}$ in FIG.1 (left panel) and FIG.1 (middle panel). It is also shown that the estimated asymmetry using ”SIDIS1” fit is very close to the asymmetry using ”SIDIS2” fit in AFTER@LHC and RHIC2 ($\rm 2<y_{1},y_{2}<3$). While in RHIC1 and RHIC2 ($\rm 3<y_{1},y_{2}<4$) the ”SIDIS1” fit results are larger than those from ”SIDIS2”. The sign of the asymmetry for ”BV-a” and ”BV-b” fit parameters depends on the magnitude of $\rm N_{d}$ and the relative magnitude of $\rm N_{u}$ and $\rm N_{d}$, and these have opposite sign which can be observed in TABLE.1. The negative sign magnitude of $\rm N_{d}$ is larger compared to $\rm N_{u}$ as a result the ”BV-b” asymmetry is negative. That is to say, the modeling of GSF strongly decides the magnitude and sign of the asymmetry.

In FIG.5 and Fig.6, we present the asymmetry predictions obtained with the two GSF fits, BV-b and SIDIS1, for all the three CM energies considered. It should be noted that in FIG.5 and 6, the y distribution peaks are in negative region for AFTER@LHC CM energy. This is due to the fact that AFTER@LHC is a fixed target experiment and we have taken $\rm y_{cm}$ to be positive in the unpolarized beam direction. In FIG.5 and Fig.6, the y-distributions of the SSA lie from left to right as CM energy increases from 115 to 500 GeV, duing to the responding rapidity cuts. This is in contrast to RHIC1 and RHIC2 curves, where we have used the convention followed by PHENIX experiment, in which rapidity is considered to be positive in the forward hemisphere of the polarized proton. We find the largest asymmetry values are of about 20% for the $\rm q_{T}$-distribution with BV-b fit and 11% for the y-distribution with SIDIS1 fit. For the $\rm q_{T}$-distributions with BV-b fit, the asymmetry becomes larger as the value of $\rm q_{T}$ increases. Generally, the asymmetry increases when the CM energy decreases. Whereas, it is desired to notice in FIG.5 and 6 that the asymmetry for 115 GeV has the lowest value than all the others in $\rm q_{T}$, $\rm y$ and $\rm\Delta y$-distributions. This is because the energy increases and more strict rapidity cuts are applied, the unpolarized cross section decreases more quickly than the polarized one. Specially, one can find from FIG.6 that 200 GeV (RHIC1) case gives the largest asymmetries in all three distributions. Moreover, in the case of $\rm q_{T}$-asymmetries, shown in FIG.5 (left panel) and 6 (left panel), we find that the functional form of the $\rm q_{T}$ dependence remains the same up to an overall factor that depends on $\rm\sqrt{s}$ and the rapidity range. This is also a reflection of the factorized $\rm q_{T}$ dependence that we have assumed for the TMDPDFs.

IV SUMMARY AND DISCUSSIONS

In this paper, we have calculated the single spin asymmetry (SSA) in the double $\rm J/\psi$ production. Within the NRQCD framework, the color singlet state $\rm{}^{3}S_{1}^{(1)}$ and color octet state $\rm{}^{3}S_{1}^{(8)}$ contributions to the double $\rm J/\psi$ production are considered. We have given out the SSA as functions of $\rm p_{T}$, $\rm y$ and $\rm\Delta y$. Sizable asymmetry is obtained as functions of them in the kinematic range $\rm\hat{s}=M^{2}\gg\Lambda_{QCD}$, $\rm 0<q_{T}<1\ \text{GeV}$, $\rm 4<p_{2T}<6\ \text{GeV}$. Typically we have considered $\rm-2.8<y_{1},y_{2}<0.2$ for $\rm\sqrt{s}=115$ GeV, $\rm 2<y_{1},y_{2}<3$ and $\rm 3<y_{1},y_{2}<3.8$ for $\rm\sqrt{s}=200$ GeV, and $\rm 3<y_{1},y_{2}<4$ for $\rm\sqrt{s}=500$ GeV. We find the asymmetry can reach more than 10 percent using BV-b fit. The results using the other three fit parameters are also not small. The sizable asymmetry indicates that the double $\rm J/\psi$ production in proton proton collision is a competitive process to probe the gluon Sivers function over a wide kinematic region accessible at the RHIC and AFTER@LHC.

Acknowledgements.

Xuan Luo thanks professor Sergey Baranov and professor Asmita Mukherjee for very useful discussions. Hao Sun is supported by the National Natural Science Foundation of China (Grant No.11675033), and by the Fundamental Research Funds for the Central Universities (Grant No. DUT18LK27).

Appendix A Square of the amplitude for $\rm g+g\rightarrow J/\psi+J/\psi$ process

The amplitude squares of gluon-gluon to $\rm J/\psi$ pair can be calculated by the FORM package Kuipers et al. (2013) straightforwardly and we have cross checked with paper Li et al. (2009)Qiao (2002). The amplitude squares of $\rm{}^{3}S_{1}^{(1)}$, $\rm{}^{3}S_{1}^{(1)}$ and $\rm{}^{3}S_{1}^{(8)}$, $\rm{}^{3}S_{1}^{(8)}$ are given below

[TABLE]

[TABLE]

where m is the mass of $\rm J/\psi$; $\rm|R(0)|$ is the magnitude of its radial wave function at origin; $\rm\hat{s}$, $\rm\hat{t}$, $\rm\hat{u}$ are the usual Mandelstam variables.