Associated production of a quarkonium and a Z boson at one loop in a quark-hadron-duality approach

In view of the large discrepancy about the associated production of a prompt $J/\psi$ and a $Z$ boson between the ATLAS data at $\sqrt{s}=8$ TeV and theoretical predictions for Single Parton Scattering (SPS) contributions, we perform an evaluation of the corresponding cross section at one loop accuracy (Next-to-Leading Order, NLO) in a quark-hadron-duality approach, also known as the Colour-Evaporation Model (CEM). This work is motivated by (i) the extremely disparate predictions based on the existing NRQCD fits conjugated with the absence of a full NLO NRQCD computation and (ii) the fact that we believe that such an evaluation provides a likely upper limit of the SPS cross section. In addition to these theory improvements, we argue that the ATLAS estimation of the Double Parton Scattering (DPS) yield may be underestimated by a factor as large as 3 which then reduces the size of the SPS yield extracted from the ATLAS data. Our NLO SPS evaluation also allows us to set an upper limit on $\sigma_{\rm eff}$ driving the size of the DPS yield. Overall, the discrepancy between theory and experiment may be smaller than expected, which calls for further analyses by ATLAS and CMS, for which we provide predictions, and for full NLO computations in other models. As an interesting side product of our analysis, we have performed the first NLO computation of $d\sigma / dP_T$ for prompt single-$J/\psi$ production in the CEM from which we have fit the CEM non-pertubative parameter at NLO using the most recent ATLAS data.

E-mail: Jean-Philippe.Lansberg@in2p3.fr, huasheng.shao@cern.ch Abstract: In view of the large discrepancy about the associated production of a prompt J/ψ and a Z boson between the ATLAS data at √ s = 8 TeV and theoretical predictions for Single Parton Scattering (SPS) contributions, we perform an evaluation of the corresponding cross section at one loop accuracy (Next-to-Leading Order, NLO) in a quarkhadron-duality approach, also known as the Colour-Evaporation Model (CEM). This work is motivated by (i) the extremely disparate predictions based on the existing NRQCD fits conjugated with the absence of a full NLO NRQCD computation and (ii) the fact that we believe that such an evaluation provides a likely upper limit of the SPS cross section. In addition to these theory improvements, we argue that the ATLAS estimation of the Double Parton Scattering (DPS) yield may be underestimated by a factor as large as 3 which then reduces the size of the SPS yield extracted from the ATLAS data. Our NLO SPS evaluation also allows us to set an upper limit on σ eff driving the size of the DPS yield. Overall, the discrepancy between theory and experiment may be smaller than expected, which calls for further analyses by ATLAS and CMS, for which we provide predictions, and for full NLO computations in other models. As an interesting side product of our analysis, we have performed the first NLO computation of dσ/dP T for prompt single-J/ψ production in the CEM from which we have fit the CEM non-pertubative parameter at NLO using the most recent ATLAS data.

Introduction.
With the advent of the LHC, the observation of the associated production of a quarkonium and a vector boson became possible. Pioneering studies of CDF at the Fermilab-Tevatron [1,2] were motivated by the search for a charged H ± decaying in a pair of Υ + W ± for instance. In 2014, ATLAS succesfully observed for the first time the simultaneous production of J/ψ + W [3] and J/ψ + Z [4]. As for now, the similar processes involving bottomonia have never been observed 1 As its customary with quarkonium physics, these new measurements brought to light significant -to say the least-tensions with the existing theories. Whereas it is very likely that such simultaneous productions may come from two independent parton scatteringsalso called Double Parton Scattering (DPS) -the ATLAS collaboration concluded that the usual production from Single Parton Scattering (SPS) was also relevant and as a matter of fact was, according to their analysis, significantly above existing theoretical predictions (see [15,16] for J/ψ + Z and [17,18] for J/ψ + W ). In the case which interests us here, namely J/ψ + Z, the SPS yield extracted by ATLAS [4] is more than five times larger 1 Another class of reactions which is of interest is that of the production of a quarkonium + a photon.
It is has been proposed to constrain the quarkonium-production mechanisms [5][6][7][8][9], to measure different characteristics of the gluon content of the proton [10,11] as well as to probe the H 0 coupling to the heavyquarks [12,13]. This was the motivation of the sole experimental study of this channel by ATLAS [14] but it mainly focused on the search for a signal from a H 0 decay rather than to the QCD continuum (no cross-section extraction was performed) which interests us here.
(corresponding to 2-σ deviation) than the largest of the theory evaluations from NRQCD (with Colour-Singlet (CS) and/or Colour-Octet (CO) contributions) 2 .
Other quarkonium-associated-production channels have also been investigated. 30 years after the pionneering analyses of NA3 [20,21], J/ψ-pair production has been analysed by the LHCb [22], CMS [23] and ATLAS [24] collaborations at the LHC as well as by the D0 collaboration [25] at the Tevatron. They are all compatible with CS contributions only (known up to NLO accuracy [26][27][28]) at small rapidity separations, ∆y, whereas they point at a significant DPS contributions for increasing ∆y, compatible with a σ eff below 10 mb. We guide the reader to [29] for a detailed discussion of these different results and to [30,31] for recent LO NRQCD analyses. Υ + J/ψ production has also been observed by the D0 collaboration [32] with a claim that the yield is highly dominated by DPS contributions (see [33] for a complete and up-to-date discussion of the theoretical aspects of such a reaction). Similar conclusions have been drawn by LHCb for J/ψ+charm [34] and Υ+charm [35] production although compatible with σ eff larger than 10 mb. It is worth emphasising that the D0 and ATLAS J/ψ-pair analyses [24,25] are the only for quarkonia where DPS and SPS were separated based on kinematical variables.
Motivated by the discrepancy uncovered by ATLAS, we perform here the very first complete evaluation of the SPS yield at NLO for the production of a J/ψ associated with a Z boson in pp collisions under the assumption of quark-hadron duality which, in the case of quarkonium production, is referred to as the Colour-Evaporation Model (CEM). Such a computation has in fact never been performed at LO. This allows us to question the size of the DPS yield assumed by ATLAS based on their W + 2-jet analysis [36] and to propose a solution for the present puzzle with a DPS yield roughly three times larger (yet compatible with the 1-σ upper value set by ATLAS).
As a side product of this complete analysis, we provide another novel NLO analysis bearing on single-J/ψ production, namely that of the single-J/ψ P T spectrum at NLO via the computation at NLO (α 4 s ) of the P T spectrum of a J/ψ recoiling against a parton, which we use to constrain the sole non-pertubative CEM parameter in a consistent way using the same P T region for single-J/ψ data as the one for J/ψ + Z. Usually, such a parameter is fit from the total J/ψ yield which is known at NLO (α 3 s ). The structure of the paper is as follows. In section 2, we explain the purpose of using a simple approach as the CEM in the context of heavy-quarkonium production. In particular, we highlight its limitation as well as its interest for the present study. In section 3, we discuss its implementation in the NLO framework which we have used; we detail the computation of the single-J/ψ P T spectrum at NLO and compare the CEM non-perturbative parameter which we obtained here with previous works. In section 4, we 2 We noticed that the theory numbers quoted in [4] have probably been misconverted into the ratio R which is compared to the data (see below). We have made thus sure that our computation of R is indeed compatible with what is implied in the ATLAS definition [19]. As a check, we have recomputed the cross sections and R in NRQCD at Leading Order (LO) and included all the possible feed-downs. Despite the wrong quoted numbers by ATLAS, what we found however ends up to be similar than the quoted ones and this does not therefore change the conclusion that SPS contributions are small and that the most optimistic NRQCD upper values are not larger than one sixth of the DPS-subtracted ATLAS measurements. present our results for prompt J/ψ + Z. We start with a discussion of the P T -integrated yield, then discuss the P T dependence of both SPS and DPS contributions which finally allows us to explain why we believe that the ATLAS data in fact allows for a DPS yield three times as large as they initially assumed. Section 5 gathers our conclusions.

Why the Colour-Evaporation Model ?
In parallel to the aforementioned advances in the study of quarkonium associated production, more precise data for single-quarkonium production are flowing from the LHC, including cross-section measurements at larger P T , improved measurements of feed-down fractions, more reliable multi-dimensional polarisation measurements (see [37] for a recent review) and, last but not least, the first measurement of the production cross-section of the spin-singlet η c [38]. At the great surprise of some, such η c data happened to be very well described by postdictions [39][40][41] with the sole CS contributions (also called CSM) leaving nearly no room for CO contributions in this channel. Such a constraint, translated to the J/ψ case via Heavy-Quark-Spin Symmetry (HQSS), could only be met by assuming a rather small value of the LDMEs for O J/ψ ( 1 S [8] 0 ) , which in turn induces, in order to reproduce the large-P T J/ψ spectra at the LHC and the Tevatron, a somewhat large value for O J/ψ ( 3 S [8] 1 ) , still in the ballpark of the NRQCD Velocity-Scaling Rules (VSR) [40]. However, these constraints are completely at odds with the assumption made in some recent works [42,43] of a very large O J/ψ ( 1 S [8] 0 ) in order to obtain an unpolarised J/ψ yield at large P T . It is also going against the trend obtained with the global-fit approach of [44] including hadroproduction data at lower P T as well as γp, γγ and e + e − data (which however fails to describe the polarisation of large-P T J/ψ). As noted in [39], the η c data points at a O J/ψ ( 1 S [8] 0 ) 10 times smaller than their initial fit. These observations "led [the authors of [39] ] to conclude that either the universality of the LDMEs is in question or that another important ingredient to current NLO NRQCD analyses has so far been overlooked." For instance, NRQCD, for a reason thatwe do not have uncovered yet, may not be reliable for polarisation observables.
Along the same lines, it was shown in [45] that the NRQCD universality was severely challenged by the P T -integrated cross sections. In the same work, the CEM was shown to describe reasonnably well the world data (see also next section). In view of this and the discussion above, it appears to us as legitimate to consider the CEM 3 as one such models on the market allowing one to explore the possible production mechanisms in a rather unbiased way. The main advantage of using the CEM is that it is possible to perform a complete one-loop computation whereas NLO NRQCD computation involving all the relevant channnels are much more complex (without showing necessarily more stability).
Beside such encouraging comparisons of the CEM P T -integrated yields, most of the Tevatron, RHIC and LHC data for P T -differential cross sections have been confronted to the CEM predictions (see [37] for a representative selection). If one disregards the spectrum at low P T where a phenomenological parton-k T -smearing procedure is usually applied to reproduce the cross section, the CEM usually tends to overshoot the data for increasing P T . In order to address this issue, different mechanisms [47][48][49] have been considered (see [50] for a brief overview) but none was the object of a consensus. In fact, the origin of the discrepancy is obvious and arises from the appearance at α 3 s of leading-P T topologies scaling like P −4 T , just as those associated with the 3 S [8] 1 octet states in NRQCD. Whereas these were originally thought to solve the ψ(2S) surplus at the Tevatron, recent NRQCD fits to the LHC and Tevatron data indicate that they have to be somewhat damped down. In NRQCD, it happens through a partial cancellation between both channels which show leading-P T contributions at NLO, namely the 3 S [8] 1 and 3 P [8] J states, opening the possibility for a dominant 1 S [8] 0 contribution in agreement with a softer P T spectrum 4 . This is precisely why the η c data are troublesome since they tend to constrain the importance of the 1 S states. In the CEM, owing to the simplicity of the model, no such cancellation can happen. The fragmentation contributions are obviously there, at any non-trivial order where the computation is carried out.
Given these discrepancies, one may wonder whether using the CEM makes sense. Our answer is "yes" as long as complete NLO predictions are not available in NRQCD and, most importantly, as long as the different NLO NRQCD fits on the market give extremely large uncertainties for quarkonium-associated-production channels. Even in some cases, some fits yield to negative cross sections [9], which illustrates our current lack of understanding of the quarkonium hadronisation. Having a simple model at hand can then be a very useful tool to investigate tensions between data and experiments case by case 5 .
In view of all these arguments, we believe the CEM to be in fact currently the best model to investigate the ATLAS excess for J/ψ + Z production since a complete NLO NRQCD computation is lacking and since the CEM would probably provide an upper theory limit given the predominance of the fragmentation channels for this process.

The CEM and the P T -integrated yields
The CEM can be seen as the application [51,52] of the quark-hadron-duality principle to quarkonium production. The cross sections for quarkonium production are obtained by integrating the cross section to produce a QQ pair in an invariant-mass region where its hadronisation into a quarkonium is likely. In practice, one considers that it can occur 4 The 1 S [8] 0 channel at NLO in principle also exhibits fragmentation-like topologies which however do not result in a P −4 T scaling since the radiated gluon from the heavy-quark pair is not soft. 5 However, if it happens in the future, after confronting the CEM to a number of new observables, that the tensions between the PT -differential cross section of single-quarkonium and the CEM are the only serious ones, it could mean that the underlying assumption of the CEM may not be too far from reality. In such a case, it could be worth devoting some efforts to investigate whether there are justifications and means to solve these tensions about the PT dependence, like an alteration of the PT spectrum by the feed-downs or a kinematical bias between the computed PT of the pair and the one of the quarkonium eventually produced. Yet, as for now, we will consider the CEM as a simple and robust model which probably tends to produce a little too many quarkonia at increasing PT when fragmentation channels are open. between 2m Q and the threshold to produce open-heavy-flavour hadrons, 2m H . One subsequently multiplies this partial heavy-quark cross section by a phenomenological factor accounting for the probability that the pair eventually hadronises into a given quarkonium state, P Q . All this amounts to consider P Q can be paralleled to the LDMEs in NRQCD, but for the fact that its size can be guessed if one assumes a simplified statistical-hadronisation scenario. Owing to the simplicity of the model, the direct or prompt yields are obtained from the same computation but with a different overall factor. First, one expects [53] that one ninth -one CS QQ configuration out of 9 possibleof the open-charm cross section in this invariant-mass region eventually hadronises into a "stable" quarkonium. In the case of J/ψ, beside this factor 9, it was argued [53], with a LO computation, that a simple statistical counting giving where the sum over i runs over all the charmonium states below the DD threshold, could describe the existing data in the late 90's. NLO fits were then performed by Vogt in [54] on data up to √ s = 62 GeV with P direct J/ψ lying between 1.5 % and 2.5 %. This simple statistical counting rule works remarkably well for J/ψ, whereas it would not work for P -waves as discussed in [50,55] as well as for ψ(2S). This shows the limit of the model in accounting for differences in the hadronisation of different quarkonium states. For the Υ, the corresponding quantity is fit with a similar size, between 2 % and 5 %. Following the state-counting argument, one would however expect a smaller number than for J/ψ. At this stage, it is important to reiterate that Eq. (3.2) ignores phase-space constraints in the hadronisation process. What the CEM really predicts is that P direct Q should be process independent.
A NLO comparison with data from fixed-target as well as from colliders is shown on Fig. 1 for the J/ψ case. Curves with different choices of the charm-quark mass and of the PDFs are shown. First we see that the mass dependence is nearly completely absorbed in the fit value of P prompt J/ψ . We stress that m Q impacts the evaluation of the cross-section principally via the range of integration over the pair invariant mass. Yet, the energy dependence is rather similar for both mass choices. One also sees the significant dependence on the PDF set. This is due to the rather low scale µ R of this process and the very low x values (x M ψ / √ s) reached. Such an effect does not appear at higher x and at larger scales as for J/ψ + Z.
A thorough study of the connection between heavy-flavour and quarkonium production in the CEM can be found in [59]. Along the lines of this analysis, we will use m c = 1.27 GeV. Results with m c = 1.5 GeV are sometimes slightly different but never such as to modify the physics conclusion.   [56][57][58] as a function of √ s (see also [45]).

The single-J/ψ P T spectrum at NLO
We turn now to the discussion of the P T spectrum of single J/ψ's. As mentioned above, the CEM is expected to predict too hard a P T spectrum. Yet, in all Tevatron and LHC computations [54,60], the hard-scattering matrix element, which is employed, is at one loop for inclusive heavy-quark production, namely α 3 s . In fact, it is based on the well-known MNR computation [61] using the specific invariant-mass cut of Eq. (3.1). At this order, the sole graphs contributing to the production of the heavy-quark pair (with or without invariant-mass cut) with a finite P T are those from 2 → 3 processes. This de facto excludes any loop contribution. As such, all these existing computations at finite P T of the pair are effectively Born-order/tree-level computation from gg[qq] → (QQ)g or gq → (QQ)q. It is therefore legitimate to wonder whether the resulting P T spectrum could be affected by large QCD α 4 s corrections, effectively NLO and not NNLO for this quantity. In particular, one could wonder whether the data can be better described at NLO and whether P NLO J/ψ is different than P LO J/ψ ? In addition to be a significant advance in the CEM usage, such a computation is in fact relevant for our study since the J/ψ measured by ATLAS with a Z are in the range where the single-J/ψ data starts to deviate from the CEM prediction at α 3 s and where such NLO corrections could matter both in the determination of the CEM parameter and in the hard-part computation for J/ψ + Z itself.
Given the straightforward connection between the CEM and heavy-quark production, such a computation is in fact not too demanding with modern tools of automated NLO frameworks, at the minimum cost of some slight tunings. In particular, we have used Mad- Graph5_aMC@NLO [62] to perform our (N)LO CEM calculations for J/ψ + a recoiling parton with a finite P T (and then of course for J/ψ + Z) 6 . As already stated above, we have taken m c = 1.27 GeV, while we checked that m c = 1.5 GeV would only marginally change our results provided that the non-perturbative CEM parameter is chosen coherently. As what regards the parton distribution function (PDF), we have used the NLO NNPDF 2.3 PDF set [57] with α s (M Z ) = 0.118 provided by LHAPDF [66]. In this case, since the heavy-quark mass dependence is de facto absorbed in the CEM parameter, the main theoretical uncertainties 7 are coming from the renormalisation µ R and factorisation µ F scale variations which account for the unknown higher-order corrections. In practice, we have varied them within 1 2 µ 0 ≤ µ R , µ F ≤ 2µ 0 where the central scale µ 0 is the transverse mass of the J/ψ in J/ψ + parton and the mass of Z boson M Z in J/ψ + Z. It has been shown [16] that the scale M Z is close to where the result is the most stable for J/ψ + Z at NLO in the CSM.
Since we wish to eventually use P NLO,prompt J/ψ for our study of the ATLAS J/ψ +Z data, we have restricted our study to the kinematical condition of the latest ATLAS single-J/ψ data at √ s = 8 TeV [67] corresponding to 11.4 fb −1 and which are of course much more precise than previous data sets.
Fitting this set with m c = 1.27 GeV, we obtain P LO,prompt J/ψ = 0.014 ± 0.001 and P NLO,prompt J/ψ = 0.009 ± 0.0004. From these, we can deduce that the K factor affecting the P T slope is close to 1.6. As previously discussed, the CEM yields start to depart from the data when P T increases (see Fig. 2). This is happening in the region where fragmentation contributions are dominant as in J/ψ +Z which makes us believe that the CEM will indeed provide an upper limit on SPS J/ψ + Z computations.
We note that the behaviour is similar at LO and NLO and that the uncertainties at NLO are smaller than at LO. We recall that the LO results (i.e. for prompt J/ψ + a recoiling parton) would coincide with the P T spectrum obtain from NLO code for inclusive prompt J/ψ results [54,60]. This illustrates the added value of this first NLO CEM study of the P T spectrum of single J/ψ using a one-loop evaluation of J/ψ + a recoiling parton. Similar studies for other quarkonia will be the object of a separate work. For completion, let us add that the χ 2 dof is 1.57 at NLO and 0.51 at LO, essentially because the LO scale uncertainties are larger.
The value of P is about a factor of 2-3 smaller than that obtained from the P T -integrated total yields (see Fig. 1). The trend is opposed to that of NRQCD where the LDME fit values from the P T -differential yield systematically overshoot that fit from the P T -integrated total yields [45,68]. Thus even for the CEM, the universality P prompt J/ψ seems to be challenged. In the following, we will naturally use the value of P fit to the P T -differential yields.

The ATLAS comparison with theory
Let us now move on to the process of interest of this analysis, namely J/ψ + Z production at the LHC where the J/ψ is promptly produced, thus not from a b-hadron decay.
Such a process has been studied in the past in NRQCD under the SPS mechanism. Only two NLO analyses exist, that of Ref. [15] considered both CS and CO contributions, but strictly speaking is not a complete NLO NRQCD analysis since the 3 P [8] J transitions were disregarded. It would only be complete if the corresponding LDME was negligibly small but, as discussed above, this would contradict single-J/ψ data. Following this study, another one considering only CS contributions appeared [16]. It corrected a mistake in [15], presented first polarisation results and argued, owing to the scale dependence of the process, that a reasonable choice for the scales is rather M Z than M 2 J/ψ + P 2 T as used in [15]. When the ATLAS collaboration released its study, they attempted to compare their data to these predictions. However they did not directly compare them to their yield since they noted that a non-negligible part of the yield was probably from DPS contributions. This conclusion was motivated by two observations. First the distribution of the events as a function of the azimuthal angle between both detected particles, ∆φ, was showing a plateau close to 0, whereas a dominant SPS yield would show a peak at π, i.e. for backto-back events. Second, under the simple assumption that the DPS yield comes from two uncorrelated scatterings, they evaluated it with the rudimentary pocket formula with σ eff extracted from their W + 2-jet analysis [36], namely 15 ± 3(stat.) +5 −3 (sys.) mb and by applying their cuts (see Table 1). Doing so, they assumed that a significant, but non-dominant, DPS yield was to be expected and that their yield could not simply be compared to the SPS predictions above. In order to extract the genuine SPS yield, the only one containing novel information on the quarkonium production mechanisms, they subtracted what they believed to be the DPS yield with σ eff 15 mb.
To be more precise, the data-theory comparison of [4] was done with the ratio of the J/ψ+Z yield over that for Z in order to reduce the uncertainty from the Z observation. The comparison with theory required them to evaluate σ(Z) under their kinematical conditions (see Table 1), where they used σ(pp → Z → e + e − ) = 533.4 pb. After the selection of the prompt J/ψ's in the sample and the subtraction of the expected DPS yield (see above), they obtained : = (45 ± 13 stat ± 6 syst ± 10 DPSsub ) × 10 −7 , whereas the CS based predictions are around (1 − 5) × 10 −7 and the most optimistic NRQCD-based predictions (with CS and CO contributions) reaches 9 × 10 −7 . As announced in our introduction, the experimental results are indeed at least five times larger (corresponding to 2-σ deviation) than the largest available theoretical predictions on the market. Larger LDMEs for the 3 S [8] 1 transition, in particular, could reduce the gap a little but at the cost of a discrepancy with the -much more precise-single-J/ψ data. In addition, introducing a nonzero 3 P [8] J contribution would probably interfere negatively with the dominant 3 S [8] 1 one. Instead of playing further with these parameters, we have thus decided to analyse the process in the CEM up to NLO accuracy, which provides an upper theory value above which any other similar evaluations would probably be unrealistic, except from a totally overlooked partonic reaction.
same, although the central scale we take now is µ 0 = M Z [16] instead of the transverse mass of J/ψ [15].
To evaluate prompt R J/ψ+Z , we used σ(pp → Z → e + e − ) = 427 pb 8 . Our results for the NLO CEM SPS contributions are shown in Table 2. They are on the order of 8 × 10 −7 and similar to the most optimistic NRQCD ones but still far for the ATLAS experimental value. It seems that the solution to the puzzle is not to be found in this direction. The K factor for the hard part we found is 2.   Table 2: Comparison for the cross-section ratio prompt R J/ψ+Z between the CEM predictions for the SPS yields, our DPS extraction and the experimental results at √ s = 8 TeV. The theoretical uncertainty for the NLO SPS is from the renormalisation and factorisation scales. [All quantities are in units of 10 −7 ]. 8 We have employed MadGraph5_aMC@NLO [62] to calculate pp → Z → e + e − by taking into account the NLO QCD corrections. The spin-correlated decay Z → e + e − is done with the help of the module MadSpin [70]. We have also matched the NLO calculation of pp → Z → e + e − to the parton showers provided by Pythia8.1 [71] via the MC@NLO method [72]. The cross section σ(Z)Br(Z → e + e − ) in the Z selection condition of Table. 1 is 427 pb at the fixed order NLO without spin correlations in the Z decay. It varies by 20% when the spin correlation and the parton-shower effects are taken into account, i.e. 505 − 520 pb. Apart from questionning the reliability of the ATLAS measurements or from invoking new physics contributions, the only other possibility left to solve the gap is to question the size of the DPS yield subtracted by ATLAS. As just said, they opted for a σ eff close to 15 mb based on their W + 2-jet analysis. Yet, recent quarkonium-related analysis [24,25,29] have pointed at values smaller than 10 mb, which would in turn induce a larger DPS yield to be subtracted.
We have therefore decided to simply fit σ eff to the total "inclusive" ATLAS J/ψ + Z yield along with the NLO CEM one for the SPS contribution and we obtained σ eff = 4.7 mb (see Table 2 and Fig. 3). Although, in principle, we could have followed the procedure we used for J/ψ +J/ψ production in [29] to estimate the DPS contribution, we have decided to rely on the ATLAS computation, apart from the normalisation obviously. In other words, our DPS yield is 3 times larger than the one employed by ATLAS, i.e. corresponding to a R of 54 × 10 −7 vs 18 × 10 −7 . In order to have an estimation of the DPS yield for the ATLAS and CMS fiducial regions, which are not given in [4], we have used a Crystal Ball function fit (see [29] for technicalities) to calculate the acceptance of J/ψ production from DPS. This allows us to compare our results to the ATLAS measurement in their fiducial region and to predict the total yields within the CMS fiducial region.
Using these numbers, we can also derive an upper limit on σ eff (corresponding to the smallest acceptable DPS yield) by subtracting the 1-σ higher value of the NLO CEM yield from the 1-σ lower value of the ATLAS measurements. This gives 7.1 mb. Even combining both these extreme limits, the ATLAS yield points at a somewhat large DPS. If now one assumes no SPS at all, one can extract a lower value for σ eff as low as 3.2 mb.

The P J/ψ T dependence
Obviously, a DPS yield three times larger (with σ eff = 4.7 mb) would significantly change the azimuthal distribution from which the ATLAS collaboration concluded for a significant, but non-dominant DPS yield. Under our assumption, it becomes now the dominant contribution, about five times as large as the SPS one. Before showing that it does not create any tension with this azimuthal distribution, we however need to discuss the cross section of a function of P J/ψ T , for a reason which will become clear in the next section.  As what regards the SPS CEM cross section, the P J/ψ T spectrum is computed exactly as above with just one less integral, whereas the P T dependence of the DPS yield is exactly the one of ATLAS from Table 5 of [4] with a simple rescaling of the normalisation due to the change in σ eff . The resulting dependence (in form of the ratio R) is shown on Fig. 4. The sum of the SPS and DPS yield in gray gives a reasonable account of the ATLAS yield, with a slight gap opening at large P T . The fact that the agreement is good at low P T is however just a consequence of the fit of σ eff . Such a distribution is rather a consistency check than a test.
However, this helps to clearly illustrate how the low P J/ψ T yield is heavily dominated by the DPS contributions (as is the total yield) and that the high P J/ψ T yield is exclusively from SPS contributions. This a key observation for our discussion in the next section.

The azimuthal distribution
This difference of the P J/ψ T spectrum indeed has an unexpected consequence on the interpretation of the azimuthal distribution of the ATLAS events since it was done with event S (8.5, 10) 10.6 (10,14) 21.0 (14,18) 6.2 (18,30) 9.8 (30, 100) 8.6 (8.5, 100) 56.1 counts, without efficiency correction 9 . It happens that the ATLAS efficiency is much higher for the last bin in P J/ψ T than for the first bin, up to 3 times in fact. This is visible from the statistical uncertainties in Fig. 4 which remain more or less constant whereas the yield is admittedly much smaller in the last bin. In turn, the events used for the ∆φ distribution results from a biased sample which is strongly enriched in high P J/ψ T events. As we discussed in the previous section, high P J/ψ T events are essentially of SPS origin thus mostly populating the ∆φ ∼ π region. Our claim is that the fact that the peak is visible is only due to that and not to a large SPS yield in general. In other words, such a distribution (unless corrected for efficiency in the future) cannot be used to discuss the ratio DPS/SPS integrated in P J/ψ T , nor to reliably extract the DPS yield or σ eff . The lower limit obtained by ATLAS DPS by assuming that the first bin in ∆φ ∼ π was fully from DPS events, 5.3 mb, is slightly different than ours since it ignores the information from data at ∆φ away from 0. Without proper information on the DPS/SPS ratio in each bin (which can be process dependent 10 ), this remains an assumption to consider that the ∆φ ∼ 0 bin is entirely fed by DPS events.
A way to check our hypothesis is to approximately fold our DPS and SPS P J/ψ T spectra with the ATLAS efficiency, to plot and sum the DPS and SPS ∆φ distributions. Since the ATLAS efficiency as a function of P J/ψ T is not publically released, we have used the following simple makeshift, i.e. to derive its P J/ψ T dependence from the corrected yield dependence (from Fig. 4) and the raw number of events in each bins derived from the statistical uncertainty quoted by ATLAS (see also Fig. 4). Such an estimation however also requires the knowledge of the background size, B, which was not given, but which we assumed to be suppressed with respect to S (i.e. the true J/ψ + Z events) as B/S ∝ (P Approximately knowing S (see Table 3) and the theoretical SPS/DPS fraction in each P J/ψ T bin, we were able to stack our theoretical events SPS and DPS events, bin by bin in P J/ψ T , in the ∆φ plot with their specific ∆φ distributions (flat for DPS, peaked at ∆φ π 9 The efficieny was however checked to be constant in ∆φ [73]. 10 For process like J/ψ-pair production, it is well know that it very much depends on the PT cuts [26,31,74]. for the SPS following our NLO CEM computation). The resulting distribution is shown on Fig. 5. They look essentially the same at LO and NLO and demonstrates that increasing the DPS yield by a factor of 3 does not create any tension with the observed ATLAS event ∆φ distribution if the efficiency corrections are approximately accounted for in the theory evaluations. We are of course eager for an updated analysis by ATLAS to avoid such a complicated procedure to compare the theory to the experimental ∆φ distribution.

4.5
The P Z T dependence Before concluding, we also give our predictions for the P Z T dependence of the yield. For the SPS yield, it follows from our NLO CEM evaluation. For the DPS yield, it follows from the simple observation that the absence of correlation between the J/ψ and the Z kinematics in DPS events 11 allows us to state that the normalized distribution 1 σ dσ dP Z T should be identical for both the single-Z production and the (DPS) J/ψ + Z production. In fact, this could be used to check the dominance of the DPS yield in some part of the phase space. The normalised distribution of 1 σ dσ dP Z T for inclusive Z production can accurately be predicted with the help of MadGraph5_aMC@NLO, which can be evaluated at NLO and matched to parton-shower routines. Fig. 6 shows the P Z T dependence of R for both the ATLAS and CMS acceptances. One notices that the NLO and LO SPS spectra are different for large P Z T values. This shows that the leading topologies to generate high-P T Z bosons are not present at LO. One also observes that the DPS yield dependence is softer than the NLO SPS one and similar to the LO SPS one.

Conclusions
In this paper, we have performed a theoretical re-analysis of the associated production of a prompt J/ψ with a Z boson at the LHC in view of the thought-provoking results of ATLAS [4] at 8 TeV LHC. Despite its cross section on the order of a few femtobarns, this process has the potential to tell us much about the quarkonium-production mechanisms and on the DPS physics in the case of a highly asymmetric system. For the first time, we have performed a NLO calculation of the hadroproduction of a J/ψ + a recoiling parton in the CEM with the help of MadGraph5_aMC@NLO and fixed the CEM non-perturbative parameter P NLO,prompt J/ψ by fitting the most recent and precise ATLAS data. With the extracted P NLO,prompt J/ψ , we have presented the first theoretical calculation of the prompt J/ψ + Z SPS production at NLO in the CEM which directly follows from the quark-hadron-duality principle in the context of quarkonium production. We do believe that our calculation can be considered as a conservative upper value of the SPS yield. In other words, any NRQCD evaluation with a larger yield would likely be based on parameters incompatible with the existing single-J/ψ data. This allows us to state that, in the ATLAS and CMS fiducial regions 12 , SPS contribution dominate the "high" P J/ψ T region, where the DPS contributions tend to dominate the "low" P J/ψ T one. Complementary measurements by LHCb at P J/ψ T as low as a few GeV would therefore be extremely important whereas future measurements at higher P T , with larger luminosities, would normally test the quarkonium-production mechanisms.
A thorough comparison with the ATLAS data has then been made by taking into account both SPS and DPS contributions. We argue that the DPS yield subtracted by ATLAS with σ eff 15 mb was underestimated and thus the expected impact of SPS overestimated. In fact, their fitted lower value, 5.2 mb, assuming the absence of SPS contributions at ∆φ → 0 and ignoring their ∆φ event distribution, is well compatible with our argument. With our NLO SPS results, we are indeed able to fit σ eff 4.7 mb from the P T -integrated total yield and compute the other differential distributions (e.g. azimuthal distribution and transverse momentum distributions). Considering the CEM as indeed an upper theory limit for the SPS yield, we are also able to derive an upper value for σ eff corresponding to the smallest possible DPS yield in agreement with the data. On the other hand, assuming a negligible SPS contribution to the P T -integrated yield 13 , we obtain a virtual lower σ eff limit, which is more conservative than the one derived by ATLAS. The σ eff range obtained likewise is in good agreement with other quarkoniumrelated extractions (see Fig. 7) and is visibly lower than the ones extracted from jetrelated observables, pointing at a possible process dependence of σ eff . Indeed, one should keep in mind that all these extractions rely on the implicit factorisation of the pocket formula (σ DPS A+B ∝ σ A × σ B ) and there does not exist proofs of such a factorisation. In fact, factorisation-breaking effects have been discussed in a number of recent studies (see e.g. [75][76][77]).  [24,25,29,36,[78][79][80][81][82][83].
These results are quite encouraging as the agreement between our theoretical results and the ATLAS data is acceptable, given the reduced number of events, and only a small discrepancy is present in P J/ψ T distribution. This calls for further analyses by experiments. In this paper, we also provide our predictions for the ongoing CMS measurement. We are also hopeful that our theoretical evaluation will motivate further theory updates at NLO to 13 hypothesis however disfavoured by differential distributions. substantiate our claim that the DPS-subtracted yield of ATLAS can indeed be accounted by known mechanisms.