Measurement of the production cross-section of $J/\psi$ and $\psi(2$S$)$ mesons in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

Measurements of the differential production cross-sections of prompt and non-prompt $J/\psi$ and $\psi(2$S$)$ mesons with transverse momenta between 8 and 360 GeV and rapidity in the range $|y|<2$ are reported. Furthermore, measurements of the non-prompt fractions of $J/\psi$ and $\psi(2$S$)$, and the prompt and non-prompt $\psi(2$S$)$-to-$J/\psi$ production ratios, are presented. The analysis is performed using 140 fb$^{-1}$ of $\sqrt{s}=13$ TeV $pp$ collision data recorded by the ATLAS detector at the LHC during the years 2015-2018.


Introduction
Studies involving heavy quarkonia provide a unique insight into the nature of quantum chromodynamics (QCD) near the boundary of the perturbative and non-perturbative regimes.However, despite a long history of research, quarkonium production in hadronic collisions still presents significant challenges to both theory and experiment.
In high-energy hadronic collisions, charmonium states can be produced either from short-lived QCD sources (referred to as prompt production) or from long-lived sources -decays of beauty hadrons (non-prompt production).These two sources can be distinguished experimentally by measuring the distance between the production and decay vertices of the charmonium state.Feed-down decays of higher charmonium states contribute to the production of / mesons for both of these sources, and should be taken into account when comparing with theoretical predictions (there is no significant feed-down contribution to (2S) production).While calculations within the framework of perturbative QCD (see e.g.Refs.[1,2]) have been reasonably successful in describing the non-prompt contributions, a satisfactory understanding of the prompt production mechanisms is still to be achieved.
Methods developed within the non-relativistic QCD (NRQCD) approach provide a framework for describing quarkonium production processes, leading to a variety of models differing in their accuracy and predictive power.In particular, Ref. [3] introduced a number of phenomenological parameters -long-distance matrix elements (LDMEs) -which can be extracted from fits to the experimental data, and are expected to describe the cross-sections and differential spectra for several sets of data reasonably well [4][5][6][7].However, various attempts to build a universal library of LDMEs to be used to describe a wider range of measurements such as the polarisation of quarkonia [8][9][10][11], their associated production [12,13] or the production of quarkonium in a wider range of processes (e.g.photo-and electro-production) have not been particularly successful [14][15][16][17][18].An alternative approach to the description of the hadronisation process of heavy quarkonia is offered by the Colour Evaporation Model (CEM) [19][20][21], which offers a simpler framework with fewer parameters, but has its own problems in describing the data [22].A combination of ATLAS results with cross-section and polarisation measurements from CMS [7,23,24], LHCb [6,[25][26][27][28][29] and ALICE [11,[30][31][32][33][34][35][36] now includes a variety of charmonium production characteristics in a wide kinematic range, thus providing a wealth of information for a new generation of theoretical models.
One way to add qualitatively new information is to extend the kinematic range of quarkonium production measurements.ATLAS has previously measured the inclusive differential cross-section for / production in   collisions at √  = 7 and 8 TeV [4], as well as the differential cross-sections for the production of   states [37], and for (2S) production [5].In most of these measurements, ATLAS exploited a dimuon trigger with a muon transverse momentum ( T ) threshold of 4 GeV, with the high- T reach limited mainly by the dimuon trigger's performance to about 100 GeV: at higher  T values the angular resolution of the muon trigger system is not sufficient to separate the two almost collinear muons.This paper describes a measurement of / ((2S)) meson production via decay in the dimuon channel, at √  = 13 TeV for meson transverse momenta of 8-360 GeV (8-140 GeV), which is a much broader range than in previous measurements.This was made possible by the use of two different triggers.Production of / and (2S) at low  T , between 8 and 60 GeV, is measured using a dimuon trigger requiring a pair of muons to each pass a  T threshold of 4 GeV, while at high  T a single-muon trigger with a  T threshold of 50 GeV was used.This allowed measurements to be performed for transverse momenta as high as 360 GeV for / and 140 GeV for (2S).The measurements include the double-differential cross-sections for production of the two vector charmonium states (separately for the prompt and non-prompt production mechanisms), the non-prompt fraction for each state, and the prompt and non-prompt (2S)-to-/ production ratios.
The paper is organised as follows.A brief description of the ATLAS detector is given in Section 2. The event selection and the analysis strategy are explained in Section 3, followed by a description of the systematic uncertainties affecting the measurement in Section 4. Results and comparisons with theoretical calculations are presented in Section 5, followed by a summary in Section 6.

The ATLAS detector
The ATLAS experiment [38] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4 coverage in solid angle. 1 It consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer.The inner tracking detector covers the pseudorapidity range || < 2.5.It consists of silicon pixel, silicon microstrip, and transition radiation tracking detectors.Metal/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity.A steel/scintillator-tile hadron calorimeter covers the central pseudorapidity range (|| < 1.7).The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to || = 4.9.The muon spectrometer surrounds the calorimeters and is based on three large superconducting air-core toroidal magnets with eight coils each.The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector.The muon spectrometer has a system of precision chambers for tracking and fast detectors for triggering.A two-level trigger system is used to select events.The first-level trigger is implemented in hardware and uses a subset of the detector information to accept events at a rate below 100 kHz.This is followed by a software-based trigger that reduces the accepted event rate to 1 kHz on average depending on the data-taking conditions.An extensive software suite [39] is used in data simulation, in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

Event selection
Data for this analysis were taken during the LHC proton-proton collision runs at √  = 13 TeV in the years 2015 to 2018.For lower values of the transverse momentum  T of the dimuon system, between 8 and 60 GeV, a dimuon trigger was used, requiring a pair of muons to each pass a  T threshold of 4 GeV.This trigger ran unprescaled during 2015 data-taking, collecting an integrated luminosity of 2.6 fb −1 .In the high  T range between 60 and 360 GeV, a single-muon trigger with a  T threshold of 50 GeV was used, unprescaled throughout the full Run 2 data-taking, providing a total integrated luminosity of 140 fb −1 .The selected events were required to contain a pair of oppositely charged muons of high quality (using the tight identification requirements defined in Ref. [40]), with  T > 4 GeV and || < 2.4.In the low  T range the two muons were required to match the two trigger objects of the dimuon trigger, while in the high  T range at least one of the muons was required to have  T > 52.5 GeV and match the trigger object.The two 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe.The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards.Cylindrical coordinates (, ) are used in the transverse plane,  being the azimuthal angle around the -axis.
The pseudorapidity is defined in terms of the polar angle  as  = − ln tan(/2).Angular distance is measured in units of ID tracks attributed to the muons were fitted to a common vertex, and the dimuon invariant mass   was required to satisfy 2.6 <   < 4.2 GeV.The transverse distance    between the primary vertex and the dimuon vertex was used to calculate the meson's pseudo-proper decay time , where  T is the reconstructed transverse momentum of the dimuon system, and  is the speed of light.The primary vertex is chosen as the reconstructed collision interaction vertex whose  coordinate is nearest to the point of closest approach of the dimuon system's trajectory to the beam axis.If an event has more than one selected dimuon candidate, all candidates are retained and treated independently.

Cross-section determination
The phase space of the measurement is divided into 34 intervals in dimuon  T covering the range from 8 to 360 GeV, and 3 intervals in absolute rapidity 2 || with boundaries at 0, 0.75, 1.5 and 2.0, thus producing 102 analysis bins overall.In each (  T , ) bin, a two-dimensional unbinned maximum-likelihood fit to the distribution of dimuon candidates in invariant mass   and pseudo-proper decay time  of the  meson is performed to obtain the raw yields  P,NP  for prompt (P) and non-prompt (NP)  mesons, where  = /, (2).The raw yields are then corrected to account for the geometrical acceptance A (), the trigger and reconstruction efficiencies  trig and  reco , and the trigger and reconstruction correction scale factors  trigSF and  recoSF , averaged over that bin.Several low  T bins are divided into narrower sub-bins to obtain raw yields in finer granularity, which are then corrected and summed to give the final yield in the corresponding analysis bin.This procedure helps to reduce measurement biases due to modelling assumptions in the regions of phase space with large statistical power.
The prompt (P) and non-prompt (NP) double-differential production cross-sections for  = /, (2) in each analysis bin are calculated as where Δ T and Δ are bin widths in  T and rapidity, and ∫ Ld is the corresponding integrated luminosity.Bin migration effects are discussed in Section 3.4.The acceptance A () is defined as the probability that a  state with (true) momentum within an analysis bin survives the following acceptance selections imposed on the two muons (assuming  T (1) >  T (2)) in the two   T ranges: The acceptance calculation is performed using Monte Carlo (MC) generator-level kinematic variables, with resolution effects taken into account at the efficiency correction stage.An isotropic angular distribution of muons in the  decay frame is assumed.Since the spin alignment of the  states may affect the acceptance, a number of non-isotropic spin-alignment scenarios are used to calculate correction factors for the measured cross-sections (see Appendix).For a given spin-alignment scenario, systematic uncertainties due to the acceptance calculation are small (see Section 4).Changing to a different spin-alignment scenario, however, can lead to noticeable changes in the cross-sections and other measured quantities, including variations as a function of  T (see Section 5).
The reconstruction efficiency  reco and trigger efficiency  trig are calculated using samples of fully simulated / and (2S) events, including appropriate trigger information.Correction scale factors  recoSF and  trigSF account for the differences between simulated and real data.See Section 3.4 for more details.
The non-prompt fractions for  = /, (2S) are defined as Finally, (2S)-to-/ production ratios are defined separately for the prompt and non-prompt production mechanisms as In calculating these quantities, the event yields, efficiencies, and acceptance corrections are used in accord with Eq. ( 1); uncertainties in the fraction and ratio measurements partially cancel out.

Fit model
The fit model's probability distribution function  (  , ) contains seven terms, with fractions   , describing four signal contributions and three types of background.Terms  = 1, 2 describe prompt and non-prompt / signal respectively; terms  = 3, 4 correspond to prompt and non-prompt (2S) signal.Term  = 5 describes the prompt background, where non-resonant dimuons are produced at the primary vertex (e.g.Drell-Yan pairs).Term  = 6 describes single-sided non-prompt background, mainly for dimuon continuum events where the two muons originate from the (cascade) decay of a single -hadron, while term  = 7 describes the double-sided part of the non-prompt continuum, where the two muons originate from different -hadrons, yielding a secondary vertex which may appear on either side of the beamline.For  = 2-7, each term is factorised into a function   (  ) of dimuon mass   and a function ℎ  () of pseudo-proper decay time , where the latter is convolved with a decay time resolution function (): The term with  = 1, which describes the prompt / peak, has a similar structure but allows for some correlations between   and , as described below.
The decay time resolution function () is parameterised as a sum of three Gaussian functions,  A ,  B and  C , with the relative weights of the first two,  A and  B , treated as free parameters: Based on MC studies with fully simulated signal samples, the means of the three Gaussian functions are fixed to zero, and the widths are linked by  B = 2 A and  C = 4 A , where  A = 0.04 ps is fixed to the smallest value found in test fits.
The parameterisations of the functions   (  ) and ℎ  () are summarised in Table 1.The mass lineshapes of the / and (2S) peaks,   (  ) for  = 1 to 4, are parameterised as weighted sums of two Gaussian functions and a Crystal Ball function [41], which are the same for the prompt and non-prompt components.
Based on MC studies, the weights are common to / and (2S), while the ratios of the peak positions and the widths are fixed to the ratio  of the masses [42] of the two states.Parameters of the Crystal Ball function were kept the same for both prompt and non-prompt / and also (2S), as verified from the MC studies.The lifetime distributions ℎ 1 () and ℎ 3 () of prompt / and (2S), respectively, are parameterised as delta functions, while for non-prompt (2S), a single-sided exponential function is used for ℎ 4 ().Since the non-prompt / sample is larger and has a wider observed  range, its decay time distribution ℎ 2 () is described by a superposition of two single-sided exponential functions with slopes related by  1 =  2 , where the constant  = 1.4 is obtained from test fits using real data.All these exponential functions are convolved with the resolution function ().
In the  = 1 term, describing the prompt / peak, the product of the narrowest Gaussian term in  1 (  ) and the narrowest Gaussian term in () was replaced by a bivariate Gaussian function in   and  with a correlation coefficient  = 0.3, to take into account the observed correlation between the measured values of these quantities.The effect of this correlation was found to be negligible for other terms.
Parameterisations for the background terms are selected using both the experience gained from similar analyses at lower energies [4] and physics considerations.In the prompt background term,  = 5, the mass distribution is modelled by a second-order polynomial, while the non-prompt mass distributions for  = 6 and 7 are parameterised as exponential functions, with independent parameters.The decay time distribution is a delta function for the prompt term  = 5, a single-sided exponential function for the main non-prompt term  = 6, and a symmetric double-sided exponential function for the last term,  = 7.Each of these is also convolved with the decay time resolution function ().
Table 1: Parameterisation of the fit model.Here ,  ,  and  denote Gaussian, Crystal Ball, exponential and second-order polynomial functions, respectively, with different indices corresponding to different parameters.The parameterisation of the  = 1 term is modified as described in the text and shown in Eq. ( 4).
i Type P/NP   () ℎ  () Fits are performed in each (sub-)bin using an unbinned maximum-likelihood method, in the dimuon mass range from 2.7 to 4.1 GeV and the decay time range between −1 and 11 ps.Twenty of the 29 parameters are determined from the fit, with the rest fixed to predetermined values obtained from test fits using samples of MC simulated / and (2S) signal events.Uncertainties due to variations of the fit model and assumptions about the fixed parameters are estimated during the studies of systematic uncertainties described in Section 4. In particular, it is found that a reliable determination of the cross-sections for prompt and non-prompt production of the (2S) meson (with significance better than 5) is only possible up to  T = 140 GeV, mainly due to poorer mass resolution and a lower signal-to-background ratio at higher  T .For  T > 140 GeV the yield of (2S) was fixed to a constant fraction (0.07) of the yield of /, by extrapolating the results from the fits at lower transverse momenta.
Figure 1 shows the mass and pseudo-proper decay time projections of the fits in several sample bins, together with the associated pull distributions.The quality of the fits, assessed by calculating a two-dimensional  2 value, is found to be good in all (sub-)bins.The main parameters determined from the fits are the prompt and non-prompt yields of / and (2S) states.The cross-sections, non-prompt fractions and production ratios were then calculated using Eqs.( 1)- (3).taking into account correlations between fit parameters.The results for all measured quantities are presented in Section 5.

Efficiency corrections
As shown in Eq. ( 1), the yields obtained from two-dimensional maximum-likelihood fits in each (sub-)bin are subject to acceptance corrections, followed by the corrections for reconstruction and trigger efficiencies obtained from / and (2S) MC simulations, and by the correction scale factors that account for differences between data and MC simulation.
MC samples used for efficiency determinations were produced either by the Pythia 8 generator [43] with the A14 set of tuned parameters [44], or by a custom particle gun generator producing single  states with a given distribution, followed by their decay into a dimuon final state.The generated events were passed through the full ATLAS detector simulation [45] based on Geant4 [46], and were reconstructed using the same software as the real data.The reconstruction efficiency  reco in each analysis (sub-)bin is defined as the ratio where  reco is the reconstructed yield within the bin boundaries defined in terms of reconstructed variables, with fiducial cuts applied to reconstructed variables, while  true is the true (generated) yield within the bin boundaries defined in terms of true variable values, with fiducial cuts applied to the true variable values.This definition takes into account detector resolution smearing of the kinematic variables used to define the fiducial cuts and bin boundaries, and hence includes bin migration effects between neighbouring bins.
The trigger efficiency  trig , also obtained using MC simulations, is defined fully in terms of reconstructed variables as where  trig is the number of triggered events among the reconstructed events.Since this measurement used two different triggers, the trigger efficiency was calculated accordingly.Further correction scale factors,  recoSF and  trigSF , are applied to account for any differences between data and MC events at reconstruction level and trigger level (see Eq. ( 1)).These are evaluated in dedicated tag-and-probe studies with auxiliary triggers, using a mixture of / →  +  − and  →  +  − decays (see Ref. [40] for details).

Systematic uncertainties
Systematic effects from a variety of sources were studied, and appropriate corrections and uncertainties were assigned to all measured quantities.The systematic uncertainties can be broadly grouped into those related to reconstruction, trigger, and acceptance corrections, and those related to the fit model.
In order to assess the systematic uncertainties related to the fit model, the fit model was varied in several ways.As mentioned in Section 3.3, in the nominal fits some of the parameters describing the lineshapes of the signal peaks in the mass and lifetime domains were fixed to the values obtained from signal MC studies.These include the values of parameters  and  of the Crystal Ball function, the constant  linking the positions of the / and (2S) mass peaks, and the factor  relating the slopes of the two exponential functions describing the distribution of / decay times, .These parameters were allowed to float, one at a time, and the fits were repeated.Some variations covered changes in the decay time resolution parameterisation, with widths of the three Gaussian functions changed independently.In other variations, alternative parameterisations were chosen for the mass dependence of the background terms, and the fits were repeated.In one of the variations the correlation parameter  from Eq. 4 was set to 0. Finally, in the highest  T bins,  T > 140 GeV, the fixed fraction 0.07, relating the (2S) and / yields, was varied by 0.01 to cover the observed range at lower  T values, and the fits were run again.After each rerun, changes in the measured yields were recorded.The outcome of this process, in each analysis bin, was a number of measurements of each yield, scattered around the result of the nominal fit for that yield.It was assumed that the probability distribution of these variations around the nominal value was uniform between the smallest and largest measured values, and, therefore, the corresponding systematic uncertainty was evaluated as the standard deviation of this uniform distribution.
Acceptance-related systematic uncertainties are governed by the size of the samples used to generate the corresponding acceptance maps.The chosen size ensured that these uncertainties are small relative to the total systematic uncertainties.Changes in acceptance due to different spin-alignment hypotheses were treated separately (see Appendix).
Systematic uncertainties due to trigger and reconstruction efficiency corrections include the uncertainties on the scale factors [40] as well as those due the sizes and shapes of the MC samples.These were added in quadrature.
The total systematic uncertainty in each bin is calculated by summing in quadrature the uncertainties from the above-mentioned sources.The total uncertainty is calculated as the sum in quadrature of the total systematic and statistical uncertainties.
An additional systematic uncertainty comes from the determination of the integrated luminosity.In the  T < 60 GeV region, where only 2015 data contributes to this measurement, the integrated luminosity has a 1.13% uncertainty [47], while for  T ≥ 60 GeV the uncertainty in the combined 2015-2018 integrated luminosity is 0.83%, as obtained using the LUCID-2 detector [48].
As an illustration, in Figure 2 the total, statistical, and systematic uncertainties, together with the main individual contributions to the systematic uncertainty, are shown for the central rapidity slice for the differential cross-sections of prompt / and non-prompt (2S) mesons, as well as for the non-prompt fractions of / and (2S) mesons, as functions of  T .For the cross-sections, apart from a few high  T bins, the uncertainties are largely systematic, while for the non-prompt fractions and (2S)-to-/ production ratios, statistical errors dominate in many bins because the systematic uncertainties partially cancel out.Error levels are similar between prompt and non-prompt cross sections.Comparable results are obtained in other rapidity slices.

Results
The measured double-differential cross-sections for prompt and non-prompt / production in the nominal isotropic spin-alignment scenario are presented in Figures 3(a  Figure 3: Differential cross-sections for (a) prompt and (b) non-prompt production of / mesons.For visual clarity, a scaling factor of 1, 10, or 100 is applied to the rapidity slices 0.00 ≤ || < 0.75, 0.75 ≤ || < 1.5, and 1.5 ≤ || < 2.0, respectively.For each data point, the horizontal bar spans the  T range covered by that bin, with the horizontal position of each point representing the mean  T in that bin.The vertical uncertainty range (obscured by the marker for some values) combines both the statistical (the inner bar) and total uncertainty.Uncertainties related to spin alignment or integrated luminosity are not included.Data up to 60 GeV were taken with a dimuon trigger with integrated luminosity 2.6 fb −1 ; data above 60 GeV were taken with a single-muon trigger with integrated luminosity 140 fb −1 .presented in Figures 5(a) and 5(b).Finally, the (2S)-to-/ production ratios are presented in Figures 6(a) and 6(b) for the prompt and non-prompt production mechanisms, respectively.
While the non-prompt fractions shown in Figure 5 increase steadily with  T up to about 100 GeV, they are almost constant for both / and (2S) in the high  T range, which suggests similar  T -dependences for the prompt and non-prompt differential cross-sections at very high transverse momenta.

Acceptance and spin alignment corrections
The transition between the low- T dimuon trigger and the high- T single-muon trigger at  T = 60 GeV presents a particular challenge because of the sharp change in event kinematics.The corresponding changes in the acceptance and efficiency correction factors are significant and could lead to discontinuities in the measured distributions.Figure 4: Differential cross-sections for (a) prompt and (b) non-prompt production of (2S) mesons.For visual clarity, a scaling factor of 1, 10, or 100 is applied to the rapidity slices 0.00 ≤ || < 0.75, 0.75 ≤ || < 1.5, and 1.5 ≤ || < 2.0, respectively.For each data point, the horizontal bar spans the  T range covered by that bin, with the horizontal position of each point representing the mean  T in that bin.The vertical uncertainty range (obscured by the marker for some values) combines both the statistical (the inner bar) and total uncertainty.Uncertainties related to spin alignment or integrated luminosity are not included.Data up to 60 GeV were taken with a dimuon trigger with integrated luminosity 2.6 fb −1 ; data above 60 GeV were taken with a single-muon trigger with integrated luminosity 140 fb −1 .
Since the spin alignment of  states may be different for the prompt and non-prompt production mechanisms, additional correction factors may be needed for all measured distributions.In order to quantitatively assess the possible impact of  spin alignment, correction factors are calculated for a variety of scenarios.It was found that the dependence on the polar angle  in the helicity frame of the  state causes the largest variation, so the angular dependence of  →  +  − decays is assumed to be ∝ (1 +   cos 2 ).The correction factors are shown in Figures 7(a) and 7(b) for the differential cross-sections and non-prompt fractions respectively, where the values   = ±0.20 are chosen to reflect the approximate level of experimental knowledge [7,36,49] and theoretical understanding [50][51][52] of this parameter.The correction factors are shown for prompt / in the central rapidity range, but were found to be essentially the same for / and (2S), for the prompt and non-prompt production mechanisms, and also for the three rapidity regions.
The potential bias due to the spin-alignment assumption is especially noticeable at the  T = 60 GeV transition, and indeed a step can be seen at this point in the / non-prompt fraction in Figure 5(a), reflecting a possible issue in the spin-alignment modelling.Correction factors for some other values of   are presented in the Appendix.These can be used for further studies, if more precise spin-alignment data Figure 5: Non-prompt production fraction of (a) / and (b) (2S) mesons.For visual clarity, a vertical shift of 0, 0.2, or 0.4 is applied to the rapidity slices 0.00 ≤ || < 0.75, 0.75 ≤ || < 1.5, and 1.5 ≤ || < 2.0, respectively.For each data point, the horizontal bar spans the  T range covered by that bin, with the horizontal position of each point representing the mean  T in that bin.The vertical uncertainty range (obscured by the marker for some values) combines both the statistical (the inner bar) and total uncertainty.Uncertainties related to spin alignment or integrated luminosity are not included.Data up to 60 GeV were taken with a dimuon trigger with integrated luminosity 2.6 fb −1 ; data above 60 GeV were taken with a single-muon trigger with integrated luminosity 140 fb −1 .and/or improved modelling become available in the future.

Theory comparison: prompt production
Model calculations of prompt production of charmonium are usually based on perturbative QCD for the production of the  c pair, and differ in the mechanism of hadronisation and formation of the bound state with specific quantum numbers.
The predictions of a model using the non-relativistic QCD approach to charmonium production crosssections at next-to-leading order (NLO NRQCD) [53], using predetermined LDMEs [54,55], are shown in comparison with our measurements of the / and (2S) production cross-sections in the top panels of Figures 8(a) and 8(b) respectively.The predictions of the model largely overlap with the data points within the theoretical uncertainties, which include variations of the renormalisation, factorisation and NRQCD scales.However, the predictions seem to overestimate the cross-sections at high  T .
One generalisation of the NRQCD approach is a model which aims to improve the description by taking into account the transverse degrees of freedom of the initial gluons in the colliding protons (  -factorisation model) [56,57].The predictions of this model, obtained with the PEGASUS event generator [58] using the LDMEs determined in Ref. [59], are compared with our measured results as shown in the middle panels of Figures 8(a  Figure 6: The (2S)-to-/ production ratio for the (a) prompt and (b) non-prompt production mechanisms.For visual clarity, a vertical shift of 0, 0.02, or 0.04 is applied to the rapidity slices 0.00 ≤ || < 0.75, 0.75 ≤ || < 1.5, and 1.5 ≤ || < 2.0, respectively.For each data point, the horizontal bar spans the  T range covered by that bin, with the horizontal position of each point representing the mean  T in that bin.The vertical uncertainty range (hidden by the marker for some values) combines both the statistical (the inner bar) and total uncertainty.Uncertainties related to spin alignment or integrated luminosity are not included.Data up to 60 GeV were taken with a dimuon trigger with integrated luminosity 2.6 fb −1 ; data above 60 GeV were taken with a single-muon trigger with integrated luminosity 140 fb −1 .account for variations of the renormalisation scale.The predictions of the model tend to underestimate the cross-sections at low  T , while in the high  T region the range of comparison is limited by the availability of the transverse-momentum-dependent parton distribution function (PDF) of the gluon [60].
A different approach is used by the Improved Colour Evaporation Model (ICEM) [61], which assigns a fixed fraction of the  c production cross-section below the open charm threshold to individual charmonium states.Comparisons of ICEM predictions with the parameter values and their uncertainties previously determined from fits to LHCb data at 7 TeV [10,62] are shown in the bottom panels of Figures 8(a) and 8(b) for / and (2S) respectively.The model seems to predict somewhat harder  T spectra than observed in the data for both / and (2S), and tends to underestimate the cross-section for (2S).

Theory comparison: non-prompt production
Theoretical calculations of non-prompt charmonium production are based on perturbative QCD for the production of a  b quark pair, their hadronisation into a pair of -hadrons, and their subsequent decay into a charmonium state with specific quantum numbers.Predictions from one such model, based on fixed-order-next-to-leading-log (FONLL) QCD calculations [1,2] were obtained using the web-based tool [63] with default parameter values, and are shown in comparison with our measurements in the top panels of Figures 9(a  both the renormalisation scale and the charm quark mass.Agreement is good at lower  T values, but the FONLL model predicts higher cross-sections for / at the high  T end.
Another set of predictions, based on the next-to-leading-order QCD calculation in the general-massvariable-flavour-number scheme (GM-VFNS) [64] are shown in the middle panels of Figures 9(a Parameter values were determined in Ref. [54,65], with uncertainties originating from renormalisation scale dependence.These predictions lead to similar results, but the deviation from data at the highest  T values is somewhat more pronounced, especially in the / case. Finally, the NRQCD model with   -factorisation can also be used to predict the  T distributions of vector charmonia through the non-prompt production mechanisms [58,66] (see bottom panels of Figures 9(a) and 9(b)).Where available, the shapes of  T distributions are reproduced fairly well for both / and (2S), but the limitations of the transverse-momentum-dependent model for the gluon PDF show up at even lower charmonium  T values.Also, in this model the cross-section for (2S) non-prompt production at low  T is somewhat underestimated.
Overall, none of the models considered here is able to describe the data over the whole measured range of transverse momenta.The general trend shown by all theoretical models is a slower-than-observed decrease of the cross-section with  T , which could be related to insufficiently accounting for PDF evolution and/or possible dependence of LDMEs on transverse momentum.In any case, these measurements should help refine theoretical models of hadronic production of quarkonium at the highest available energies and at transverse momenta well beyond 100 GeV.

Summary
This paper describes a measurement of the double-differential production cross-sections of / and (2S) charmonium states in   collisions at √  = 13 TeV, performed through their decays into dimuons and using 140 fb −1 of data collected by the ATLAS detector at the LHC during Run 2. The cross-sections for each of the two states are measured separately for prompt and non-prompt production mechanisms.The non-prompt fractions for each state are also measured, along with the (2S)-to-/ production ratios.The rapidity range of the measurement is || < 2. For (2S) meson the transverse momentum range is 8−140 GeV, while for / state the results cover a much wider transverse momentum range, from 8 GeV to 360 GeV, extending well beyond the range of previous measurements.In the high  T range the results show similar  T -dependences for the prompt and non-prompt differential cross-sections, with the non-prompt fractions being nearly constant for both / and (2S) states.
The results are compared with a number of theoretical predictions, which describe the data with varying degrees of success.The extended  T reach of this measurement provides important fresh input for future tuning of theoretical models.

Figure 1 :
Figure 1: Mass (left) and pseudo-proper decay time (right) projections of the fit result for selected analysis (sub-)bins.The values of  2 /d.o.f. for the corresponding 2-dimensional fits are, from top to bottom: 1.09, 0.89, 1.20 and 0.93.

Figure 2 :
Figure 2: Total, statistical, and systematic uncertainties (in %) as functions of  T for the differential (a) prompt / and (b) non-prompt (2S) cross-sections, and for the non-prompt fractions of (c) / and (d) (2S), in the rapidity slice 0.00 ≤ || < 0.75.The main components of the systematic uncertainties are also shown.
) and 3(b), respectively.The same quantities for (2S) are shown in Figures4(a) and 4(b).The non-prompt production fractions for / and (2S) are ) and 8(b) for / and (2S) respectively, where the theoretical uncertainties only ) and 9(b) for / and (2S) respectively.Here the uncertainties cover variations of

Figure 7 :
Figure 7: Spin-alignment hypothesis correction factors for the / (a) differential cross-section and (b) non-prompt production fraction, for a number of spin-alignment scenarios.The correction factors are approximately the same for / and (2S), for the prompt and non-prompt production mechanisms, and also for the three rapidity regions.The discontinuities at  T = 60 GeV are due to the transition from a low- T dimuon trigger to a high- T single-muon trigger, and the corresponding change in event acceptance.

Figure 8 :Figure 9 :
Figure 8: Ratios of various theoretical predictions (described in the text) to the data points from this measurement, for the prompt production of (a) / and (b) (2S) in the central rapidity region.In each  T bin, the shaded area represents the ratio of the theoretical prediction to the measured value, with the vertical spread showing the uncertainties of the respective model.Error bars on the black circles show fractional uncertainties of this measurement.