Pinning down the linearly-polarised gluons inside unpolarised protons using quarkonium-pair production at the LHC

We show that the production of J/psi or Upsilon pairs in unpolarised pp collisions is currently the best process to measure the momentum-distribution of linearly-polarised gluons inside unpolarised protons through the study of azimuthal asymmetries. Not only the short-distance coefficients for such reactions induce the largest possible cos 4 phi modulations, but analysed data are already available. Among the various final states previously studied in unpolarised pp collisions within the TMD factorisation approach, di-J/psi production exhibits by far the largest asymmetries, up to 50% in the region studied by the ATLAS and CMS experiments. In addition, we use the very recent LHCb data at 13 TeV to perform the first fit of the unpolarised transverse-momentum-dependent gluon distribution.

Introduction.-Probably one of the most striking phenomena arising from the extension of the collinear factorization -inspired from Feynman's and Bjorken's parton model-to Transverse Momentum Dependent (TMD) factorization [1][2][3][4] is the appearance of azimuthal modulations induced by the polarization of partons with nonzero transverse momentumeven inside unpolarized hadrons. In the case of gluons in a proton, which trigger most of the scatterings at high energies, this new dynamics is encoded in the distribution h ⊥ g 1 (x, k 2 T ) of linearly-polarized gluons [5]. In practice, they generate cos 2φ (cos 4φ) modulations in gluon-fusion scatterings where single (double) gluon-helicity flips occur. They can also alter Transverse Momentum (TM) spectra, such as that of a Brout-Englert-Higgs H 0 boson [6,7], via double gluon-helicity flips. Beside containing crucial information to better understand the dynamics of gluons confined inside protons, their knowledge is also important to advance the precision of kinematicaldistribution predictions at the LHC.
In this Letter, we show that di-J/ψ production, which among the quarkonium-associated-production processes has been the object of the largest number of experimental studies at the LHC and the Tevatron [8][9][10][11][12], is in fact the ideal process to perform the first measurement of h ⊥ g 1 (x, k 2 T ). It indeed exhibits the largest possible azimuthal asymmetries in regions already accessed by the ATLAS and CMS experiments where such modulations can be measured. Along the way of our study, we perform the first extraction of f g 1 (x, k 2 T ) -its unpolarized counterpart-using recent LHCb data.
Its structure is as follows. We first recall the necessary concepts of TMD factorization in gluon-induced scatterings of unpolarized pp collisions. After discussing the applicability of TMD factorization to vector-quarkonium-pair (QQ) production, we outline the main formulae relevant for this case. This prepares our discussion of the extraction of f g 1 and of our predictions of the azimuthal asymmetries induced by the linearly-polarized gluons.
TMD factorization for gluon-induced scatterings.-TMD factorization extends collinear factorization by accounting for the parton TM, generally denoted by k T . It is particularly relevant in scatterings with large momentum transfers where the effect of k T -on the order of the proton mass-are nevertheless visible. This can usually be achieved when a pair of particles is observed with a large invariant mass but with a small TM. Di-Q production at the LHC falls in this category.
In practice, the gluon TMDs in an unpolarized proton of momentum P and mass M p are defined through the hadron correlator Φ µν g (x, k T ) [5,13], parametrized in terms of two independent TMDs, the unpolarized distribution f g 1 (x, k 2 T ) and the distribution of linearly-polarized gluons h ⊥ g 1 (x, k 2 T ), as where the gluon four-momentum k is decomposed as k = xP+ k T + k − n [n is any light-like vector (n 2 = 0) such that n · P 0], k 2 T = −k 2 T and g µν T = g µν − (P µ n ν + P ν n µ )/P·n. Under TMD factorization (see e.g. Fig. 1) and up to corrections suppressed by powers of the observed system TM over its invariant mass, the cross section for any gluon-fusion process (here g(k 1 ) + g(k 2 ) → Q(P Q,1 ) + Q(P Q,2 ) ) can be expressed as the convolution of a partonic short-distance contri- where s = (P 1 + P 2 ) 2 , x {1,2} = P QQ ·P {2,1} /P 1 ·P 2 with P QQ = P Q,1 + P Q,2 and dR is the phase space element of the outgoing particles. The partonic hard-scattering amplitude M for the gluon-initiated subprocess can be calculated in perturbative QCD through a series expansion in α s [14].
Owing to process-dependent Wilson lines in the definition of the correlators which they parametrize, the TMDs are in general not universal. Physics wise, these Wilson lines describe the non-perturbative interactions of the active parton -the gluon in our case-with soft spectator quarks and gluons in the nucleon before or after the hard scattering. For the production of di-leptons, γγ, di-Q or boson-Q pairs via a Color-Singlet Transitions (CST) [15][16][17] i.e. for purely colorless final states-in pp collisions, only initial-state interactions (ISI) between the active gluons and the spectators can occur. Mathematically, these ISI can be encapsulated [18] in TMDs with past-pointing Wilson lines -the exchange can only occur before the hard scattering. Such gluon TMDs correspond to the Weiszäcker-Williams distributions relevant for the low-x region [19,20]. Besides, in lepton-induced production of colorful final states, like heavy-quark pair or dijet production [21,22], to be studied at a future Electron-Ion Collider (EIC) [23], only final-state interactions (FSI) take place.
Yet, since f g 1 and h ⊥ g 1 are time-reversal symmetric (T -even), contrary to other TMDs [24,25] like the gluon distribution in a transversally polarised proton also called the Sivers function [26], TMD factorization tells us that one in fact probes the same distributions in both the production of colorless systems in hadroproduction and of colorful systems in leptoproduction. Extracting them in different reactions is thus essential to test this universality property of the TMDs -akin to the well-known sign change of the quark Sivers effect [18,27]-, in order to validate TMD factorization.
Di-Q production & TMD factorization.-For TMD factorization to apply, di-Q production should at least satisfy both following conditions. First, it should result from a Single-Parton Scattering (SPS). Second, FSI should be negligible, which is satisfied when quarkonia are produced via CSTs [14].
Double-parton-scatterings (DPSs) leading to di-J/ψ are not problematic as they only become significant at large ∆y (or large pair invariant masses) -as noted for other observables [28][29][30][31][32][33][34]. For small ∆y, in particular in the D0, CMS and ATLAS samples where a TM cut was applied on the J/ψ, the SPSs are dominant. As for the dominance of CSTs to the SPS yield, we guide the reader to [35][36][37]. The latter is only challenged in regions where the DPS mechanism is anyhow dominant (large ∆y) which we will not consider. We recall that, within Non-Relativistic QCD (NRQCD) [38] (see [39][40][41] for reviews), the CSTs are leading-order (LO) in the v 2 expansion-v being the heavy-quark velocity in the Q rest frame-and have been studied up to next-to-leading (NLO) accuracy in α s [42][43][44] in collinear factorization. The feed down from excited states is also not problematic: J/ψ+χ c production is suppressed [35] and J/ψ + ψ can be treated exactly like J/ψ + J/ψ. As for di-Υ, the CST should be even more dominant and we could not find any argument why the DPS/SPS ratio could be larger. To conclude, we stress that the current study is in fact the first one in TMD factorization for vector QQ.
Using the notations of [45], the structure of the TMD cross section for QQ production reads where dΩ = dcos θ CS dφ CS is expressed in terms of Collins-Soper (CS) angles [46] and where M QQ , Y QQ and P QQT are the invariant mass, the rapidity and the TM of the pair -the latter two to be measured in the hadron c.m.s.. In the CS frame, the Q direction is along e = (sin θ CS cos φ CS , sin θ CS sin φ CS , cos θ CS ). Let us stress that the overall factor is a phase space factor specific to the mass of the final-state particles and the analyzed differential cross sections, and that the hard factors F i do not depend neither on Y QQ nor P QQT . In addition, let us note that -away from thresholdcos θ CS ∼ 0 corresponds to ∆y ∼ 0 in the hadron c.m.s., that is our preferred region to avoid DPS contributions. The TMD convolutions in Eq. (3) are defined as The explicit expressions for the weights in Eq. (3) are identical for all the gluon-induced processes and can be found in [45].
The short-distance coefficients F i -The factors F i are calculable process by process and we refer to [45] for details on how to obtain them from specific combinations of the helicity amplitudes of the partonic scattering under consideration. As such, they can be derived from the uncontracted amplitude given in [47]. Their full expressions, which are of limited interest for the following discussions, can be found as supplemental material (SM)[48] along with a short derivation that, for any process, .] Yet, both the limits of large and small QQ mass, M QQ , are very interesting. Indeed, when M QQ is much larger than the quarkonium mass, M Q , one finds that, for cos θ CS = c θ → 0, where N = 2 11 3 −4 π 2 α 4 s |R Q (0)| 4 , R Q (0) being the Q radial wave function at the origin. As in collinear factorization, the Bornorder cross section scales as α 4 s . Computing the one-loop QCD corrections would introduce the same α s corrections to all the F i (see [49] for the Drell-Yan case). At threshold, M QQ → 2M Q , one gets: Beside the cos θ CS dependence, which indicates that F 2,3 are suppressed near ∆y ∼ 0, one also observes that F 2 (F 3 ) scales like M −4 QQ (M −2 QQ ) relative to F 1 and F 4 . In other words, the modification of the P QQT dependence due to the linearlypolarized gluons encoded in F 2 vanishes at large scales [From Eq. (6), one notes that it is however already very small even at threshold.] whereas the cos 4φ CS modulation (double helicity flip) takes over the cos 2φ CS one (single helicity flip).
The fact that F 4 → F 1 , for c θ → 0 away from the threshold, is the most important result of this study and is, to the best of our knowledge, a unique feature of di-J/ψ and di-Υ production. As such, and thanks to the possibility to re-analyze collected di-J/ψ data, it is indeed the ideal one to extract the linearly-polarized gluon distributions. The previously studied γγ [50], H 0 +jet [51], Q+γ [52], Q+γ or Q+Z [45] processes show significantly smaller values of F 4 /F 1 , thus correspondingly smaller values of the cos 4φ CS modulation for a given magnitude of h ⊥g 1 . Knowing the F i and an observed differential yield, one can thus extract the various TMD convolutions of Eq. (4) from their azimuthal (in)dependent parts. When the cross section is integrated over φ CS , the contribution from F 3,4 drops out from Eq. (3) and only depends on C f g 1 f g 1 and C w 2 h ⊥g 1 h ⊥g 1 . To go further, we define cos nφ CS [for n = 2, 4] weighted differential cross sections normalized to the azimuthally independent term as: It is understood that cos nφ CS computed in a range of M QQ , Y QQ , P QQT or cos θ CS is the ratio of corresponding integrals. Using Eq. (3), one gets in a single phase-space point: The TM spectrum.-Before discussing the expected size of the azimuthal asymmetries for di-Q production, let us have a closer look at the TM dependence of Eq. (3), which is entirely encoded in the four convolutions defined in Eq. (4). These are process independent quantities, unlike the F i . Since gluon TMDs are still unknown, we need to resort to models in order to predict their impact.
As proposed in e.g. [53], one can assume a simple Gaussian dependence on k 2 T for f g 1 , namely f g is the common collinear gluon distribution. Concerning h ⊥g 1 , we know that it is constrained by the model-independent positivity bound [5]: holding for any value of x and k 2 T . As noted in [6], this bound is satisfied by the following form h ⊥g with r < 1 and we take r = 2/3 which maximizes the second moment in k T of h ⊥g 1 . With such a choice, all the TMD convolutions can analytically be computed. They are gathered as supplemental material. It is nevertheless also customary to build up a model of the TMD h ⊥g 1 by saturating its bound, namely h ⊥g . Such a choice is for instance motivated by recent studies in the highenergy (low-x) limit (see e.g. [19,54]) where this bound is reached. The corresponding convolutions can easily be calculated numerically and their behavior are also gathered as supplemental material for the interested reader. Having two models at hand is very convenient, not only to obtain a range of the modulations, but also to evaluate how efficient a given measurement is to constrain the TMD h ⊥g 1 by comparing its uncertainty to our range. In the following, "Model 1" refers to the Gaussian form with r = 2/3 and "Model 2" to the form saturating the positivity bound.
Apart from the functional form, the TMD convolutions, C[w f g] -from which follows the P QQT spectra-only depend on k 2 T . For the first time, we are able to use experimental data to fix it via C[ f 1 f 1 ] from the P QQT spectrum recently measured by the LHCb Collaboration at 13 TeV [12] (see Fig. 2a). The LHCb study was made without any TM cuts, thus near threshold where M QQ ∼ 2M Q . As such, we can safely neglect any contribution from h ⊥g 1 given the size of F 2 /F 1 near threshold (cf. Eq. (6)). We further note that, for TMD Ansätze with factorized dependences on x and k 2 T , the normalized P QQT spectrum does not depend on any other kinematical variables. By fitting the P QQT spectrum up to M QQ /2, we obtain k 2 T = 4.9±0.8 GeV 2 . To the best of our knowledge, this is the first time that such a process-independent quantity is experimentally determined in a pure gluon-induced process with a colorless final state, for which TMD factorization applies. The fact that the TMD curve undershoots the data for P QQT M QQ /2 is expected, as it leaves room for hard finalstate radiations not accounted for in TMD factorization as well as DPS contributions exactly where they can appear.
In this fit, and for the rest of study, we have neglected the scale evolution effects of the TMDs. These are believed to increase k 2 T [6,55]. We have checked that increasing k 2 T as log(M 2 QQ ) [M QQ being the hard scale] induces changes in cos nφ CS which systematically remains smaller than switching from Model 1 and 2. As such, we are confident that performing a full study with evolution would not change our conclusions relative to the -large-size of cos nφ CS for di-J/ψ production and its uniqueness to experimentally determine h ⊥g 1 . On the other hand, a reliable global analysis of TM spectra of H 0 and di-J/ψ production in the LHCb, CMS and ATLAS acceptance should account for these evolution effects owing to the large scale differences between such various data samples. This is left for a future study.
In our above extraction of k 2 T , we have neglected the influence of h ⊥g 1 on the P QQT spectrum. We note that the situation is analogous to Q + γ [52], Q + γ or Q + Z [45] with a negligible impact of h ⊥g 1 on the TM spectra but significantly different from that for di-photon [50], single η c [56], di-η c [57] and H 0 +jet [51] production. Data nonetheless do not exist yet for any of these channels. Unfortunately, the CMS di-Υ sam-  2. (a) The normalized P QQT dependence of the di-J/ψ yield obtained with a Gaussian f g 1 with k 2 T fit to the normalized LHCb data at 13 TeV [12] [The data in the gray zone were not used for the fit since TMD factorization does not apply there]. (b,c) cos nφ CS for n = 2, 4 computed for | cos θ CS | < 0.25 for both our model of h g 1⊥ , for 3 values of M QQ (8, 12 and 21 GeV) relevant respectively for the LHCb [12], CMS [10] and ATLAS [11] kinematics. The spectra are plotted up to M QQ /2. Our results do not depend on Y QQ . ple [58] is not large enough (40 events) to check our k 2 T fit. With 100 fb −1 of 13 TeV data, their sample should allow for such a check.
Azimuthal dependences.-Having fixed the functional form of the TMDs and k 2 T and having computed the factors F i , we are now ready to provide predictions for the azimuthal modulations through cos nφ CS as a function of P QQT , cos θ CS or M QQ . Fig. 2b & 2c show cos nφ CS (n = 2, 4) as a function of P QQT for both our models of h ⊥ g 1 for 3 values of M QQ , 8, 12 and 21 GeV for | cos θ CS | < 0.25. These values are relevant respectively for the LHCb [12], CMS [10] and ATLAS [59] kinematics. Still to keep TMD factorization applicable, we have plotted the spectra up to M QQ /2 . Let us also note that with our factorized TMD Ansätze, cos nφ CS do not depend on Y QQ . Indeed, the pair rapidity only enters the evaluation of dσ via the momentum fractions x 1,2 in the TMDs. It thus simplifies in the ratios.
Compared to previously studied cases, the size of the expected azimuthal asymmetries is particularly large, e.g. for P 2 QQT k 2 T . cos 4φ CS even gets close to 50 % in the P QQT region probed by CMS and ATLAS for | cos θ CS | < 0.25; this is probably the highest value ever predicted for a gluonfusion process which directly follows from the extremely favorable hard coefficient F 4 -as large as F 1 . Such values are truly promising to extract the distribution h ⊥g 1 of linearlypolarized gluons in the proton which appears quadratically in cos 4φ CS . In view of these results, it becomes clear that the kinematics of CMS and ATLAS, with naturally larger invariant masses, are better suited with much larger expected asymmetries than that of LHCb, not far from threshold. cos 2φ CS , which would allow one to lift the sign degeneracy of h ⊥g 1 in cos 4φ CS , is on the order of a few per cent for | cos θ CS | < 0.25 (Fig. 2b). Yet, we have seen that F 3 vanishes for small cos θ CS (Eq. (5)). It would thus be expedient to extend the range of | cos θ CS | pending the DPS contamination. Indeed, in view of recent di-J/ψ phenomenological studies [35,60,61], one expects the DPSs to become dominant at large ∆y while these cannot be treated along the lines of our analysis. To ensure the SPS dominance, it is thus judicious to avoid the region ∆y > 2, and probably ∆y > 1 to be on the safe side. Even though the relation between ∆y -measured in the hadronic c.m.s.-and cos θ CS is in general not trivial, it strongly simplifies when P 2 QT (M 2 Q , P 2 QQT ), such that cos θ CS = tanh ∆y/2 [62]. Up to | cos θ CS | ∼ 0.5, the sample should thus remain SPS dominated in particular with the CMS and ATLAS P QT cuts. In fact, in a bin 0.25 < | cos θ CS | < 0.5, cos 2φ CS nearly reaches 30% . On the contrary, cos 4φ CS exhibits a node close to cos θ CS ∼ 0.3.
[See plots as supplemental material]. As such, measuring cos 4φ CS for | cos θ CS | < 0.25 and 0.25 < | cos θ CS | < 0.5 would certainly be instructive. If our models for h ⊥g 1 are realistic, this is definitely within the reach of CMS and ATLAS, probably even with data already on tape.
Conclusions.-We have found out that the short-distance coefficients to the azimuthal modulations of di-J/ψ and di-Υ yields equate the azimuthally independent terms, which renders these processes ideal probes of the linearly-polarized gluon distributions in an unpolarized proton, h ⊥ g 1 . This is even more true that experimental data already exist -more will be recorded in the near future-and it only remains to analyze them along the lines discussed above, by evaluating the ratios cos 2φ CS and cos 4φ CS . In fact, we have already demonstrated the relevance of the LHC data for di-J/ψ production by constraining, for the first time, the TM dependence of f g 1 at a scale close to 2M ψ .
Let us also note that similar studies can be carried out at fixed-target set-ups where luminosities are large enough to detect J/ψ pairs. The COMPASS experiment with pion beams may also record di-J/ψ events as did NA3 in the 80's [63,64]. Whereas single-J/ψ production may partly be from quarkantiquark annihilation, di-J/ψ production should mostly be from gluon fusion and thus analyzable along the above discussions. Using the 7 TeV LHC beams [65] in the fixed-target mode with a LHCb-like detector [66][67][68], one can expect 1000 events per 10 fb −1 , enough to study a possible x dependence of k 2 T as well as to look for azimuthal asymmetries generated by h g 1⊥ . Such studies could also be complemented with targetspin asymmetry studies [69][70][71], to study the gluon Sivers function f ⊥g 1T as well as the gluon transversity h g 1T or the distribution of linearly-polarized gluons in a transversely polarized proton, h ⊥g 1T paving the way for an in-depth gluon tomography of the proton.
Acknowledgements.-We thank D. Boer, M. Echevarria and H.S. Shao for useful comments and L.P. Sun for discussions about [47]. The work of J.P.L. and F.S. is supported in part by the French IN2P3-CNRS via the LIA FCPPL (Quarkonium4AFTER) and the project TMD@NLO.
The work of C.P. is supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 647981, 3DSPIN). The work of M.S. is supported in part by the Bundesministerium für Bildung und Forschung (BMBF) grant 05P15VTCA1.

SUPPLEMENTAL MATERIAL
A. Bounds on the short-distance coefficients Instead of using their definition in terms of the helicity amplitudes [45], one can express the short-distance coefficients F i using the different combinations of the partonic hard-scattering amplitude elements in the CS frame. One then gets that they only depend on the x and y elements of M : It is then straightforward from Eq. (1) to infer the following bound : F ( ) 2,3,4 ≤ F 1 . These expressions are valid for any gluon fusion process computed within the framework of TMD factorization.
B. Complete expressions of the short-distance coefficients F i As for the F i we have