Measurement of the azimuthal anisotropy of $\Upsilon$(1S) and $\Upsilon$(2S) mesons in PbPb collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

The second-order Fourier coefficients ($v_2$) characterizing the azimuthal distribution of $\Upsilon$(1S) and $\Upsilon$(2S) mesons arising from PbPb collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV are studied. The $\Upsilon$ mesons are reconstructed in their dimuon decay channel, as measured by the CMS detector. The data set corresponds to an integrated luminosity of 1.7 nb$^{-1}$. The scalar product method is used to extract the $v_2$ coefficients of the azimuthal distribution. Results are reported for the rapidity range $|y|$ $\lt$ 2.4, with the transverse momentum 0 $\lt$ $p_\mathrm{T}$ $\lt$ 50 GeV/$c$, and in three centrality ranges of 10-30%, 30-50% and 50-90%. In contrast to the J/$\psi$ mesons, the measured $v_2$ values for the $\Upsilon$ mesons are found to be consistent with zero.


Introduction
The primary motivation for studying high-energy heavy ion collisions is to better understand the hot and dense matter produced in these interactions [1]. Lattice quantum chromodynamics (LQCD) calculations indicate that at sufficiently high temperatures an environment is created for a crossover to occur from hadronic matter to a strongly interacting system of deconfined quarks and gluons [2]. This deconfined medium is called the "quark-gluon plasma" (QGP) [3]. One of the most prominent signatures of QGP formation is that the production of quarkonia, the bound states of a heavy quark and its antiquark, is suppressed with respect to expectations from scaling the yields in proton-proton collisions by the number of binary nucleon-nucleon (NN) collisions. This suppression arises because the quarkonia binding is weakened by color screening caused by the surrounding partons in the medium [4], and therefore the extent of the quarkonia suppression is expected to be sequentially ordered by the binding energies of the quarkonia states. In heavy ion collisions, the majority of heavy quarks are produced from hard scattering processes at an early stage and therefore can interact with the QGP throughout its entire evolution.
Because of the binding energy dependence of the screening, the bottomonium states (Υ(1S), Υ(2S), Υ(3S), χ b , etc.) are particularly useful probes to understand the space-time evolution of the QGP. The sequential suppression of the yield of Υ(nS) states was first observed by CMS [5, 6] at a NN center-of-mass energy √ s NN = 2.76 TeV. More recently, results with improved statistical precision have been reported by both the ALICE Collaboration [7] and CMS Collaboration [8,9] at √ s NN = 5.02 TeV. The suppression of the Υ(1S) meson has also been studied at √ s NN = 200 GeV at RHIC [10]. The available experimental data, spanning from 0.20 to 5.02 TeV, have provided new insight into the thermal properties of the QGP [11].
The screening due to the QGP can also result in an azimuthal asymmetry in the observed yields of quarkonia. In non-central heavy ion collisions, the produced QGP has a lenticular shape in the transverse plane. Consequently, the average path length for quarkonia traveling through the medium depends on the direction taken with respect to this shape, with a larger suppression in the direction of the longer axis [12][13][14][15]. The anisotropic distribution of particles can be characterized by the magnitudes of the Fourier coefficients (v n ) of the azimuthal correlation of particles, through the relation dN dφ Here φ is the angle of the particles and Ψ n is the angle characterizing the symmetry plane of the nth harmonic based on the charged particle density distribution [16]. By studying the azimuthal distribution of the quarkonia, it is possible to develop a more comprehensive understanding of the dynamics of their production. TeV, respectively. Significant evidence for an azimuthal asymmetry in the particle yields is found for both experiments. However, the J/ψ asymmetry can also arise from the late-stage recombination of uncorrelated charm quark pairs inside the QGP, which acts in the opposite direction as color screening. The relative contribution of recombination to the inclusive cross section of bottomonium states is expected to be smaller than for charmonium states [20][21][22][23]. For this reason, the azimuthal anisotropy of bottomonia better reflects the screening that occurs during its passage through the QGP medium. The ALICE Collaboration has measured the v 2 for Υ(1S) mesons at forward rapidity (2.5 < y < 4) using 5.02 TeV PbPb collisions [24]. track from the nearest primary vertex must be less than 0.3 and 20 cm in the transverse and the longitudinal directions, respectively. To assure high efficiency, the individual muons were required to have p µ T > 3.5 GeV/c, with |η µ | < 2.4. Pairs of oppositely charged muons are fitted with a common vertex constraint and used in this analysis if the confidence level of the fit is greater than 1%.

Analysis
The Υ meson signals and combinatorial backgrounds for selected dimuon pairs are separated using fits to the invariant mass (m inv ) distribution. The range of the studied mass spectrum is 8-14 GeV/c 2 , which covers both Υ(1S) and Υ(2S) resonance peaks, as well as the Υ(3S) mesons.
Acceptance and efficiency correction factors are used as event-by-event weights when the invariant mass distribution and the average v 2 values are computed. Acceptance is defined as the probability for both the decay muons originated from an Υ meson to be within p µ T > 3.5 GeV/c, |η µ | < 2.4. It is determined as a function of p T using PYTHIA 8.212 [30] Monte Carlo (MC) generator with tune CP5 [31], where Υ mesons from 3 S 1 , 3 P J , 3 D J bottomonium states form via the color-singlet [32, 33] and color-octet [34] mechanisms. The MC events are then reweighted to match the p T spectra to the PbPb data. The efficiencies for reconstruction, muon identification, and trigger selection are calculated by a full simulation of the CMS detector using GEANT4 [35]. Each Υ event is embedded in a HYDJET 1.9 [36] PbPb event to reproduce the hydrodynamic background environment. The average acceptance and efficiency factors are 53.7% (56.7%) and 56.9% (60.0%) for Υ(1S) (Υ(2S)), respectively. A systematic uncertainty is assigned based on the difference found between the MC results and experimental data. The MC and experimental data are compared using the tag-and-probe method with single muons obtained from J/ψ meson decays [37,38]. The efficiency ratio of data over simulation is applied as a weight for the derivation of the final dimuon efficiency.
The v 2 values of Υ candidates are determined using the scalar product (SP) method, in a similar way as the one employed to determine the elliptic flow of prompt D 0 mesons in PbPb collisions [39]. Using this method, the particle asymmetries associated with different subevents are characterized by Q-vectors that are determined for different η ranges using either tracker or HF calorimeter data [25]. The Q-vectors are given by Q 2 = ∑ M k=1 ω k e 2iφ k , where the sum is over the multiplicity of particles for the tracker or the number of towers for the HF calorimeters, φ is the azimuthal angle of the particle or the tower, and ω is the weight given by the p T of a particle in the tracker or the transverse energy deposited in a tower. A three subevent method [25,39] is used, employing the tracker, with |η| < 0.75, and the two HF calorimeters that cover the ranges 3 < η < 5 (HF+) and −5 < η < −3 (HF-). The Q-vector of an Υ candidate is defined as Q 2,Υ = e i2φ where φ is the azimuthal angle of the candidate. The v 2 coefficient is then given by Here, Q 2A and Q 2B correspond to subevents based on the HF detectors, and Q 2C is based on the tracker. The subscripts A and B refer to either HF+ or HF-, depending on the rapidity of the Υ candidate. In order to avoid autocorrelations and to reduce nonflow effects, the η gap between Υ candidates and the subevent detector is required to be at least 3 units. For this reason, HF+ is selected for A (B) when the Υ candidate is produced at negative (positive) rapidity.
The brackets denote averages over all Υ candidate events, taking the real part of the Q-vector products. The denominator in Eq. (2) is a resolution correction that depends on the particle multiplicity in the HF detector referenced by the A subscript [25]. The lower panel of Fig. 1 shows an example of the resulting v 2 distribution as a function of the Υ candidate invariant mass. The separation of signal and background is done in two steps. First, the yields of Υ signals are extracted from the invariant mass distribution without using the v 2 information. The purpose of this step is to find the best parameters for the probability distribution describing the Υ mass peaks. Unbinned log-likelihood fits to the mass distribution are used to obtain the best fit to the data. The analysis is done using a single rapidity range with |y| < 2.4 and, thus, uses Υ candidates obtained in both the barrel and endcap regions of CMS. For this reason, the signal mass distribution (Sig) for each Υ meson is formed by the sum of two Crystal Ball (CB) functions [40], to account for the different mass resolutions of the two regions. For both CB functions, the mass and the radiative tail parameters are constrained to be equal since these are not sensitive to the detector resolution. The Υ(1S) mass is taken as a free parameter of the fit, allowing for a possible scaling error in the momentum calibration for the reconstructed tracks. The mean (m 0 ) and width (σ) parameters for the excited states (Υ(2S) and Υ(3S)) are found by scaling the Υ(1S) values by the ratio of published masses [41]. Other parameters for the excited states are constrained to be identical to those of the Υ(1S) mesons. The normalization parameters N sig,1S , N sig,2S and N sig,3S are free parameters of the fit. More details about the parameters for the invariant mass fit can be found in Ref. [9].
The background mass spectrum (Bkg) is modeled by a function obtained by multiplying a realvalued error function and an exponential function. The error function is used to reproduce a kinematic enhancement found at low mass, resulting from the single muon p T threshold. The exponential component is motivated by the exponentially falling structure of the combinatorial background in the high mass region. There are four parameters characterizing the shape of background, the µ and σ parameters of the error function, the decay length of the exponential function, and an overall normalization parameter.
In the second step, the invariant mass and v 2 distributions are fitted simultaneously using binned χ 2 fits. The signal function parameters are fixed to the values obtained from the previous step, although the normalizations N sig,1S , N sig,2S , and N sig,3S and the background function parameters are left free. The v 2 dependence on m inv is taken as is the v 2 value for the Υ(iS) mesons and is assumed to be independent of m inv . The α i (m inv ) coefficients are the fractions that the Υ(iS) states occupy as a function of invariant mass, as determined from the fit. They are defined as The background v 2 value, v Bkg 2 (m inv ), is modeled as a second-order polynomial function of the invariant mass. Figure 1 shows an example of a simultaneous fit to the mass (upper panel) and v 2 distributions (lower panel).

Systematic Uncertainties
The systematic uncertainties in this analysis come from various sources, including the modeling of signal and background probability functions, as well as acceptance and efficiency corrections (muon identification and trigger). The uncertainties for the various sources are studied as functions of p T and centrality. However, within the limits of the evaluation, no clear p T or centrality trend is found. For the uncertainty associated with the signal, an alternative function is formed by adding one CB and one Gaussian function having equal mass parameters. The differences in the v 2 values from the nominal one are taken as the uncertainties. The difference is typically less than 0.0025 for the Υ(1S) mesons and is equal to 0.014 for the Υ(2S) mesons. The uncertainties resulting from holding the final fit parameters constant are studied by releasing each parameter, one by one, and redoing the fit with the additional free parameter. The largest deviation of the resulting v 2 values from the nominal result is taken as the uncertainty. The difference is found typically in the range 0.0006 to 0.005 for the Υ(1S) mesons and is 0.0093 for the Υ(2S) mesons. To determine the uncertainty related to the background modeling, the functions for both the invariant mass distribution and v The uncertainty resulting from the acceptance and efficiency corrections is evaluated by comparing the results with and without weighting the simulated p T spectra. By virtue of the cylin-drical symmetry of the detector, the effect of these corrections is minor. An additional systematic uncertainty for the efficiency is assigned based on the tag-and-probe method. The difference of the v 2 values with and without applying the tag-and-probe correction is assigned as the uncertainty and is found to be in the range 0.002-0.012. Finally, uncertainties resulting from the hadronic event selection are considered. The collision event filter is varied, thus allowing for a migration of Υ(1S) and Υ(2S) mesons across centrality boundaries. The effect on the measured v 2 value is taken as a systematic uncertainty. For Υ(1S) the differences are typically below 0.0008 and for Υ(2S) the uncertainty is 0.0067. Individual systematic uncertainties are expected to be uncorrelated and are therefore added in quadrature for the final assigned values. The absolute systematic uncertainty in the v 2 values for the Υ(1S) meson range in most intervals from 0.002 to 0.015, which is smaller than the statistical uncertainty. The largest uncertainty corresponds to the most peripheral and lowest p T interval measured, which is dominated by the uncertainty related to the v Bkg 2 (m inv ) modeling and hadronic event selection. The uncertainty for the Υ(2S) v 2 value, which is obtained for the integrated interval only, is 0.037.

Results
The p T integrated results are shown in Fig. 2 [42,43], Hong and Lee [44,45], Bhaduri et al. [46], and Reygers et al. [47]. All results are for the rapidity range |y| < 2.4. The vertical bars denote statistical uncertainties, and the rectangular boxes show the total systematic uncertainties.
In Fig. 2 (right panel), the p T dependence of Υ(1S) meson v 2 values is measured for the 10-90% centrality interval. The v 2 values are consistent with zero in the measured p T range, except for the 6 < p T < 10 GeV/c interval that shows a 2.6 σ deviation from zero. The result is also compared with theoretical predictions from five different approaches. The green shaded area (Du,

Rapp [22]
) is a calculation using the kinetic rate equation framework for the time evolution of bottom quarks. In this model, the in-medium effect is applied with a lattice QCD based equation of state. Note that it is calculated for the centrality interval 20-40%. The red line (Yao et al. [42,43]) is computed from a real-time simulation of heavy quarks using coupled Boltzmann transport equations. The dissociation and recombination rates are calculated from potential non-relativistic QCD for this model. The dashed violet curve (Hong, Lee [44,45]) is another kinetic model prediction in which the heavy quark collisional terms are obtained from Bethe-Salpeter amplitude and hard thermal loop resummation. A diffusion constant D(2πT) = 6, where T is the temperature of the medium, was used. With this parameter, the nuclear modification factor, which is the ratio of the quarkonium yield in a nucleus-nucleus collision to the corresponding yield in a proton-proton (pp) collision scaled by the number of binary collisions, is consistent with the CMS result in Ref. [9]. The dashed brown line (Bhaduri et al. [46]) shows the result of 3+1d quasiparticle anisotropic hydrodynamics framework. The initial conditions for temperature and shear viscosity to entropy density ratio are tuned to the identified hadron spectra and flow harmonics in 5.02 TeV PbPb data [48]. The gray shaded area (Reygers et al. [47]) is another hydrodynamics calculation based on the blast-wave fits using the data of lighter particles, which was scaled down by a factor of 1/10. The v 2 values are computed in three centrality intervals (10-20%, 40-50%, 60-70%) and the gray envelope is constructed to minimally cover all results. Despite their slightly different p T dependence, all models in the figure derive small v 2 values for p T < 10 GeV/c, which is consistent with the experimental result. For the highest p T range, however, the kinetic model (Hong, Lee) and the blast-wave model begin to deviate from the observed behavior. In particular, the blast-wave model predicts a v 2 value much larger than found. In this model, the elliptic flow of Υ(1S) is expected to sharply rise from p T = 9 GeV/c, following the mass ordering of v 2 (p T ) by lighter particles used for the fit [47]. The p T differential results for Υ(1S) mesons in each centrality interval are shown in Fig. 3 and the v 2 values are found to be consistent with zero within uncertainties. A similar result has been obtained by the ALICE Collaboration [24], where the measurement was done in a complementary rapidity region (2.5 < y < 4), as shown in Fig. 4. The CMS data points in this figure were analyzed using the same p T range as the ALICE result to allow for a comparison where only the rapidity ranges of the two measurements differ. The ALICE Collaboration also measured inclusive J/ψ v 2 values (also shown in Fig. 4) with the same kinematic conditions as their Υ(1S) results and there find a significant and finite v 2 value. Together, the CMS and ALICE results indicate that the geometry of the medium has little influence on the Υ(1S) yields [22,49] and that recombination is not a dominant process in the production of this meson. The results also indicate that the path-length dependence of Υ(1S) suppression is small.
The current results support the assumption of a high melting temperature of the Υ(1S) meson. Consequently, the dissociation can only take place in the earliest stages of the collision. This makes the Υ(1S) meson less sensitive to the path length than other quarkonia states.

Summary
The v 2 coefficients for Υ(1S) and Υ(2S) mesons are measured in PbPb collisions at a nucleonnucleon center-of-mass energy of 5.02 TeV. Results are reported for the rapidity range |y| < 2.4, in the transverse momentum interval 0 < p T < 50 GeV/c, and in three centrality classes of 10-30%, 30-50%, and 50-90% for the Υ(1S) meson, while the centrality interval 10-90% is used for the Υ(2S) meson. The v 2 values are observed to be compatible with zero for the Υ(1S) meson in the measured kinematic and centrality intervals. The observation contrasts with the positive v 2 values reported for J/ψ mesons, suggesting different medium effects for charmonia and bottomonia. This confirms, with higher statistical precision and with an extended p T range, the result from the ALICE experiment measured in a complementary rapidity region. The measured values of v 2 for p T < 10 GeV/c are consistent with the theoretical predictions based on the kinetic model, the hydrodynamics framework, and the blast-wave model. For the high p T region, the Υ(1S) result provides an important constraint for models. The v 2 value found for Υ(2S) mesons, which is being reported for the first time, is also consistent with zero. As there are expected to be differences in the various processes through which the QGP affects Υ(1S) and Υ(2S) mesons, these measurements provide new inputs for the study of bottomonia production in heavy ion collisions.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. [  [26] I. Selyuzhenkov and S. Voloshin, "Effects of nonuniform acceptance in anisotropic flow measurements", Phys. Rev. C 77 (2008) 034904, doi:10.1103/PhysRevC.77.034904, arXiv:0707.4672.