Associated production of a dilepton and a $\Upsilon(J/\psi)$ at the LHC as a probe of gluon transverse momentum dependent distributions

TL;DR

This paper explores how the associated production of a dilepton and a quarkonium state at the LHC can serve as a probe for gluon transverse momentum dependent distributions, linking measurable cross sections to gluon TMDs.

Contribution

It provides a detailed theoretical framework connecting dilepton-quarkonium production to gluon TMDs, including a phenomenological analysis of TMD extraction feasibility.

Findings

Cross sections relate to unpolarised and linearly polarised gluon TMDs.

Feasibility analysis supports potential for experimental TMD extraction.

Framework applicable to LHC energies for gluon structure studies.

Abstract

We discuss the impact on the study of gluon transverse momentum dependent distributions (TMDs) of the associated production of a lepton pair and a $\Upsilon$ (or a $J/\psi$) in unpolarised proton-proton collisions, $pp\to {\cal Q} \, \ell \bar{\ell} X$, at LHC energies, where one can assume that such final states are dominantly induced by gluon fusion. If the transverse momentum of the quarkonium-dilepton system - namely, the transverse momentum imbalance of the quarkonium state and the lepton pair - is small, the corresponding cross sections can be calculated within the framework of TMD factorisation. Using the helicity formalism, we show in detail how these cross sections are connected to the moments of two independent TMDs: the distribution of unpolarised gluons, $f_1^g$, and the distribution of linearly polarised gluons, $h_1^{\perp g}$. We complete our exhaustive derivation of these general relations with a phenomenological analysis of the feasibility of the TMD extraction, as well as some outlooks.

Equations

$$\mathrm{d}{\sigma }_{\mathrm{TMD}}$$

$${\Gamma }^{\mu \nu }(x,{\mathbf{k}}_T)=\frac{1}{2x}\left(-g_T^{\mu \nu }{f}_1^g(x,{\mathbf{k}}_T^2)+\frac{k_T^{\mu }{k}_T^{\nu }+\frac{1}{2}{\mathbf{k}}_T^2{g}_T^{\mu \nu }}{M^2}{h}_1^{\perp g}(x,{\mathbf{k}}_T^2)\right).$$

$${\epsilon }_{{\lambda }_a}^{\mu }({\bar{k}}_a)=\left(0,-\frac{{\lambda }_a}{\sqrt{2}},-\frac{i}{\sqrt{2}},0\right)$$

$$-g_T^{\mu \nu }=\sum _{{\lambda }_{a/b}}{\epsilon }_{{\lambda }_{a/b}}^{\mu }({\bar{k}}_{a/b}){\epsilon }_{{\lambda }_{a/b}}^{\nu \ast }({\bar{k}}_{a/b}).$$

$$\mathrm{d}{\sigma }_{\mathrm{TMD}}=\frac{(2\pi {)}^4}{S^2}{\mathrm{dPS}}_n\sum _{{\lambda }_a,{\bar{\lambda }}_a,{\lambda }_b,{\bar{\lambda }}_b=\pm 1}\frac{1}{(N_c^2-1{)}^2}\sum _{a,b;I}{\mathcal{A}}_{{\lambda }_a,{\lambda }_b;I}^{ab}({\bar{k}}_a,{\bar{k}}_b;\{P_i\}){\mathcal{A}}_{{\bar{\lambda }}_a,{\bar{\lambda }}_b;I}^{ab\ast }({\bar{k}}_a,{\bar{k}}_b;\{P_i\})\times$$

$$\int {\mathrm{d}}^2{\mathbf{k}}_{aT}\int {\mathrm{d}}^2{\mathbf{k}}_{bT}{\delta }^{(2)}({\mathbf{k}}_{aT}+{\mathbf{k}}_{bT}-{\mathbf{q}}_T){\Gamma }_{{\bar{\lambda }}_a,{\lambda }_a}(x_a,{\mathbf{k}}_{aT},{\zeta }_a,\mu ){\Gamma }_{{\bar{\lambda }}_b,{\lambda }_b}(x_b,{\mathbf{k}}_{bT},{\zeta }_b,\mu )+\mathcal{O}(q_T/Q).$$

$${\Gamma }_{{\bar{\lambda }}_a,{\lambda }_a}(x_a,{\mathbf{k}}_{aT})$$

$${\Gamma }_{{\bar{\lambda }}_b,{\lambda }_b}(x_b,{\mathbf{k}}_{bT})$$

$$\mathrm{d}{\sigma }_{\mathrm{TMD}}$$

$$\mathcal{C}[wfg]\equiv \int {\mathrm{d}}^2{\mathbf{k}}_{aT}\int {\mathrm{d}}^2{\mathbf{k}}_{bT}{\delta }^{(2)}({\mathbf{k}}_{aT}+{\mathbf{k}}_{bT}-{\mathbf{q}}_T)w({\mathbf{k}}_{aT},{\mathbf{k}}_{bT},{\mathbf{q}}_T)f(x_a,{\mathbf{k}}_{aT}^2)g(x_b,{\mathbf{k}}_{bT}^2).$$

$$w_2=\frac{2({\mathbf{k}}_{aT}\cdot {\mathbf{k}}_{bT}{)}^2-{\mathbf{k}}_{aT}^2{\mathbf{k}}_{bT}^2}{4M^4}$$

$$w_{3a}=\frac{{\mathbf{k}}_{aT}^2{\mathbf{q}}_T^2-2({\mathbf{q}}_T\cdot {\mathbf{k}}_{aT}{)}^2}{2M^2{\mathbf{q}}_T^2}$$

$${\hat{F}}_1({\bar{k}}_a,{\bar{k}}_b;\{P_i\})$$

$${\hat{F}}_2({\bar{k}}_a,{\bar{k}}_b;\{P_i\})$$

$${\hat{F}}_{3a}({\bar{k}}_a,{\bar{k}}_b;\{P_i\})$$

$${\hat{F}}_{3b}({\bar{k}}_a,{\bar{k}}_b;\{P_i\})$$

$${\hat{F}}_4({\bar{k}}_a,{\bar{k}}_b;\{P_i\})$$

$${\mathcal{A}}_{{\lambda }_a,{\lambda }_b}^{gg\to \mathcal{Q}\ell \bar{\ell }}=\sum _{j={\gamma }^{\ast },Z}i{C}_j{\mathcal{A}}_{{\lambda }_a,{\lambda }_b}^{gg\to \mathcal{Q}j;\mu }\frac{g_{\mu \nu }-P_{B\mu }{P}_{B\nu }/M_B^2}{M_B^2-{\Delta }_j}\bar{u}(l,\sigma ){\gamma }^{\nu }(a_j+b_j{\gamma }_5)v(\bar{l},\bar{\sigma }).$$

$${\mathcal{F}}_{{\lambda }_a{\lambda }_b;{\bar{\lambda }}_a{\bar{\lambda }}_b}$$

$$L_{j\bar{j}}^{\nu \sigma }=4(a_j{a}_{\bar{j}}+b_j{b}_{\bar{j}})\left(l^{\nu }{\bar{l}}^{\sigma }+l^{\sigma }{\bar{l}}^{\nu }-\frac{1}{2}{M}_B^2{g}^{\nu \sigma }\right)+4(a_j{b}_{\bar{j}}+b_j{a}_{\bar{j}}){ϵ}^{\nu \sigma \eta \kappa }{l}_{\eta }{\bar{l}}_{\kappa }.$$

$$\int \mathrm{d}\Omega L_{j\bar{j}}^{\nu \sigma }$$

$${\mathcal{F}}_{{\lambda }_a{\lambda }_b;{\bar{\lambda }}_a{\bar{\lambda }}_b}$$

$$\frac{\mathrm{d}{\sigma }_{\mathrm{TMD},\mathrm{LO}}^{pp\to J/\psi [{\rm Y}]\ell \bar{\ell }X}}{{\mathrm{d}}^{4}q\mathrm{d}{M}_B^2\mathrm{d}\Omega }$$

$${\hat{F}}_i=\frac{4{\alpha }_s^2|R_0(0){|}^2\Lambda M_B^2}{3{\pi }^3{x}_a{x}_b{S}^2{Q}^3{M}_{\mathcal{Q}}{N}_c(N_c^2-1{)}^2}\sum _{j,\bar{j}={\gamma }^{\ast },Z}\frac{{\tilde{C}}_j{\tilde{C}}_{\bar{j}}(a_j{a}_{\bar{j}}+b_j{b}_{\bar{j}})}{(M_B^2-{\Delta }_j)(M_B^2-{\Delta }_{\bar{j}}{)}^{\ast }}\frac{{\hat{f}}_i(\alpha ,\beta ,\theta )}{\hat{D}(\alpha ,\beta ,\theta )},$$

$$\hat{D}(\alpha ,\beta ,\theta )$$

$${\hat{f}}_1(\alpha ,\beta ,\theta )$$

$${\hat{f}}_2(\alpha ,\beta ,\theta )$$

$${\hat{f}}_{3a}(\alpha ,\beta ,\theta )$$

$${\hat{f}}_4(\alpha ,\beta ,\theta )$$

$$a_{1,0}$$

$$a_{1,2}$$

$$a_{1,4}$$

$$a_{2,0}$$

$$a_{3,0}$$

$$a_{4,0}$$

$$N^{(0)}\equiv \int \mathrm{d}\Omega \frac{\mathrm{d}{\sigma }_{\mathrm{TMD},\mathrm{LO}}^{pp\to J/\psi [{\rm Y}]\ell \bar{\ell }X}}{{\mathrm{d}}^{4}q\mathrm{d}{M}_B^2\mathrm{d}\Omega }=\frac{\mathrm{d}{\sigma }_{\mathrm{TMD},\mathrm{LO}}^{pp\to J/\psi [{\rm Y}]\ell \bar{\ell }X}}{{\mathrm{d}}^{4}q\mathrm{d}{M}_B^2}$$

$$N^{(2)}\equiv \int \mathrm{d}\Omega \cos 2\varphi \frac{\mathrm{d}{\sigma }_{\mathrm{TMD},\mathrm{LO}}^{pp\to J/\psi [{\rm Y}]\ell \bar{\ell }X}}{{\mathrm{d}}^{4}q\mathrm{d}{M}_B^2\mathrm{d}\Omega }$$

$$N^{(4)}\equiv \int \mathrm{d}\Omega \cos 4\varphi \frac{\mathrm{d}{\sigma }_{\mathrm{TMD},\mathrm{LO}}^{pp\to J/\psi [{\rm Y}]\ell \bar{\ell }X}}{{\mathrm{d}}^{4}q\mathrm{d}{M}_B^2\mathrm{d}\Omega }$$

$${\hat{F}}_1(Q,\alpha ,\beta )$$

Full text

††journal: Nuclears Physics B

1 Introduction

Three-dimensional momentum distributions of gluons in the nucleon – the so-called gluon transverse momentum dependent distributions (TMDs) – have attracted much attention recently [1, 2]. Theoretically, gluon TMDs appear in factorisation formulae that explicitly take the transverse momentum of gluons into account (TMD factorisation), see e.g. [3, 4, 5, 6]). These formulae apply to cross sections that are differential in the transverse momentum of the final state, ${\bm{q}}_{T}$, in a kinematical region where it is much smaller than the hard scale of the process $Q$, $|{\bm{q}}_{T}|\ll Q$. Typically, the hard scale refers to the virtuality of the exchanged gauge boson in lepton-nucleon collisions or the invariant mass of the final state in hadron-hadron collisions. Our understanding of the evolution with the hard scale of the gluon TMDs has recently significantly improved, see e.g. [7, 8, 9, 10].

Several processes have been identified as sensitive probes of gluon TMDs. Probably, the theoretically cleanest one is the production of a pair of (almost) back-to-back heavy quark and antiquark or of a di-jet in lepton-nucleon collisions [11, 12, 13]. A measurement of such processes could however only be performed at a future Electron-Ion-Collider. On the other hand, reactions initiated by two protons can also provide insights on the gluon TMDs at the LHC or at RHIC. For instance, one possibility is to look at di-photon production, still with a large azimuthal separation [14]. This process, however, suffers from additional contributions from quark-induced channels, at RHIC energies in particular. As such, it is probably not the cleanest probe of gluon TMDs in hadron collisions that one could think of. Furthermore, the experimental detection of direct photons requires a specific isolation procedure which may be difficult to implement in a realistic measurement.

A handier probe of gluon TMDs is certainly to be found among the production of quarkonium states [15, 16, 17, 18, 19, 20, 21], built up of either charm or bottom quarks, since they are often dominantly produced through gluon fusion at proton-proton colliders and some of them, like the spin triplet vector states, are relatively easy to detect in their di-muons channels. As for now, the access to gluon TMDs has been investigated in [15, 6] with single inclusive $\eta_{c,b}$ or $\chi_{c,b}$-production in proton collisions. One drawback of such single-particle analyses is that they are restricted to low transverse momenta, typically below half the quarkonium mass, which makes such experimental studies particularly challenging. Furthermore, in the case of $\chi_{c,b}$-production [22], TMD factorisation may not hold because of infrared divergences specific to the $P$-wave production111For review, the readers is referred to  [26, 27, 28]. Such issues do not appear for spin singlet $S$-wave state production, which however has so far only been studied for the $\eta_{c}$ down to $P_{T}\simeq 6$ GeV by the LHCb collaboration [23].

This restriction on the usable phase space for TMD factorisation to apply can be avoided by looking at two-particle final states. One example is the associated production of a $J/\psi$ or a $\Upsilon$ with a direct photon [24]. Like for heavy-quark pair and di-jet electroproduction or di-photon hadroproduction, the large scale $Q$ is set by the invariant mass of the system, which can be large when both the quarkonium and the photon are produced almost back to back with large individual transverse momentum, yet with a small transverse momentum for the pair (its imbalance). This is a rather convenient configuration to be studied experimentally. Moreover, in the case of quarkonium + photon production [24], one can enrich the event sample in colour-singlet contributions, which are purely from gluonic interaction, by isolating the quarkonium, since it has a non-zero transverse momentum. This makes it a golden-plated probe to extract gluon TMDs inside unpolarized protons at the LHC. However, along the lines of [25], it may be that colour-octet contributions to quarkonium production associated with a colour singlet particle, like a SM boson ($\gamma$, $W^{\pm}$, $Z^{0}$ and $H^{0}$), could also be treated in the TMD factorisation framework. One of the crucial aspects yet to be fully understood in this case is whether an imbalance to be measured in the final state can be related to the transverse momentum of the initial partons. For this to be true, final-state gluon emissions should admittedly be suppressed at least to a tractable extent.

For unpolarised colliding protons at the LHC, there are two relevant gluon TMDs: the distribution of unpolarised gluons, $f_{1}^{g}$, and the distribution of linearly polarised gluons, $h_{1}^{\perp g}$. The latter is of particular interest as this distribution flips the helicity of the gluon entering the partonic cross section. As a result, the linear polarisation of the gluon manifests itself in two ways: a modification of the transverse-momentum dependence of the cross section in a characteristic way, and an azimuthal modulation of the cross section. As a matter of fact, the linear polarisation may serve as a general new tool in particle physics. Examples of the usefulness of the linear gluon polarisation have been discussed in the context of $H^{0}$ boson production in [29, 30], as well as for $H^{0}$+jet production [25].

Like all other TMDs, $f_{1}^{g}$ and $h_{1}^{\perp\,g}$ are affected by the presence of initial and final state interactions, whose effects are encoded in the Wilson lines needed for their gauge-invariant definition. TMD factorisation may therefore fail for some processes, as we already mentioned above. Moreover, TMDs may become process dependent even in those cases where factorisation can be proven. The gluon TMDs appearing in all the proton-proton scattering reactions where only initial state interactions are present, like the ones under study here, correspond to the so called Weiszäcker-Williams distributions in the small $x$ region. They can be related to gluon TMDs extracted, for example, in heavy quark pair and dijet production in deep-inelastic electron-proton scattering processes [13]. In particular, we expect to probe the same $f_{1}^{g}$ and $h_{1}^{\perp\,g}$ distributions in all such processes, because they are $T$-even TMDs. This is a very important property that still needs to be confirmed by experiments. We refer to [21] for further details.

In this paper, we will explore the relevance of final states consisting of a heavy quarkonium, like $\Upsilon$ or $J/\psi$, produced with a dilepton, be it from a $Z$ boson or from a virtual photon, in the kinematical configuration already mentioned above, such that their transverse momenta are almost back to back. In such a case, we will show how they can help extracting information about the linear polarisation of gluons. Our proposal is motivated by the fact that the detection of a dilepton may experimentally be cleaner or easier compared to the detection of a photon. For instance, the cross section for $J/\psi$ production in association with a $Z$ boson has already been studied by the ATLAS collaboration [31] and compared to theoretical predictions [32, 33, 34]. This showed that, in the ATLAS acceptance, a significant contribution from double parton scattering (DPS) is expected, which does not fall in the scope of this work. In what follows, we will assume that constraining both observed particles to be back to back makes the DPS yield small enough and we will not venture into the region of large rapidity separations where it can be dominant. Moreover, $J/\psi+\gamma$ and $\Upsilon+\gamma$ have only been investigated [35] – also by ATLAS – in the context of $H^{0}$ studies, thus at very large invariant mass where the yields are extremely small and the contributions we are after are a background to the $H^{0}$ signal. Studies in the kinematical region considered in [24] or in [36, 37] have not yet been done. $J/\psi(\Upsilon)+W^{\pm}$ would also be an option since it has also been experimentally studied – still by ATLAS [38] – but it is not dominated by gluon induced reactions [39, 40].

The paper is organised as follows: In Section 2 we discuss the gluon fusion process in the TMD factorisation approach for an arbitrary final state and analyse the general structure of the differential cross section. We then analytically calculate the partonic cross sections for each of the azimuthal structures to leading order (LO) accuracy for the final state ${\cal Q}+\ell\bar{\ell}$ where ${\cal Q}$ is a spin-triplet vector quarkonium. In Section 3, we give our numerical results for the quarkonium + $Z$ final state, while numerical results for a quarkonium + $(\ell\bar{\ell})$ final state are presented in Section 4. We draw our conclusions in Section 5.

2 Analytic Calculation within the TMD approach

In this section we consider the process $p(P_{a})+p(P_{b})\to{\cal Q}(P_{\cal Q})+\ell(l)+\bar{\ell}(\bar{l})+X$, where ${\cal Q}$ denotes a heavy quarkonium bound state, in a kinematical regime where the final state momentum $q^{\mu}\equiv P_{\cal Q}^{\mu}+l^{\mu}+\bar{l}^{\mu}\equiv P_{\cal Q}^{\mu}+P_{B}^{\mu}$ has a small transverse component, $q_{T}$, with respect to the beam axis in the proton center-of-mass (c.m.) frame. To be precise, the transverse momentum has to be much smaller than the hard scale of the process, e.g. the invariant mass $Q$ ($q^{2}=Q^{2}$) of the final-state, namely $q_{T}\ll Q$. This is the regime where TMD factorisation can be applied.

Furthermore, we assume that the underlying production mechanism of the heavy quarkonium + lepton pair is due to gluon interactions only. As such, this final state can be considered as a probe for gluon TMDs in proton collisions. In Ref. [24], we have shown that, at the LHC, this is indeed the case for the associated production of a heavy quarkonium and a real photon, instead of a lepton pair.

A second assumption we make here is that the heavy quarkonium is produced directly as a colour-singlet state. For the production of $\Upsilon-\gamma$ at the LHC, this assumption is valid, but not necessarily for $J/\psi-\gamma$ [24]. Where needed, such an assumption can be ensured by isolating the quarkonium (as done in [35]).

2.1 General Structure of the Cross Section

In a first step we present the fully differential cross section in the TMD approach (cf. Refs. [14, 29, 30, 11, 24, 15, 9, 6]) in a general form as

[TABLE]

where $a$ and $b$ are colour indices, $N_{c}=3$ is the number of colours, and the sum $\sum_{I}$ denotes the summation over those indices that can identify the particles in the final state, like their helicity. The differential phase space factor reads $\mathrm{dPS}_{n}=\prod_{i=1}^{n}\frac{d^{3}p_{i}}{(2\pi)^{3}2E_{i}}$, with $E_{i}=\sqrt{{\bm{p}}_{i}^{2}+m_{i}^{2}}$. Moreover, $\mathcal{A}$ denotes the hard scattering amplitude for an arbitrary gluon-induced process $gg\to p_{1}+...+p_{n}$, with $n$ colour-singlet particles of momenta $P_{1},...,P_{n}$ in the final state. It is factorised from the soft part contained in the integral in the second line and can be perturbatively calculated. The gluon momenta entering this amplitude are approximated as $\bar{k}_{a/b}^{\mu}=x_{a/b}P_{a/b}^{\mu}$. The longitudinal momentum fractions are set to $x_{a}=q\cdot P_{b}/P_{a}\cdot P_{b}$ and $x_{b}=q\cdot P_{a}/P_{a}\cdot P_{b}$, with $q=P_{1}+P_{2}+...+P_{n}$.

The formula in (2) describes an arbitrary gluon-induced process with a colour-singlet final state in the TMD formalism. Here, the $\tilde{\Gamma}$’s denote the non-perturbative gluonic TMD matrix elements that were properly defined in $y_{T}$ coordinate space in Ref. [9], including the renormalisation scale $\mu$ as well as an additional scale $\zeta$. The evolution in $\zeta$ of the TMD correlator $\tilde{\Gamma}$ is governed by the Collins-Soper evolution equation (cf. Refs. [3, 9]). When going from equations (2) to (2), a simple Fourier transform w.r.t. $y_{T}$ is performed. Hence, $\tilde{\Gamma}$ and $\Gamma$ are related via a Fourier transform. In particular, $\Gamma$ depends on the longitudinal gluon momentum fraction $x$ and the gluon transverse momentum ${\bm{k}}_{T}$.

The TMD correlator $\Gamma$ can be parametrised in terms of gluon TMDs (cf. Refs. [1, 42]). For an unpolarised nucleon one finds two structures of the following form,

[TABLE]

The TMD distribution $f_{1}^{g}$ can be interpreted as the distribution of unpolarised gluons in an unpolarised nucleon, while the function $h_{1}^{\perp g}$ is the distribution of linearly polarised gluons [1]. In addition, we have introduced the transverse projector $g_{T}^{\mu\nu}=g^{\mu\nu}-P^{\mu}n^{\nu}-P^{\nu}n^{\mu}$, where $P$ is the nucleon momentum and $n$ is an adjoint light-cone vector such that $n^{2}=0$ and $P\cdot n=1$. Moreover, in (3), $M$ denotes the nucleon mass. We note that the parametrisation (3) is slightly different to the one in Ref. [9] where the $M^{2}$ in (3) is replaced by $\frac{1}{2}{\bm{k}_{T}^{2}}$. Also, we have an additional factor $1/x$ in (3).

Getting back to the expression for the cross section in (2), it is useful to consider it in the c.m. frame of both incoming protons, with $P_{a/b}$ along the positive (negative) $z$-direction, e.g., $P^{\mu}_{a/b}=\sqrt{S/2}n^{\mu}_{\pm}$, with $n^{\mu}=(1,0,0,\pm 1)/\sqrt{2}$. In this frame, we can work with simple polarisation vectors of the gluons (which are approximated to be collinear to the proton momenta in the amplitudes $\mathcal{A}$ in (2)). One can easily find that these polarisation vectors acquire the following form in terms of the gluon helicities $\lambda_{a/b}=\pm 1$,

[TABLE]

These polarisation vectors lead to the common polarisation sum in the Feynman gauge,

[TABLE]

In addition, the polarisation vectors are transverse, i.e., $\bar{k}_{a}\cdot\varepsilon_{\lambda_{a}}(\bar{k}_{a})=\bar{k}_{b}\cdot\varepsilon_{\lambda_{a}}(\bar{k}_{a})=\bar{k}_{a}\cdot\varepsilon_{\lambda_{b}}(\bar{k}_{b})=\bar{k}_{b}\cdot\varepsilon_{\lambda_{b}}(\bar{k}_{b})=0$. Of course, the polarisation vectors (4) are only determined up to a phase, therefore other realisations are possible as well.

Since the contraction of Minkowski indices $\mu,\nu,\rho,\sigma$ in (2) is to be understood as a contraction in the transverse space, we can insert the polarisation sums (5) and rewrite (2) in terms of the helicity amplitudes,

[TABLE]

Here, we introduced the notation $\mathcal{A}_{\lambda_{a},\lambda_{b}}\equiv\varepsilon^{\mu}_{\lambda_{a}}(\bar{k}_{a})\ \varepsilon^{\nu}_{\lambda_{b}}(\bar{k}_{b})\ \mathcal{A}_{\mu\nu}$ for the helicity amplitudes, and the helicity gluon correlator $\Gamma_{\bar{\lambda},\lambda}\equiv(\varepsilon^{\mu}_{\bar{\lambda}})(\bar{k})\,\varepsilon^{\nu\,\ast}_{\lambda}(\bar{k})\,\Gamma_{\mu\nu}$. For an unpolarised nucleon the parameterisation in terms of the gluon helicities then takes the following form,

[TABLE]

with ${\bm{k}}_{aT}=(k_{ax},k_{ay})$ and ${\bm{k}}_{bT}=(k_{bx},k_{by})$. It is evident from (7,8) that the gluon TMD $f_{1}^{g}$ conserves the gluon helicity at the non-perturbative level, while the the distribution of linearly polarised gluons, $h_{1}^{\perp g}$, flips it. We can insert the parameterisations (7,8) into (6) and write the differential cross section in the general form

[TABLE]

where the different transverse momentum convolutions of gluon TMDs are abbreviated as

[TABLE]

Since the distribution of linearly polarised gluons carries gluon transverse-momentum-dependent prefactors in the parameterisations (7,8), these prefactors emerge again in the convolutions in (9) in the weighting factors $w$. To be specific we have

[TABLE]

We find five different structures in the cross section in (9): the first one is given by a convolution of two unpolarised gluon TMDs, while the second and fifth are given by the convolutions of two distributions of linearly polarised gluons. In the latter case the helicities of each gluon are flipped. We separated both these structures because the second typically provides azimuthally isotropic contributions, while the fifth leads to azimuthal modulations (around the beam axis) of the differential cross section. The third and fourth structures are mixed convolutions of an unpolarised and a linearly polarised gluon TMD. Here, only one of the gluon helicities is flipped.

The factors $\hat{F}_{i}$ in (9) can be calculated perturbatively since they are defined in terms of the partonic amplitudes $\mathcal{A}$ in the following way,

[TABLE]

The results in (9) - (12) are the main ones of this subsection, and are valid for an arbitrary process induced by the interaction of two gluons in the initial state, where a colour-singlet final state is produced. In the following, we will analyse a specific hadronic reaction in which a heavy quarkonium state and a lepton pair are produced back to back in transverse space.

2.2 The subprocess $gg\to{\cal Q}\,\ell\,\bar{\ell}$

2.2.1 The helicity structure

We will now compute the amplitude $\mathcal{A}$ for the process $gg\to{\cal Q}\,\ell\,\bar{\ell}$ and use these results to determine the explicit expressions of the $\hat{F}_{i}$ in (9). The first observation we make is that the lepton pair is produced either through the decay of a virtual photon or a $Z$ boson. Hence, we can decompose the amplitude further and write

[TABLE]

With the decomposition (13) at hand we can calculate the amplitude squared,

[TABLE]

We note that any explicit dependence on the single lepton momenta is hidden in the last line, where we have introduced the Lorentz-invariant leptonic tensor

[TABLE]

We now integrate over the momentum $\bar{l}$ of the antilepton, keeping in mind that $\mathrm{d}^{4}l\,\mathrm{d}^{4}\bar{l}\,\delta(l^{2})\,\delta(\bar{l}^{2})=\mathrm{d}^{4}P_{B}\,\mathrm{d}^{4}\bar{l}\,\delta(\bar{l}^{2})\,\delta((P_{B}-\bar{l})^{2})$. Since the last line of (14) is Lorentz-invariant, this is conveniently done in the center-of-mass frame of the lepton pair where the pair momentum takes the simple form $P_{B}^{\mu}=(M_{B},0,0,0)$. The delta functions are easily integrated out, and we find $\mathrm{d}^{4}l\,\mathrm{d}^{4}\bar{l}\,\delta(l^{2})\,\delta(\bar{l}^{2})\to\frac{1}{8}\,\mathrm{d}^{4}P_{B}\,\mathrm{d}\Omega$, with $\Omega$ the solid angle which determines the direction of the antilepton. Using an explicit form of the antilepton momentum $\bar{l}^{\mu}=\frac{M_{B}}{2}(1,\overset{\to}{e})$ with $\overset{\to}{e}=(\sin\theta\,\cos\phi,\sin\theta\,\sin\phi,\cos\theta)$ in this specific dilepton c.m.-frame we can perform the solid angle integration and obtain

[TABLE]

We insert this result into the quantity $\mathcal{F}$ in (14) and perform the contraction with the numerators of the weak boson-propagators. This leads to a factor $-M_{B}^{2}(g_{\mu\rho}-P_{B\mu}P_{B\rho}/M_{B}^{2})$ coming from $M_{B}^{2}\sum_{\lambda_{B}=0,\pm 1}\varepsilon_{\mu,\lambda_{B}}^{\ast}(P_{B})\,\varepsilon_{\rho,\lambda_{B}}(P_{B})$ in the numerator, which we have indentified as the sum over the helicities of the polarisation vectors of the virtual electroweak boson. Note that, since the virtual electroweak boson is massive (formally carrying a mass $P_{B}^{2}=M_{B}^{2}$), three polarisations are necessary. Utilising the polarisation sum in this way allows us to work with helicity amplitudes when calculating the process $gg\to{\cal Q}\,\gamma^{\ast}/Z$.

Finally, we also need to modify the remaining phase space $\mathrm{d}^{4}P_{\cal Q}\,\mathrm{d}^{4}P_{B}\,\delta(P_{\cal Q}^{2}-M_{\cal Q}^{2})\theta(P_{\cal Q}^{0})$ where $M_{\cal Q}$ is the mass of the heavy quarkonium state. We do so by considering the relative momenta $q^{\mu}=P_{\cal Q}^{\mu}+P_{B}^{\mu}$ and $\Delta q^{\mu}=P_{\cal Q}^{\mu}-P_{B}^{\mu}$ and note that $\mathrm{d}^{4}P_{B}\,\mathrm{d}^{4}P_{\cal Q}=\frac{1}{16}\mathrm{d}^{4}q\,\mathrm{d}^{4}\Delta q$. It is most convenient to analyse these factors in a c.m.-frame of the heavy quarkonium and the virtual electroweak gauge boson, such as the Collins-Soper (CS) frame. In this frame the gluon momenta are along the $z$-axis whereas the heavy quarkonium and electroweak gauge boson momenta have an explicit representation in this frame, $P_{\cal Q}^{\mu}=\left(\sqrt{\Lambda^{2}+M_{\cal Q}^{2}},\Lambda\,\overset{\to}{e}\right)$ and $P_{B}^{\mu}=\left(\sqrt{\Lambda^{2}+M_{B}^{2}},-\Lambda\,\overset{\to}{e}\right)$, with $\overset{\to}{e}=(\sin\theta\,\cos\phi,\sin\theta\,\sin\phi,\cos\theta)$ as before, where $\theta$ and $\phi$ are the Collins-Soper angles. From considering $q^{2}=(P_{\cal Q}+P_{B})^{2}\equiv Q^{2}>0$ we find that $\Lambda=\sqrt{\lambda(Q^{2},M_{\cal Q}^{2},M_{B}^{2})}/(2Q)$, with $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2xz-2yz$. Hence, the relative momenta take the explicit form $q^{\mu}=(Q,0,0,0)$ and $\Delta q^{\mu}=((M_{\cal Q}^{2}-M_{B}^{2})/Q,2\Lambda\overset{\to}{e})$ in the CS frame, and we conclude that $\mathrm{d}^{4}P_{\cal Q}\,\mathrm{d}^{4}P_{B}\,\delta(P_{\cal Q}^{2}-M_{\cal Q}^{2})\theta(P_{\cal Q}^{0})\to\Lambda/(4Q)\,\mathrm{d}^{4}q\,\mathrm{d}M_{B}^{2}\,\mathrm{d}\Omega$. Collecting all the results above, we rewrite (14) as

[TABLE]

2.2.2 The amplitude $\mathcal{A}^{gg\to{\cal Q}\,(Z/\gamma^{\ast})}$

As a final step, we need to calculate the amplitude $\mathcal{A}^{gg\to{\cal Q}\,(Z/\gamma^{\ast})}$. We will do so to leading order accuracy in perturbative QCD and assume a colour-singlet heavy quarkonium state. The leading order diagrams are shown in Fig. 1. In the following we restrict ourselves to heavy quarkonium states without orbital angular momentum quantum numbers, in particular to a $J/\psi$ or $\Upsilon$ state. In the colour-singlet model one typically neglects relative momenta of the heavy quark-antiquark pair in the hard part such that the wave function of the heavy quarkonium state shrinks to the origin, and the vertex of the transition of both heavy quarks forming a $J/\psi$ or $\Upsilon$ $\--$ the green blobs in Fig. 1 $\--$ reduces to a vertex $1/\sqrt{4\pi N_{c}M_{\cal Q}}\,R_{0}(0)\,(P\!\!\!/_{\cal Q}-M_{\cal Q})\,\varepsilon\!\!\!/_{\lambda_{\cal Q}}(P_{\cal Q})$, where $R_{0}(0)$ is the radial wave function of the heavy quarkonium state at the origin in the quarkonium rest frame, and $\varepsilon_{\lambda_{\cal Q}}(P_{\cal Q})$ the polarisation vector of the $J/\psi$ or $\Upsilon$ (cf. Ref. [15]). Note that we have again three polarisations $\lambda_{\cal Q}=0,\pm 1$ for a massive spin-1 particle. We also note that the coupling of the electroweak gauge boson depends on the flavor of the quark, i.e., we have a quark-gauge boson vertex of the form $\gamma^{\mu}(a_{q}+b_{q}\gamma_{5})$, with $a_{c}=-5/3+(8/3)m_{W}^{2}/m_{Z}^{2}$, $b_{c}=-1$ for a charm quark coupling to a $Z$-boson, $a_{b}=7/3-(4/3)m_{W}^{2}/m_{Z}^{2}$, $b_{b}=+1$ for a bottom quark coupling to a $Z$-boson, and $a_{c}=a_{b}=1$, $b_{c}=b_{b}=0$ for a quark coupling to a photon.

The leading order diagrams can then be calculated in a straightforward way, and we will not elaborate on the details of the calculation. Once we have obtained the analytic expressions for the helicity amplitudes $\mathcal{A}^{gg\to J/\psi[\Upsilon]\,(Z/\gamma^{\ast})}$ we insert them into (17). It is then only a minor step to extract the $\hat{F}_{i}$-prefactors in (12).

2.2.3 The azimuthal dependence of the differential cross section

To go further, we decompose, in the CS frame, the differential cross section in terms of factors $\hat{F}_{i}$ :

[TABLE]

where the $\hat{F}_{i}$ are functions of the $\theta$ CS-angle, the parameters $\alpha\equiv M_{\cal Q}/Q$ and $\beta=M_{B}/Q$. The prefactors $\hat{F}_{i}$ have the form

[TABLE]

where $\tilde{C}_{\gamma^{\ast}}=-4\pi\alpha_{\mathrm{em}}e_{q}$ and $\tilde{C}_{Z}=m_{Z}^{2}G_{F}a_{q}/(2\sqrt{2})$, and $\Lambda=\frac{Q}{2}\sqrt{\lambda(1,\alpha^{2},\beta^{2})}$. The functions $\hat{f}_{i}$, $\hat{D}$ acquire the following form in terms of the auxiliary variables $A\equiv 1+\alpha^{2}-\beta^{2}$ and $B\equiv\sqrt{A^{2}-4\alpha^{2}}=\sqrt{\lambda(1,\alpha^{2},\beta^{2})}\geq 0$,

[TABLE]

as well as

[TABLE]

with

[TABLE]

Note that in the limit of real photon production, that is, $M_{B}\to 0$ or $\beta\to 0$, we recover the results of the $\hat{F_{i}}$ that were found in Ref. [24]. We find that the partonic prefactor $\hat{f}_{2}$ vanishes in the real photon limit $\--$ a unique feature of this particular final state.

Eqs. (18-22) constitute the main analytical results of this work. It is however instructive to analyse the cross section that is integrated over the CS angles including a possible azimuthal weighting factor. These weighting factors may enable us to disentangle the various azimuthal contributions in (18). For example, analoguously to Ref. [24], we can deduce the following weighted cross sections from (18),

[TABLE]

Note that we use the same symbols $\hat{F}_{i}$ for the integrated or weighted partonic prefactors, i.e., $\hat{F}_{1,2}(Q,\alpha,\beta)=2\pi\int\mathrm{d}\theta\,\hat{F}_{1,2}(Q,\alpha,\beta,\theta)$ and $\hat{F}_{3,4}(Q,\alpha,\beta)=\pi\int\mathrm{d}\theta\,\hat{F}_{3a/b,4}(Q,\alpha,\beta,\theta)$. The $\theta$-integration can be performed analytically, and we find the following integrated partonic prefactors utilising Eqs. (19-22),

[TABLE]

[TABLE]

where we have used the definition

[TABLE]

The advantage of the quantities $N^{(i)}$ is that they are $\--$ in principle $\--$ Lorentz-invariant after having integrated out the CS-angles. Hence, they can be evaluated directly in the hadron c.m.-frame, i.e., the lab frame at the LHC.

Below we use these formulae to estimate the size of the effects from linearly polarised gluons that one can expect at the LHC by measuring the associated heavy-quarkonium + $Z$ - final state.

3 Numerical Predictions for an associated quarkonium + $Z$ - final state

We can numerically compare the relative size of both contributions from unpolarised and linearly polarised gluons to the angular-integrated cross section $N^{(0)}$ in (23) by considering the LO ratios $\hat{F}_{2,3,4}(Q,\alpha,\beta)/\hat{F}_{1}(Q,\alpha,\beta)$ from Eqs. (24 - 27). In this section, we focus on the production of a (quasi)-real $Z$-boson, i.e., dileptons with an invariant mass around the $Z$-pole mass $M_{B}\simeq m_{Z}$. To this end, we consider a $M_{B}$-bin around $m_{Z}$ with a bin size of, say, 2 GeV. Hence, we define cross sections for real $Z$-boson production in the following way,

[TABLE]

The structures of the corresponding quantities $N^{(i)}$ in (23) remain valid for real $Z$-boson production, but with integrated partonic prefactors

[TABLE]

The LO ratios $\hat{F}_{2,3,4}(Q,\alpha,\beta)/\hat{F}_{1}(Q,\alpha,\beta)$ for real $Z$-boson productions are shown in Fig. 2. In the left plot we present our result for an associated $\Upsilon$ state with mass $m_{\Upsilon}=9.46\,\mathrm{GeV}$. In the right plot, we have shown it for a $J/\psi$ state with mass $m_{J/\psi}=3.1\,\mathrm{GeV}$. First of all, we observe that the ratios are rather small in general: the ratio $F_{2}/F_{1}$ is about half of a percent at most for $\Upsilon$-production and even smaller ($\leq 10^{-3}$) for $J/\psi$ production. One can expect that the convolution $\mathcal{C}[w_{2}\,h_{1}^{\perp g}\,h_{1}^{\perp g}]$ in (23) from linearly polarised gluons does not exceed in size the convolution $\mathcal{C}[f_{1}^{g}\,f_{1}^{g}]$ [29, 8]. For example, the models of Ref. [8] indicate that the (scale-dependent) ratio $R=\mathcal{C}[w_{2}\,h_{1}^{\perp g}\,h_{1}^{\perp g}]/\mathcal{C}[f_{1}^{g}\,f_{1}^{g}]$ is at most about $2/3$ for a small scale $Q\sim 3\,\mathrm{GeV}$, but typically (much) smaller for larger scales. Hence, it is not unreasonable to neglect the contribution from linearly polarised gluons for the quantity $N^{(0)}$ in (23) and approximate to good accuracy,

[TABLE]

In a next step we evaluate the quantities $N_{Z}^{(i)}$ for real $Z$-boson production in the hadron c.m.-frame where $P^{\mu}_{a/b}=(\sqrt{S}/2)(1,0,0,\pm 1)$ and $q^{\mu}=(\sqrt{Q^{2}+{\bm{q}}_{T}^{2}}\,\cosh Y,{\bm{q}}_{T},\sqrt{Q^{2}+{\bm{q}}_{T}^{2}}\,\sinh Y)$. Here, $Y$ denotes the rapidity of quarkonium-dilepton pair, i.e., the rapidity of the final state. Also, $\mathrm{d}^{4}q=Q\,\mathrm{d}Q\,\mathrm{d}Y\,\mathrm{d}^{2}{\bm{q}}_{T}$. We then investigate the following distributions (ratios) that where already proposed in Ref. [24],

[TABLE]

Note that, for these ratios, the transverse momentum of the final state ${\bm{q}}_{T}$ $\--$ the quarkonium-dilepton pair’s transverse momentum imbalance $\--$ is restricted to be smaller than $Q/2$ in order to roughtly fulfill the TMD factorisation requirement $|{\bm{q}}_{T}|\ll Q$.

In order to estimate the size of the effects of the linearly polarised gluons hidden in the observables $S_{Z}^{(2)}$ and $S_{Z}^{(4)}$, it is instructive to first calculate the ratios of the perturbative prefactors $\hat{F}_{3}^{Z}/\hat{F}_{1}^{Z}$ and $\hat{F}_{4}^{Z}/\hat{F}_{1}^{Z}$. We show these ratios in Fig. 2 as well, plotted vs. the invariant mass $Q$ for a $\Upsilon$ and $J/\psi$ state. We find rather small ratios, around $10^{-4}$ for $\hat{F}_{3}^{Z}/\hat{F}_{1}^{Z}$ and $5\times 10^{-5}$ for $\hat{F}_{4}^{Z}/\hat{F}_{1}^{Z}$ at the peak around $Q=120\,\mathrm{GeV}$ for a $\Upsilon$-particle. The ratios for $J/\psi$-production are even smaller, of the order of $10^{-6}$. This is in contrast to the $\Upsilon+\gamma$ final state discussed in Ref. [24] where the corresponding ratios are about $0.05$ and $0.03$, respectively, for an invariant mass $Q=20\,\mathrm{GeV}$. This already suggests that the final state containing a heavy quarkonium state and a real $Z$-boson may not be sufficiently suitable to study the distribution of linearly polarised gluons, $h_{1}^{\perp g}$. A more quantitative statement about the feasibility of measurements of the azimuthal observables $S^{(2)}$ and $S^{(4)}$ can be given through an estimate of the convolution integrals in Eqs. (32,33,34). We use the same model Ansätze for the unpolarised gluon $f_{1}^{g}$ that were also adopted in the predictions for a $\Upsilon+\gamma$ state in Ref. [24], where two parameterisations of the unintegrated gluon distribution (UGD) have been used as an input for the TMD gluon function $f_{1}^{g}$. Although the UGD established for small-$x$ physics may not be identified in general with the unpolarised gluon TMD, such an input serves as a first numerical estimate of the size that can be expected for the distributions $S^{(i)}$. In particular, we use the Set B0 solution to the CCFM equation with an initial condition based on the HERA data from Refs. [43, 44] and the KMR parameterisation from Ref. [45] for the UGD. In Ref. [24] the saturation of the positivity bound [1] for the distribution of linearly polarised gluons was assumed, i.e., $h_{1}^{\perp g}(x,{\bm{k}}_{T}^{2})=2M^{2}/{\bm{k}}_{T}^{2}\,f_{1}^{g}(x,{\bm{k}}_{T}^{2})$. We rely on the same assumption in this work as well. In Fig. 3 the azimuthal ${\bm{q}}_{T}$-distribution $S^{(2)}$ and $S^{(4)}$ from Eqs. (33,34) are shown for real $Z$-bosons and are of the size of about $10^{-8}$, which is about four orders of magnitudes smaller than the corresponding contributions for a $\Upsilon+\gamma$ final state. The ${\bm{q}}_{T}$-integrated azimuthal observables $\int d^{2}{\bm{q}}_{T}\,S^{(2)}$ and $\int d^{2}{\bm{q}}_{T}\,S^{(4)}$ amount to roughly $0.007\%$ and $0.001\%$, respectively $\--$ three orders of magnitude smaller than for $\Upsilon+\gamma$. Hence, we conclude that it will be very difficult to access the linearly polarised gluons with a $\Upsilon+Z$ final state. Having said that, this particular final state might be a suitable candidate for studying experimentally the unpolarised gluon TMD on its own through the measurement of the azimuthally independent distribution $S^{(0)}$. The difference between this observable for a $\Upsilon+Z$ final state and a $\Upsilon+\gamma$ final state is that the invariant mass varies. In Ref. [24] the distribution $S^{(0)}$ was studied at $Q=20\,\mathrm{GeV}$ for a $\Upsilon+\gamma$ final state, while in Fig. 3 this distribution is shown for $Q=120\,\mathrm{GeV}$. The larger invariant mass is a consequence of the mass of the $Z$-boson. Hence, the ${\bm{q}}_{T}$-distribution $S^{(0)}$ is much broader for $\Upsilon+Z$. In addition, the detection of a dilepton pair at the $Z$-pole is experimentally much favorable in contrast to real photon detection because an isolation procedure is needed in the latter case.

4 Numerical prediction for quarkonium + a dilepton from an off-shell photon

In this section we repeat the steps of the previous section but we consider lepton pairs with a relatively small invariant mass $M_{B}\in[5\,\mathrm{GeV},\,7\,\mathrm{GeV}]$ between the $J/\psi$ and $\Upsilon$ families. This mass range is far away from the $Z$-pole mass $m_{Z}$, and dilepton creation from decays of virtual photons instead of $Z$-bosons should dominate.

The advantage is that we can investigate lower values of $Q$ with a minimum final state invariant mass $Q_{\mathrm{min}}=7\,\mathrm{GeV}+M_{\cal Q}$. In fact, we expect to observe some similarities with the associated $\mathcal{Q}+\gamma$ final state discussed in Ref. [24]. The theoretical formulae for this kinematic range are the same as for real $Z$-boson production of the last section, except that the integration region in (30) is $5\,\mathrm{GeV}\leq M_{B}\leq 7\,\mathrm{GeV}$.

In Fig. 4, we plot the ratios $F_{i}/F_{1}$ for associated Quarkonium - dilepton production with small dilepton masses vs. the final state invariant mass $Q$. We observe that these ratios are considerably larger (by 1-2 orders of magnitude) than for real $Z$-boson production, but still smaller than for real photon production. Overall, one can say that the associated Quarkonium - dilepton final state is rather sensitive to the dilepton mass $M_{B}$. In particular, the ratio $F_{2}/F_{1}$ indicates that the prefactor $\hat{F_{2}}$ characterising the linearly polarised gluons may be a few percent of the factor $\hat{F}_{1}$ for lower $Q$. Nevertheless, we still consider the approximation (31) to be justified.

We then calculate the $\bm{q}_{T}$-distributions $S^{(i)}$ for small dilepton masses and show our results in Fig. 5. Since the distribution $S^{(0)}$ in the upper panel does not depend on the specific final state, we obtain the same result for a final state invariant mass $Q=20\,\mathrm{GeV}$ as for real photon production in Ref. [24]. The ratios $F_{3,4}/F_{1}$ are smaller compared to real photon production as mentioned before, and this results in smaller azimuthal $\bm{q}_{T}$-distributions $S^{(2)}$ and $S^{(4)}$. We estimate the overall $\bm{q}_{T}$-integrated $\cos 2\phi$ effect from linearly polarised gluons to be about $0.5\%-0.6\%$ in the dilepton mass range $M_{B}\in[5\,\mathrm{GeV},\,7\,\mathrm{GeV}]$ for an invariant mass $Q=20\,\mathrm{GeV}$, while the $\cos 4\phi$ modulation amounts to $0.04\%-0.045\%$. Although these rates are smaller by an order of magnitude compared to real photon production, the experimental advantage of cleaner final state may outweigh the disadvantage of a smaller effect. This discussion is however beyond the scope of the present analysis.

Finally, we also investigate an intermediate dilepton mass range $M_{B}\in[20\,\mathrm{GeV},25\,\mathrm{GeV}]$. The corresponding ratios of the LO prefactors $F_{i}/F_{1}$ are shown in Fig. 6 where we observe a reduction of a factor of about 3-4 compared to the small dilepton mass range $M_{B}\in[5\,\mathrm{GeV},\,7\,\mathrm{GeV}]$, both for a $\Upsilon$ and $J/\psi$. Consequently, also the azimuthal $\bm{q}_{T}$-distributions $S^{(2)}$ and $S^{(4)}$, taken at a larger invariant final state mass $Q=40\,\mathrm{GeV}$ and shown in Fig. 7, are smaller. The overall effect amounts to $0.13\%-0.15\%$ for the $\cos 2\phi$ modulation and $0.01\%-0.012\%$ for the $\cos 4\phi$ modulation.

5 Conclusions

In this paper, we have presented the treatment of an arbitrary process induced by gluon-gluon fusion with a colour-singlet final state in the TMD approach. Using the helicity formalism, we have analysed the general structure of the fully differential, unpolarised cross section for this process in terms of TMD distributions of unpolarised and linearly polarised gluons inside an unpolarised nucleon. We have then calculated the partonic cross sections underlying the azimuthally independent and dependent structures of the hadronic cross section for a specific final state: a heavy quarkonium in a colour-singlet state and a real $Z$-boson. We have found that, in contrast to quarkonium production associated with a photon, the azimuthally dependent contributions are strongly suppressed with respect to the azimuthally independent ones. Therefore, the distribution of linearly polarised gluons in the nucleon is most probably not experimentally accessible for a quarkonium final state that is associated with a real $Z$-boson. The unpolarised gluon TMDs may however be studied by means of the azimuthally independent transverse momentum distribution of the process $pp\to\Upsilon ZX$, which can be directly measured at the LHC.

We have also investigated the associated quarkonium + dilepton production for small and medium dilepton masses. These dileptons are predominantely generated from virtual photon decays. Although the effects from linearly polarised gluons are still smaller than for the real-photon case, a TMD extraction might be done for dilepton masses between the $J/\psi$ and $\Upsilon$ masses at the LHC in view of the recent experimental studies of $\Upsilon+\Upsilon$ [46] and $J/\psi+\Upsilon$ [47] production.

Acknowledgements

We thank D. Boer, W. den Dunnen, M.G. Echevarria, T. Kasemets, H.S. Shao, A. Signori, J.X. Wang for helpful discussions. The work of CP is supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 647981, 3DSPIN). The work of JPL is supported in part by the CNRS-IN2P3 (project TMD@NLO). The work of MS is supported in part by the Bundesministerium für Bildung und Forschung (BMBF) grant 05P15VTCA1.