First measurement of quarkonium polarization in nuclear collisions at the LHC

The polarization of inclusive J/$\psi$ and $\Upsilon(1{\rm S})$ produced in Pb-Pb collisions at $\sqrt{s_{\rm{NN}}}=5.02$ TeV at the LHC is measured with the ALICE detector. The study is carried out by reconstructing the quarkonium through its decay to muon pairs in the rapidity region $2.5<y<4$ and measuring the polar and azimuthal angular distributions of the muons. The polarization parameters $\lambda_{\theta}$, $\lambda_{\phi}$ and $\lambda_{\theta\phi}$ are measured in the helicity and Collins-Soper reference frames, in the transverse momentum interval $2<p_{\rm T}<10$ GeV/$c$ and $p_{\rm T}<15$ GeV/$c$ for the J/$\psi$ and $\Upsilon(1{\rm S})$, respectively. The polarization parameters for the J/$\psi$ are found to be compatible with zero, within a maximum of about two standard deviations at low $p_{\rm T}$, for both reference frames and over the whole $p_{\rm T}$ range. The values are compared with the corresponding results obtained for pp collisions at $\sqrt{s}=7$ and 8 TeV in a similar kinematic region by the ALICE and LHCb experiments. Although with much larger uncertainties, the polarization parameters for $\Upsilon(1{\rm S})$ production in Pb-Pb collisions are also consistent with zero.


Introduction
Quarkonia, bound states of charm (c) and anticharm (c) or bottom (b) and antibottom (b) quarks, represent an important tool to test our understanding of quantum chromodynamics (QCD), since their production process involves both perturbative and non-perturbative aspects. At high energy, the creation of the heavy quark-antiquark pair is a process that can be described using a perturbative QCD approach, due to the large value of the charm and bottom quark masses (m c ∼ 1.3 GeV/c 2 , m b ∼ 4.2 GeV/c 2 [1]). However, the subsequent formation of the bound state is a non-perturbative process that can be described only by empirical models or effective field theory approaches. Among those, models based on Non-Relativistic QCD (NRQCD) [2] give the most successful description of the production cross section, as measured at high-energy hadron colliders (Tevatron, RHIC, LHC) [3][4][5][6][7][8][9][10][11][12][13][14]. In the NRQCD approach, the non-perturbative aspects are parameterized via long-distance matrix elements (LDME), corresponding to the possible intermediate color, spin and angular momentum states of the evolving quark-antiquark pair. The values of LDMEs need to be fitted on a subset of the available measurements and can be then considered as universal quantities, in the sense that they can be used in the calculation of production cross sections and other observables corresponding, for example, to different collision systems and energies. Other theory approaches, as the Color Singlet Model [15], the Color Evaporation Model [16] and the k T -factorization [17] are also used to describe the quarkonium production process.
Among the various charmonium states, the J/ψ meson, with quantum numbers J PC = 1 −− , was the first to be discovered. It is surely the most studied, also due to the sizeable decay branching ratio to dilepton pairs ((5.961 ± 0.033)% for the µ + µ − channel [1]) that represents an excellent experimental signature. While the J/ψ production cross sections are well reproduced by NRQCD-based models, it was soon realized that describing the measured polarization of this state represents a much more difficult problem [18]. The polarization, corresponding to the orientation of the particle spin with respect to a chosen axis, can be accessed via a study of the polar (θ ) and azimuthal (φ ) production angles, relative to that axis, of the two-body decay products in the quarkonium rest frame. Their angular distribution W (θ , φ ) is parameterized as W (θ , φ ) ∝ 1 3 + λ θ 1 + λ θ cos 2 θ + λ φ sin 2 θ cos 2φ + λ θ φ sin 2θ cos φ , with the polarization parameters λ θ , λ φ and λ θ φ corresponding to various combinations of the elements of the spin density matrix of J/ψ production [19]. In particular, the two cases (λ θ = 1, λ φ = 0, λ θ φ = 0) and (λ θ = −1, λ φ = 0, λ θ φ = 0) correspond to the so-called transverse and longitudinal polarizations, respectively. At leading order, the high-p T production is dominated by gluon fragmentation and therefore the J/ψ would be expected to be transversely polarized [18]. However, the results from the CDF experiment at Tevatron showed that the J/ψ exhibits a very small polarization [20,21], an observation which was impossible to reconcile with the NRQCD prediction. As of today, on the experimental side, accurate results on inclusive and prompt (i.e., removing contributions from b-quark decays) J/ψ polarization have become available at LHC energies [22][23][24][25]. They confirm that this state shows little or no polarization in a wide rapidity (up to y = 4.5) and transverse momentum region (from 2 to 70 GeV/c), with the exception of the LHCb measurements at √ s = 7 TeV [24], where the value λ θ = −0.145 ± 0.027, corresponding to a weak longitudinal polarization, was obtained in the interval 2 < p T < 15 GeV/c and 2 < y < 4.5, in the helicity frame (its definition will be given later in Sec. 3). On the theory side, a huge effort was pursued in order to move to a complete next-to-leading order (NLO) description of the J/ψ production process [26,27], and to the calculation of the polarization variables [28,29]. Further important progress includes a quantitative evaluation of the contribution of feed-down processes (J/ψ coming from the decay of χ c and ψ(2S) states) on the polarization observables [30]. It was shown that at NLO there are rather large cancellations between contributions corresponding to the different possible combinations of the spin and angular momentum of the intermediate cc states, reaching a more satisfactory description of the absence of polarization observed in the data [31]. However, those descriptions usually require the inclusion of both cross section and polarization results in the fit of the LDME, leading to a more limited predictive power on the polarization observables and to large variations in the values of the extracted LDME values, depending on the set of data used for their determination. Finally, the description of the J/ψ production in the NRQCD framework was recently extended to the low-p T region, and the polarization parameters were studied in a color glass condensate (CGC) + NRQCD formalism, obtaining a fair agreement with LHC data at forward rapidity [32]. Measurements of the polarization parameters are also available for several bottomonium states, and in particular for the ϒ(1S), ϒ(2S) and ϒ(3S) resonances, which were shown to exhibit little or no polarization at LHC energies [33][34][35]. Approaches similar to that adopted for charmonium, which also need to take into account the rather complex feed-down decay structure for these states, lead to a fair agreement with the experimental results [36].
In this Letter, we move a step forward by presenting the first measurement of J/ψ and ϒ(1S) polarization in ultrarelativistic heavy-ion interactions performed by the ALICE Collaboration by studying Pb-Pb collisions at √ s NN = 5.02 TeV. Such collisions represent an important source of information for the investigation of the phase diagram of QCD [37], and in particular for the study of the properties of the quark-gluon plasma (QGP), a state of matter where quarks and gluons are not confined inside hadrons [38]. Among the experimental observables studied in heavy-ion collisions the suppression of heavy quarkonium production is a fundamental signal, since QGP formation prevents the binding of the heavy-quark pair due to the screening of the color charge [39] and, more generally, has strong effects on the spectral functions [40]. At LHC energies, another mechanism, corresponding to the (re)generation of charmonium states in the QGP and/or when the system hadronizes, becomes relevant [41,42], in particular at low p T , due to the large charm-quark multiplicity (> 100 pairs in a central Pb-Pb collision). The presence of a deconfined system may in principle affect also the polarization of quarkonium states. In Ref. [43] the observation of a partial transverse polarization for the J/ψ was predicted in case of QGP formation, due to a modification of the non-perturbative effects in the high energy-density phase. More generally, the observed prompt J/ψ are known to be a mixture of direct production and decay products from higher-mass charmonium states (ψ(2S), χ c ). In nuclear collisions, since suppression effects are expected to affect more strongly the less bound states, the relative contribution of direct and feed-down production would change with respect to that in pp collisions, and the overall measured polarization may be different according to the potentially different polarization of the various states [44,45]. On the other hand, the contribution of the regeneration mechanism in the J/ψ formation process by recombination of uncorrelated cc pairs is likely to give rise to unpolarized production at low p T . Finally, the possible presence of polarization is known to strongly affect the acceptance for J/ψ detection in the dilepton decay (up to 20-30% in ALICE [22]), and its measurement is an important requisite for an unbiased evaluation of the absolute yields in nuclear collisions. A first measurement of ϒ(1S) polarization in Pb-Pb collisions is also presented in this Letter, even if the corresponding candidate sample is smaller by a factor ∼30, leading to larger uncertainties. For such a state, considerations similar to those discussed for the J/ψ should hold, except that the contribution of the regeneration mechanism should be negligible due to the much lower multiplicity of bottom quarks with respect to charm.
The next sections of the Letter are organized as follows. Section 2 contains a short description of the experimental apparatus and some details on the data sample used in this analysis. The analysis procedure and the evaluation of systematic uncertainties are presented in Sec. 3, while the results on the J/ψ and ϒ(1S) polarization parameters λ θ , λ φ and λ θ φ are shown in Sec. 4. The conclusions are presented in Sec. 5.

Experimental setup and data sample
The measurement described in this Letter is performed with the ALICE detector [46,47], whose main components are a central barrel and a forward muon spectrometer. The latter covers the pseudorapidity region −4 < η < −2.5 and is used to detect muon pairs from quarkonium decays [48]. The muon spectrometer includes a hadron absorber made of concrete, carbon and steel with a thickness of 10 interaction lengths, followed by five tracking stations (cathode-pad chambers), with the central one embedded inside a dipole magnet with a 3 T·m field integral. Downstream of the tracking system, an iron wall filters out the remaining hadrons as well as low-momentum muons originating from pion and kaon decays, and is followed by two trigger stations (resistive plate chambers). Another forward detector, the V0 [49], composed of two scintillator arrays located at opposite sides of the interaction point (IP) and covering the pseudorapidity intervals −3.7 < η < −1.7 and 2.8 < η < 5.1, provides the minimum bias (MB) trigger which is given by a coincidence of signals from the two sides. Among the central barrel detectors, the two layers of the Silicon Pixel Detector (SPD), with |η| < 2 and |η| < 1.4 coverage, and corresponding to the inner part of the ALICE Inner Tracking System (ITS) [50], are used to determine the position of the interaction vertex. Finally, the Zero Degree Calorimeters (ZDC) [51], located on either side of the IP at ± 112.5 m along the beam axis, detect spectator nucleons emitted at zero degrees with respect to the LHC beam axis and are used to reject electromagnetic Pb-Pb interactions.
The analysis is based on events where, in addition to the MB condition, two opposite-sign tracks are detected in the triggering system of the muon spectrometer (dimuon trigger). The dimuon trigger selects tracks each having a transverse momentum above a threshold nominally set at p µ T = 1 GeV/c, corresponding to the value for which the single-muon trigger efficiency reaches 50% [52]. The single-muon trigger efficiency reaches a plateau value of 98% at ∼ 2.5 GeV/c.
The events are further characterized according to their centrality, i.e., the degree of geometric overlap of the colliding nuclei. It is estimated by means of a Glauber model fit to the V0 signal amplitude distribution [53,54], with more central events leading to a larger signal in the V0. In this analysis, events corresponding to the most central 90% of the inelastic Pb-Pb cross section are selected, as for these events the MB trigger is fully efficient and the residual contamination from electromagnetic processes is negligible.
The results of the analysis are obtained using the √ s NN = 5.02 TeV Pb-Pb data samples collected by the ALICE experiment during the years 2015 and 2018, corresponding to an integrated luminosity L int ∼ 750 µb −1 .

Data analysis
The J/ψ and ϒ(1S) candidates are formed by combining opposite-sign muons reconstructed using the tracking algorithm described in Ref. [48]. In order to reject tracks at the edge of the spectrometer acceptance, the condition −4 < η µ < −2.5 is required. In addition, tracks must have a radial transverse position at the end of the absorber in the range 17.6 < R abs < 88.9 cm. This selection is applied to remove tracks passing through the inner and denser part of the absorber, which are strongly affected by multiple scattering. For each muon candidate, a match between tracks reconstructed in the tracking system and track segments in the muon trigger system is required.
The J/ψ polarization parameters λ θ , λ φ and λ θ φ are studied as a function of transverse momentum in the intervals 2 < p T < 4, 4 < p T < 6 and 6 < p T < 10 GeV/c. For each p T interval, a two-dimensional (2D) grid of dimuon invariant-mass spectra is created, corresponding to intervals in cos θ and φ , where θ and φ are the polar and azimuthal emission angles, respectively, of the decay products in the J/ψ rest frame, with respect to the reference axis. More in detail, the 2D grid covers the fiducial region −0.8 < cos θ < 0.8 (17 intervals), 0.5 < φ < π − 0.5 rad (8 intervals, assuming a symmetric distribution around φ = π), with the choice of the boundaries as well as the width of the intervals dictated by acceptance considerations.
The analysis is performed choosing two different reference systems for the determination of the angular variables. In the Collins-Soper (CS) frame the z-axis is defined as the bisector of the angle between the direction of one beam and the opposite of the direction of the other one in the rest frame of the decaying particle, allowing therefore an evaluation of the polarization parameters with respect to the direction of motion of the colliding hadrons. In the helicity (HE) reference frame the z-axis is given by the direction of the decaying particle in the center-of-mass frame of the collision, and therefore the polarization can be evaluated with respect to the momentum direction of the J/ψ itself. The φ = 0 plane is the one containing the two beams in the J/ψ rest frame.
For each dimuon invariant-mass spectrum, the J/ψ raw yield is obtained by means of a binned maximum likelihood fit in the interval 2.1 < m µ µ < 4.9 GeV/c 2 . The background continuum is parameterized with a Gaussian distribution whose width varies linearly with the mass or, alternatively, with a fourth degree polynomial function times an exponential. The J/ψ signal is modeled with a pseudo-Gaussian function or with a Crystal Ball function with asymmetric tails on both sides of the peak [55].
The J/ψ mass is kept free in the fits, while for each interval (i, j) in (cos θ ,φ ) the width is fixed to , scaling the resonance width extracted from Monte Carlo (MC) simulations (σ i, j,MC J/ψ ) by the ratio between the width obtained by fitting the angle-integrated spectrum in data (σ J/ψ ) and MC (σ MC J/ψ ) for the p T interval under consideration. The parameters of the non-Gaussian tails of the resonance are kept fixed to the MC values. The ψ(2S) contribution, although comparatively negligible, is also taken into account in the fits, with its width and mass fixed in each fit to those of the J/ψ according to the relations σ ψ( with the Particle Data Group (PDG) masses taken from Ref. [1]. In Fig. 1 (left) an example of a fit to the invariant-mass spectrum in the J/ψ mass region is shown. Due to the stability of the extracted J/ψ parameters (mass, width), the fits were carried out directly on the sum of the 2015 and 2018 invariant mass spectra.  Figure 1: Examples of fits to the invariant-mass distributions in the helicity reference frame. The left plot corresponds to the J/ψ mass region, while on the right a fit to the ϒ(1S) mass region is shown. The fits are performed using an extended Crystal Ball function for the resonance signals, while the background is parameterized with a variable width Gaussian.
The J/ψ raw yields as a function of the angular variables are then corrected by the product of the acceptance and detector efficiency (A × ε), which is evaluated as a function of cos θ and φ on a 2D grid via MC simulations. The J/ψ are generated according to p T and y distributions directly tuned on data [56] via an iterative procedure [57], and their decay muons are propagated inside a realistic description of the ALICE setup, based on GEANT 3.21 [58]. The misalignment of the detection elements and the timedependent status of each electronic channel during the data taking period are taken into account as well. In the J/ψ generation an isotropic distribution of decay products, corresponding to the assumption of no polarization, is adopted. Due to the choice of relatively small (cos θ , φ ) intervals, the A × ε values for each interval are quite insensitive to the specific angular distribution assumed in the generation.
The three polarization parameters λ θ , λ φ and λ θ φ are obtained through χ 2 -minimization fits of the 2D J/ψ distributions, corrected for acceptance and efficiency, according to Eq. 1. For each combination of signal and background shape used in the fit to the dimuon invariant-mass spectra, a separate evaluation of the polarization parameters is carried out and their average is taken as the best estimate. The statistical uncertainty is given by the average of the statistical uncertainties of the 2D fits, while the root mean square of the results provides the systematic uncertainty on the signal extraction, with the absolute values ranging between 0.002 and 0.039. The overall procedure described above was checked beforehand with a MC closure test. The 2D fits on the (cos θ , φ ) distributions only allow a determination of the absolute value of λ θ φ , due to the presence of sin 2θ in the corresponding term that induces an ambiguity in its sign. It is checked that the values of λ θ and λ φ are stable against the choice of the sign of the λ θ φ term. In the following the λ θ φ values corresponding to the choice of a positive sign are quoted. Figure 2 illustrates an example of the fit to the angular distributions. For better visibility, both the distribution and the fitted function are projected along one dimension.  In addition to the systematic uncertainty related to the choice of the mass shapes for signal and background, several other sources are taken into account. First, an alternative procedure for extracting the J/ψ signal is carried out, by keeping its width as a free parameter in the invariant-mass fits. The corresponding results for the polarization parameters are then obtained and the averages of the values corresponding to fixing the width or not are taken as the central values for λ θ , λ φ and λ θ φ . Half the difference between the results obtained with free or MC-anchored widths is then considered as a further systematic uncertainty related to the signal extraction. This uncertainty is found to be the leading contribution to the total absolute systematic uncertainty on the polarization parameters, and ranges between 0.001 and 0.063, the latter value corresponding to the uncertainty on λ HE θ for 2 < p T < 4 GeV/c. Another source of systematic uncertainty is related to the evaluation of the trigger efficiency. The muon trigger response function as a function of the single muon transverse momentum p µ T can be obtained via MC or with a procedure based on data [59]. Small deviations are found for p µ T < 2 GeV/c which induce an effect on A × ε for the J/ψ. Therefore, the polarization parameters are re-calculated with A × ε values weighted in such a way to account for the deviations. The variation of the polarization parameters between the different trigger efficiency estimates is taken as the related systematic uncertainty, with values ranging from 0.001 to 0.043, the highest values being found for λ HE θ in 2 < p T < 4 GeV/c. The systematic uncertainty related to the evaluation of the muon tracking efficiency is found to be negligible for this analysis, allowing a significant reduction of the total systematic uncertainty with respect to previous pp analyses [23]. Indeed, although the difference between efficiencies calculated via MC or from data [59] is of the order of 2%, a detailed investigation has shown no dependence on the angular variables and therefore no effect on the polarization parameters.
Finally, the systematic uncertainty induced by the choice of the p T and y distributions used as an input for the calculation of A × ε is evaluated testing alternative p T and y parameterizations, which are obtained by varying within their uncertainties the default distributions directly tuned on Pb-Pb data. The polarization parameters extracted with the modified values of A × ε are compared with those obtained with the default input shapes and the corresponding systematic uncertainty extracted in this way is found to range between 0.001 and 0.030, with the largest value assigned to λ HE θ for 2 < p T < 4 GeV/c. The influence of the choice of the angular distributions of the J/ψ decay products for the A × ε calculation is also investigated by means of an iterative procedure on these input distributions. The effect is found to be negligible, also due to the fact that the 2D correction procedure on the angular variables is by definition relatively insensitive to the specific choice of the corresponding distributions. A summary of the values of all the absolute systematic uncertainties, which are considered as uncorrelated as a function of p T , is reported in Table 1. The total systematic uncertainties are obtained, for each parameter and p T interval, as the quadratic sum of the values.  A similar procedure is followed for the extraction of the ϒ(1S) polarization parameters. Due to the smaller candidate sample, integrated values over the kinematic interval 2.5 < y < 4, p T < 15 GeV/c are obtained. The main difference with respect to the 2D approach followed for the J/ψ is the use of a simultaneous fit to the 1D angular distributions [23], after integration over the other variables. The requirement p µ T > 2 GeV/c, which helps reducing the combinatorial background, is included [60]. The ϒ(1S) signal extraction in the various cos θ and φ intervals is performed by means of invariant-mass fits (see the right panel of Fig. 1 for an example). The functions chosen for the resonances are the same as in the J/ψ analysis (pseudo-Gaussian or Crystal Ball), the mass value is fixed to that obtained from a fit to the integrated invariant-mass distribution, while the width for each angular interval is fixed to the MC value scaled by the ratio of the widths between data and MC for the angle-integrated distributions. The tail parameters are fixed to MC values. The small contribution from ϒ(2S) is also included in the fits [60]. The background continuum is parameterized with a Gaussian distribution whose width varies linearly with the mass or, alternatively, with a second degree polynomial function times an exponential. The systematic uncertainty on the signal extraction is calculated with the same procedure adopted for the J/ψ. An uncertainty related to the choice of the signal width has also been considered, taken as the half-difference between the results obtained with the prescription described above and using as an alternative prescription the pure MC values. The uncertainty on the trigger efficiency is negligible, due to the additional requirement on the single-muon transverse momentum which selects a p T -region where the trigger efficiency is very high and its evaluation via data and MC is consistent. Finally, the procedure for the determination of the uncertainty related to the ϒ(1S) kinematic distributions used in the MC is the same as for the J/ψ. The total systematic uncertainties for the ϒ(1S) analysis are reported in Table 3, together with the results.

Results
The polarization parameters for J/ψ inclusive production in Pb-Pb collisions at √ s NN = 5.02 TeV in the helicity and Collins-Soper reference frames are shown in Fig. 3 and the corresponding numerical values are reported in Table 2. In Fig. 3, λ θ , λ φ and λ θ φ are also compared with the LHCb [24] and ALICE [23] measurements in pp collisions at √ s = 7 and 8 TeV, respectively.  For all the p T intervals and in both reference frames the values of the polarization parameters exhibit at most slight deviations from zero. In particular, λ HE θ indicates a slight transverse polarization at low p T (∼ 2.1σ effect, calculated using the Gaussian approximation), while λ CS θ shows a weak longitudinal polarization (∼ 2.1σ ). When increasing p T , the central values of λ θ become close to zero. All values of λ φ and λ θ φ are, in absolute value, smaller than 0.1, except for λ HE θ φ , which is −0.124 at low p T and deviates from zero by ∼ 2.4σ .
When comparing with the pp results, no significant difference is found with respect to ALICE results at √ s = 8 TeV, which are compatible with zero. A significant difference is found with respect to the higherprecision LHCb results at √ s = 7 TeV, reaching 3.3σ in the interval 2 < p T < 4 GeV/c in the helicity reference frame, where pp data [24] indicate a small but significant degree of longitudinal polarization, while the Pb-Pb results favor a slightly transverse polarization. In Pb-Pb collisions at LHC energies, a significant fraction of the detected J/ψ originates from the recombination of cc pairs in the QGP phase or when the system hadronizes. Moreover, the contribution from higher-mass charmonium states decaying (the LHCb markers were shifted horizontally by +0.3 GeV/c for better visibility) in the rapidity interval 3 < y < 3.5. The error bars represent the total uncertainties for the pp results, while for Pb-Pb statistical and systematic uncertainties are plotted separately as a vertical bar and a shaded box, respectively. In the left part of the plot the polarization parameters in the helicity reference frame are reported, in the right those for the Collins-Soper frame.
to J/ψ could vary between pp and Pb-Pb due to different suppression effects for each state in nuclear collisions. Therefore, the observed hint for a different polarization in pp and Pb-Pb might be a reflection of the different production and suppression mechanisms in the two systems, but more precise data, along with quantitative theory estimates, are needed for a definite conclusion. It should also be noted that the ALICE results refer to inclusive production, while LHCb has measured prompt J/ψ. However, as discussed in Ref. [22], the size of the non-prompt component is small in the covered p T region (of the order of 15% at high p T ) and its polarization was also measured to be small by CDF (λ HE θ ∼ −0.1 [21]), implying that the net effect of this source on inclusive J/ψ polarization should be negligible.
In Table 3 the values of the ϒ(1S) polarization parameters are shown. The λ θ values are consistent with zero, with large uncertainties that prevent a firm conclusion on the absence of polarization in nuclear collisions. The λ φ and λ θ φ values are also consistent with zero. The relatively smaller uncertainties for these parameters are related to a more uniform acceptance distribution as a function of the azimuthal angular variable. Table 3: ϒ(1S) polarization parameters in the helicity and Collins-Soper reference frames measured in Pb-Pb collisions at √ s NN = 5.02 TeV in the rapidity interval 2.5 < y < 4 and for transverse momentum p T < 15 GeV/c. The first uncertainty is statistical and the second systematic.

Conclusions
The first measurement of the polarization parameters for J/ψ production in nuclear collisions at LHC energies was carried out by the ALICE Collaboration in Pb-Pb interactions at √ s NN = 5.02 TeV. The λ θ , λ φ and λ θ φ parameters were evaluated in the helicity and Collins-Soper reference frames in the rapidity interval 2.5 < y < 4 and in the transverse momentum interval 2 < p T < 10 GeV/c. All the parameter values are close to zero, with a ∼ 2.1σ indication for a small transverse polarization in the helicity frame at low p T , and a corresponding indication for a small longitudinal polarization in the Collins-Soper frame (∼ 2.1σ effect). When comparing these results with pp data taken at higher energy at the LHC, an interesting feature is a significant difference in λ HE θ with respect to the LHCb results which showed instead a small longitudinal polarization in a similar kinematic domain. This first result obtained for J/ψ in nuclear collisions and described in this Letter represents therefore a starting point for future studies connecting such features with the known differences in the production mechanisms between pp and nucleus-nucleus collisions. Results were also obtained for the first time for the ϒ(1S) polarization, integrated over p T and y, showing, within the large uncertainties of the measurement, values compatible with the absence of polarization.