Strong collective behavior is found in the azimuthal correlations of particles emitted in relativistic nucleus-nucleus (AA) collisions at the BNL RHIC [1, 2, 3, 1, 4] and at the CERN LHC [5, 6, 7, 8, 9, 10]. These correlations, which are long-range in pseudorapidity ($\eta$), suggest the formation of a strongly interacting quark-gluon plasma (QGP) that exhibits nearly ideal hydrodynamic behavior [11, 12, 13]. The azimuthal correlation structure of emitted particles is typically characterized by its Fourier components [14]. In particular, within a hydrodynamic picture, the second and third Fourier anisotropy components are known as elliptic ($v_2$) and triangular ($v_3$) flow, respectively, and reflect the QGP medium response to the initial collision geometry and its fluctuations [15, 16, 17]. In recent years, similar long-range collective azimuthal correlations have also been observed in events with high final-state particle multiplicity in proton-proton ($p$$p$) [18, 19, 20, 21], proton-nucleus ($p\text{A}$) [22, 23, 24, 25, 26, 27, 28, 29, 30], and lighter AA collisions [31, 32, 33], raising the question of whether a fluid-like QGP is created in these much smaller systems. While experimental measurements in these small systems are consistent with the hydrodynamic expansion of a tiny QGP droplet, alternative scenarios based on gluon saturation in the initial state also claim to capture the main features of the correlation data (recent reviews are provided in Refs. [34, 35]).

Heavy-flavor quarks (charm and bottom) are primarily produced at a very early stage via initial hard scattering because of their large masses. As such, they are largely decoupled from the bulk production of soft gluons and light-flavor quarks at a later stage in AA collisions, and thereby probe the properties and dynamics of the QGP through its entire evolution [36]. A strong elliptic flow ($v_2$) signal has been observed for open heavy-flavor $D^0$ mesons in both AuAu collisions at RHIC [37] and PbPb collisions at the LHC [38, 39, 40], suggesting that charm quarks may develop strong collective flow behavior. Furthermore, a recent measurement of the elliptic flow of $J/\psi$ mesons in PbPb collisions at $\sqrt{s_{{}_{NN}}}=5.02\text{TeV}$ [41] has provided additional evidence for the collective behavior of charm quarks in the QGP.

In the study of collectivity in small systems, such as that occurring in $p$$p$ or $p\text{Pb}$ collisions, a key open question is whether the strong collective behavior observed for bulk constituents in high-multiplicity events also extends to charm and bottom quarks. Long-range correlations involving inclusive muons at high transverse momentum ($p_T$) reveal a hint of heavy-flavor quark collectivity in $p\text{Pb}$ collisions [42]. Furthermore, the recent observation of a significant elliptic flow signal for prompt $D^0$ mesons in $p\text{Pb}$ collisions has provided evidence for charm quark collectivity in a small system [43]. The $v_2$ signal for $D^0$ mesons is found to be smaller than that of light-flavor hadrons at a given $p_T$, indicating that in these small systems there is a weaker collective motion for charm quarks, as compared to that of the bulk medium, than found in large AA systems. However, as the $D^0$ meson carries both a light and a charm quark, the relative contribution of these different flavor quarks to the observed $v_2$ signal is not fully constrained. Without detailed theoretical modeling, a scenario is not excluded where the $D^0$ meson $v_2$ signal is entirely carried by the light-flavor quark. The observation of an elliptic flow signal for $J/\psi$ mesons in a small system could provide more direct evidence of charm quark collectivity and could impose new constraints on the collective dynamics of heavy-quark production in such collisions. Furthermore, heavy-quark collectivity may also provide a hint of how, in small systems, hard probes interact with the QGP [36], assuming this is formed. A measurement of inclusive $J/\psi$ (combined charmonia and $J/\psi$ mesons from decay of open beauty hadrons) $v_2$ in $p\text{Pb}$ collisions was reported in Ref. [44], where positive $v_2$ coefficients were found in the range of $3<p_T<6\text{GeV}$ with center-of-mass rapidities $-4.46<y_{cm}<-2.96$ or $2.03<y_{cm}<3.53$. A recent model calculation of $J/\psi$ $v_2$ in $p\text{Pb}$ collisions suggests little $v_2$ signal arising from final-state interactions between charm quarks and the QGP medium [45].

This Letter presents the first measurement of prompt $J/\psi$ meson elliptic flow (excluding contributions from $B$ hadron decays) from long-range two-particle correlations in very high multiplicity $p\text{Pb}$ collisions at $\sqrt{s_{{}_{NN}}}=8.16\text{TeV}$. The $v_2$ harmonics for prompt $J/\psi$ mesons in the ranges $-2.86<y_{cm}<-1.86$ and $0.94<y_{cm}<1.94$ are determined over a wide $p_T$ range from 0.2 to 10$\text{GeV}$. To estimate the possible residual contribution from back-to-back jet-like correlations, the $v_2$ values are also presented after subtracting correlations obtained from low-multiplicity $p\text{Pb}$ events (denoted as $v_2^{\text{sub}}$), where jet-like correlations are assumed to dominate. The results are compared to those of the light strange-flavor $K_S^0$ and $\Lambda$ hadrons, and the open heavy-flavor prompt $D^0$ meson, which were previously reported by CMS [43] in the same $p_T$ range but in a different rapidity range of $-1.46<y_{\text{cm}}<0.54$. In order to explore possible collectivity at the partonic level, a comparison is also presented in terms of the transverse kinetic energy per constituent quark ($K{E}_T$/$n_q$, where $K{E}_T=\sqrt{m^2+p_T^2}-m$, and $n_q$ is the number of constituent quarks).

The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume, there are four primary subdetectors including a silicon pixel and strip tracker detector, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. Iron and quartz-fiber Cherenkov hadron forward (HF) calorimeters cover the range $3.0<\left|\eta \right|<5.2$. Muons are measured in the range $\left|\eta \right|<2.4$ in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers. The silicon tracker measures charged particles within the range $\left|\eta \right|<2.5$. For charged particles with $1<p_T<10\text{GeV}$ and $\left|\eta \right|<1.4$, the track resolutions are typically 1.5% in $p_T$ and 25–90 (45–150)$\mu \text{m}$ in the transverse (longitudinal) impact parameter [46]. A detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [47].

The $p\text{Pb}$ data at $\sqrt{s_{{}_{NN}}}=8.16\text{TeV}$ used in this analysis were collected in 2016, and correspond to an integrated luminosity of 186${\text{nb}}^{-\text{1}}$. The beam energies are 6.5$\text{TeV}$ for the protons and 2.56$\text{TeV}$ per nucleon for the lead nuclei. Because of the asymmetric beam conditions, particles selected in the laboratory rapidity range of $1.4<y_{\text{lab}}<2.4$ ($-2.4<y_{\text{lab}}<-1.4$) have a corresponding nucleon-nucleon center-of-mass frame rapidity range of $0.94<y_{cm}<1.94$ ($-2.86<y_{cm}<-1.86$), with positive rapidity defined in the proton beam direction. To minimize statistical uncertainties, the quoted $J/\psi$ meson $v_2$ results combine the individual values obtained for the proton and lead beam directions.

The $p\text{Pb}$ data are analyzed in different ranges of $N_{trk}^{offline}$, where $N_{trk}^{offline}$ is the number of primary charged particle tracks [46] with $\left|\eta \right|<2.4$ and $p_T>0.4\text{GeV}$. The main results are obtained with events in the high-multiplicity range $185\le N_{trk}^{offline}<250$. To select these events, dedicated triggers were developed, as discussed in Refs. [48, 49]. Events with $N_{trk}^{offline}<35$ are also used to estimate the possible contribution of residual back-to-back jet-like correlations. These lower-multiplicity events are selected online with a hardware-based trigger requiring two muon candidates in the muon detectors with no explicit momentum or rapidity threshold [50]. In the offline analysis, hadronic collisions are selected by requiring at least one HF calorimeter tower with more than 3$\text{GeV}$ of total energy in each of the two HF detectors. Events must contain a primary vertex close to the nominal interaction point of the beams, within 15 cm along the beam direction, and 0.2 cm in the plane transverse to beam direction. The $N_{trk}^{offline}$ range limits correspond to fractional inelastic cross sections from 100 to 57% for $N_{trk}^{offline}<35$, and from 0.33 to 0.01% for $185\le N_{trk}^{offline}<250$, respectively.

The offline muon reconstruction algorithm starts either by finding tracks in the muon detectors, which are then fitted together with tracks reconstructed in the silicon tracker (global muons), or by extrapolating tracks from the silicon tracker to match a hit on at least one segment of the muon detectors (tracker muons). The muon candidates are required to pass the identification criteria of the particle-flow algorithm [51], which suppresses contamination of “punch-through” hadrons misidentified as muons, based on energy deposition in the calorimeters. The soft muon selection criteria are also imposed, as defined in Ref. [52], to further improve the purity of muons.

The $J/\psi$ meson candidates are formed from pairs of oppositely charged muons, originating from a common vertex. Based on the vertex probability distributions for signal and background candidates, the probability that the dimuon pair shares a common vertex is required to be larger than 1%, lowering the background from random combinations as well as from semileptonic decays of bottom and charm hadrons. Because of the long lifetime of $b$ hadrons compared to that of $J/\psi$ mesons, the nonprompt $J/\psi$ meson component can be reduced by placing constraints on the pseudo-proper decay length [53]. This is defined by ${\ell }_{J/\psi }^{3D}=L_{xyz}{m}_{J/\psi }/\left|p_{\mu \mu }\right|$, where $L_{xyz}$ is the distance between the primary and dimuon vertices, $m_{J/\psi }$ is the Particle Data Group [54] world average value of the $J/\psi$ meson mass (assumed for all dimuon candidates), and $p_{\mu \mu }$ is the dimuon momentum. The upper limit (decreasing as a function of $p_T$) imposed on the ${\ell }_{J/\psi }^{3D}$ value is based on Monte Carlo (MC) studies with simulated event samples of pythia 8.209 [55, 56], and found to reject 75–90% (from low to high $p_T$) of nonprompt $J/\psi$ mesons, largely independent of multiplicity. The residual nonprompt $J/\psi$ meson fraction in the data is estimated to be approximately 5% across the full $p_T$ range, and its effect on the $v_2$ measurement is propagated as a systematic uncertainty, as described in Section 5.

The azimuthal anisotropy of $J/\psi$ mesons is extracted from the long-range ($\left|\Delta \eta \right|>1$) two-particle azimuthal correlations, following an identical procedure to that described in Refs. [21, 27, 43]. A two-dimensional (2D) correlation function is constructed by pairing each $J/\psi$ candidate with reference primary charged-particle tracks with $0.3<p_T<3$$\text{GeV}$ and $\left|\eta \right|<2.4$ (denoted “ref” particles), and calculating

To extract the genuine $v_2$ values of the $J/\psi$ meson signal ($v_2^S$), the contribution from background candidates ($v_2^B$) has to be subtracted from the $v_2$ values of all $J/\psi$ meson candidates, as obtained in the previous step. The procedure is to first fit the dimuon mass spectrum with a function composed of three components: two Crystal Ball functions [58] with different widths but common mean and tail parameters for the $J/\psi$ signal (the tail parameters are fixed to the values obtained from simulation), $S\left(m_{\text{inv}}\right)$, and an exponential function to model the combinatorial background, $B\left(m_{\text{inv}}\right)$. Then, the signal plus background $v_2^{S+B}\left(m_{\text{inv}}\right)$ distribution is fitted with: {linenomath}

where {linenomath}

Here, $v_2^B\left(m_{\text{inv}}\right)$ for the background $J/\psi$ candidates is modeled as an exponential function of the invariant mass, and $\alpha \left(m_{\text{inv}}\right)$ is the $J/\psi$ signal fraction obtained from the mass spectrum fit. An example of fits to the mass spectrum and $v_2^{S+B}\left(m_{\text{inv}}\right)$ in the $p_T$ interval 6.0–8.0$\text{GeV}$ for the multiplicity range $185\le N_{trk}^{offline}<250$ is shown in Fig. 1. The residual contribution of back-to-back dijets to the measured $v_2$ results is estimated from low-multiplicity $p\text{Pb}$ events and is removed from the signal after accounting for the jet yield ratio of the selected events, following a jet subtraction procedure similar to that established in Refs. [57, 21, 43]. The Fourier coefficients, $V_{n\Delta }$, extracted from Eq. (2) for $N_{trk}^{offline}<35$, are subtracted from the $V_{n\Delta }$ coefficients obtained in the high-multiplicity region, with {linenomath}

Here, $Y_{\text{jet}}$ represents the jet yield obtained by integrating the difference of the short-range ($\left|\Delta \eta \right|<1$) and long-range event-normalized associated yields for each multiplicity class. The ratio, $Y_{\text{jet}}$/$Y_{\text{jet}}(N_{trk}^{offline}<35)$, is introduced to account for the enhanced jet correlations resulting from the selection of higher-multiplicity events. For $p_T\left(J/\psi \right)<4.5\text{GeV}$, the jet yield ratio cannot be directly estimated from the two-particle azimuthal correlations, as the $J/\psi$ candidates tend to have larger $\eta$ values than the acceptance for charged particles. Therefore, the value is assumed to be the same as that for the high-$p_T$  region, where no $p_T$ dependence has been observed. It was also previously observed that the values of jet yield ratio for $D^0$ and strange particle species show little dependence on $p_T$ over the full $p_T$ range [43].

Sources of systematic uncertainties on the prompt $J/\psi$ meson $v_2$ measurement include the $J/\psi$ meson yield correction (acceptance and efficiency correction derived from pythia simulation), the nonprompt $J/\psi$ meson contamination, the background $v_2^B\left(m_{\text{inv}}\right)$ function form, the signal and background invariant mass PDF, the jet subtraction procedure, the contamination of events containing more than one $p\text{Pb}$ interaction (pileup), and the trigger bias. In this Letter, the quoted uncertainties in $v_2$ are absolute values, and are found to have no dependence on $p_T$, except those for the jet subtraction procedure. Systematic uncertainties originating from different sources are added in quadrature to obtain the overall systematic uncertainty shown as boxes in the figures.

To evaluate the uncertainties arising from the efficiency correction to the $J/\psi$ meson yield, the $v_2$ values are compared to the uncorrected ones, yielding an uncertainty of 0.008. The effect on the measured $v_2$ due to the residual contribution from nonprompt $J/\psi$ mesons is evaluated by varying the ${\ell }_{J/\psi }^{3D}$ requirement such that the nonprompt $J/\psi$ meson yield is doubled. The $v_2$ values are found not to change by more than $\pm 0.004$, which is assigned as the systematic uncertainty due to the $J/\psi$ meson yield correction. Possible differences in the rejection efficiency of nonprompt $J/\psi$ mesons between data and simulation are investigated and found to be negligible. The systematic uncertainties from the background $v_2$ function form are evaluated by comparing $v_2^B\left(m_{\text{inv}}\right)$ values based on first-, second-, and third-order polynomial fits to the background distribution. The resulting $J/\psi$ signal $v_2$ values are found to vary by less than 0.009. Systematic effects related to signal invariant mass PDF are found to be negligible by releasing, one at a time, the fixed tail parameters of the Crystal Ball functions. The variation of $v_2$, while changing the background invariant mass PDF to a second- or third-order polynomial function is also found to be negligible. In the jet subtraction procedure, the statistical precision of the jet yield ratio is limited. The $v_2^{\text{sub}}$ results are found to be consistent within $\pm 0.002$ to $\pm 0.014$ (increasing with $p_T$) when varying the jet yield ratio by its statistical uncertainty. The systematic uncertainties from the potential pileup effect and the trigger bias are taken to be the same as for inclusive charged particles in Ref. [49], where they can be established with good statistical precision. The pileup and trigger bias uncertainties are negligible compared to the other sources of systematic uncertainties, as the fraction of residual pileup events is only a few % and the trigger efficiency is close to 100%.

Figure 2 shows the $v_2$ results of prompt $J/\psi$ mesons at forward rapidities ($-2.86<y_{cm}<-1.86$ or $0.94<y_{cm}<1.94$) for high-multiplicity ($185\le N_{trk}^{offline}<250$) $p\text{Pb}$ collisions, covering a $p_T$  range from 0.2 to 10$\text{GeV}$. Results obtained separately for $J/\psi$ meson rapidity in the Pb- and p-going direction are compared, and found to be consistent within statistical uncertainties. Thus, as mentioned earlier, combined $v_2$ values are presented for the best statistical precision. The $v_2$ results for $K_S^0$ and $\Lambda$ hadrons (light, strange-flavor), and prompt $D^0$ mesons (open heavy-flavor), reported in a previous CMS publication [43] for the midrapidity region $-1.46<y_{cm}<0.54$, are also shown for comparison.

Positive prompt $J/\psi$ meson $v_2$ values are observed over a wide $p_T$ range from about 2 to 8$\text{GeV}$. The prompt $J/\psi$ meson $v_2$ results show a trend of first increasing up to $p_T\sim 4$$\text{GeV}$ and then decreasing toward higher $p_T$. This observed trend appears to be in common with the other hadron species shown. In the $p_T$ range below 5$\text{GeV}$, the $v_2$ values for $J/\psi$ and $D^0$ mesons are consistent with each other within the uncertainties, while an indication of smaller $v_2$ values for $J/\psi$ mesons than that for $D^0$ mesons is seen for $p_T>5$$\text{GeV}$. Over the full $p_T$ range, the $v_2$ signal values for both $J/\psi$ and $D^0$ hadrons are smaller than those for $K_S^0$ and $\Lambda$ hadrons. This observation is consistent with the earlier conclusion that charm quarks develop a weaker collective dynamics than light quarks in small systems [43], unlike what is seen in AA collisions. Because of experimental limitation, $v_2$ values for the prompt $J/\psi$ meson and the other meson species are not compared within the same rapidity range, possibly affecting their comparison. The rapidity dependence of $v_2$ values for charged particles in $p\text{Pb}$ collisions has been measured [59, 60], suggesting up to around 15% variation from $\left|y_{\text{lab}}\right|\sim 0$ to $2.4$.

To better study the elliptic flow signal coming purely from long-range collective correlations, the $J/\psi$ $v_2$ results are corrected for residual jet correlations. The resulting ($v_2^{\text{sub}}$) values are shown in Fig. 3 (upper) for prompt $J/\psi$ mesons as a function of $p_T$ with $185\le N_{trk}^{offline}<250$, and compared to similarly corrected $K_S^0$, $\Lambda$, and $D^0$ hadron results [43]. The effect of the correction for all particle species is most noticeable at very high $p_T$, while the overall $p_T$ dependence of the $v_2$ data remains unchanged. The $K_S^0$ mesons have a larger correction applied to their $v_2$ values (possibly because $K_S^0$ mesons are more correlated with the bulk multiplicity, and thus are biased toward stronger jet correlations due to the selection of high multiplicities) and their $v_2^{\text{sub}}$ values after the correction tend to converge to those of the prompt $J/\psi$ and $D^0$ mesons at high $p_T$.

A recent model calculation of $J/\psi$ $v_2$ in $p\text{Pb}$ collisions, based on final-state interactions between produced charm quarks and a QGP medium, suggests a very small $v_2$ signal of less than 0.01 [45]. This calculation indicates that additional contributions, e.g., those from initial-state interactions, may be needed to account for the observed $v_2$ signal of prompt $J/\psi$ mesons reported in this Letter.

Motivated by the possible formation of a hydrodynamically expanding QGP medium in small systems, the elliptic flow signals for $K_S^0$, $\Lambda$, $J/\psi$ and $D^0$ hadrons are compared as a function of transverse kinetic energy ($K{E}_T$) in Fig. 3 (lower), to account for the mass difference among the four hadron species [61, 62]. Here, the values of $v_2^{\text{sub}}$ and $K{E}_T$ are both divided by the number of constituent quarks, $n_q$, to represent the collective flow signal at the partonic level in the context of the quark coalescence model [63, 64, 65], which postulates that the elliptic flow signal of a hadron is a sum of contributions from individual constituent quark flow values. As was previously reported in $p\text{Pb}$ collisions [27, 43], a scaling of $n_q$-normalized $v_2^{\text{sub}}$ values is observed between the $K_S^0$ meson and $\Lambda$ baryon, shown in Fig. 3 (lower). This scaling between light baryon and meson species systems produced in the collision (known as the number-of-constituent-quark or NCQ scaling) was first discovered in AA colliding systems [66, 61, 62], indicating that collectivity is first developed among the partons, which later recombine into final-state hadrons. The elliptic flow signal per quark ($v_2^{\text{sub}}/n_q$) for prompt $J/\psi$ mesons at low $K{E}_T/n_q$ range is consistent with those of $K_S^0$, $\Lambda$, and prompt $D^0$ hadrons within large statistical uncertainties for the current data. There is a hint that the prompt $J/\psi$ meson data tend to fall on the same trend as those of $K_S^0$ and $\Lambda$ baryons, all of which are above the prompt $D^0$ meson data. A more definitive conclusion could be drawn with future high precision data. For $K{E}_T/n_q>1\text{GeV}$, the $v_2^{\text{sub}}/n_q$ for prompt $D^0$ and $J/\psi$ mesons are consistently below that of the $K_S^0$ meson. An indication of smaller $v_2^{\text{sub}}/n_q$ values for $J/\psi$ mesons than for $D^0$ mesons is seen for $K{E}_T/n_q>1\text{GeV}$. As $J/\psi$ mesons contain two charm quarks, while $D^0$ mesons contain a charm and a light-flavor quark, this observation would be consistent with a weaker collective behavior of heavy-flavor quarks than light quarks, possibly a consequence of the much smaller size of the collision system. Future data with improved precision will provide crucial insights to fully constrain the collective behavior of light- and heavy-flavor quarks in high-multiplicity, small systems.

In summary, the elliptic flow harmonic ($v_2$) for prompt $J/\psi$ mesons in high-multiplicity proton-lead ($p\text{Pb}$) collisions at $\sqrt{s_{{}_{NN}}}=8.16\text{TeV}$ is presented as a function of transverse momentum ($p_T$). Positive $v_2$ values are observed for prompt $J/\psi$ mesons at forward rapidity ($-2.86<y_{cm}<-1.86$ or $0.94<y_{cm}<1.94$) over a wide $p_T$ range ($2<p_T<8$$\text{GeV}$). This observation provides evidence for charm quark collectivity in high-multiplicity $p\text{Pb}$ collisions, similar to that first observed for light-flavor hadrons. The observed ordering of $v_2$ among light-flavor, open and hidden heavy-flavor hadrons at intermediate and high-$p_T$ regions (e.g., above 4$\text{GeV}$) adds support to the earlier conclusion that heavy quarks exhibit weaker collective behavior than light quarks or gluons in small systems, unlike what is found in AA collisions. For particle transverse kinetic energy per constituent quark values less than 1$\text{GeV}$, the $v_2$ of prompt $J/\psi$ mesons is consistent with prompt $D^0$, $K_S^0$ and $\Lambda$ hadrons, within current uncertainties. A model calculation based on final-state interactions between charm quarks and a QGP medium in $p\text{Pb}$ collisions significantly underestimates the measured prompt $J/\psi$ $v_2$ signal. The new prompt $J/\psi$ meson results, together with previous results for light-flavor and open heavy-flavor hadrons, provide novel insights into the dynamics of the heavy quarks produced in small systems that lead to high final-state multiplicities.