Imprints of fluctuating proton shapes on flow in proton-lead collisions at the LHC

TL;DR

This paper investigates how fluctuating proton shapes influence flow observables in proton-lead collisions at the LHC, using a combined classical Yang-Mills and hydrodynamic model constrained by HERA data.

Contribution

It demonstrates the critical role of sub-nucleon scale proton shape fluctuations in describing flow measurements, with shapes constrained by experimental data rather than free parameters.

Findings

Proton shape fluctuations significantly affect azimuthal harmonic coefficients.

Inclusion of temperature-dependent viscosities improves agreement with data.

Results include detailed flow coefficients and HBT radii for charged and identified particles.

Abstract

Results for particle production in root-s = 5.02 TeV p+Pb collisions at the Large Hadron Collider within a combined classical Yang-Mills and relativistic viscous hydrodynamic calculation are presented. We emphasize the importance of sub-nucleon scale fluctuations in the proton projectile to describe the experimentally observed azimuthal harmonic coefficients, demonstrating their sensitivity to the proton shape. We stress that the proton shape and its fluctuations are not free parameters in our calculations. Instead, they have been constrained using experimental data from HERA on exclusive vector meson production. Including temperature dependent shear and bulk viscosities, as well as UrQMD for the low temperature regime, we present results for mean transverse momenta, harmonic flow coefficients for charged hadrons and identified particles, as well as Hanbury-Brown-Twiss radii.

Figures

Equations

$$T_q({\mathbf{b}}_{\perp })=\frac{1}{2\pi B_q}{e}^{-{{\mathbf{b}}_{\perp }}^2/(2B_q)},$$

$${\pi }^{\mu \nu }=T_{\mathrm{CYM}}^{\mu \nu }-\frac{4}{3}\epsilon u^{\mu }{u}^{\nu }+\frac{\epsilon }{3}{g}^{\mu \nu }.$$

$$(\eta /s)(T)$$

$$+b(T-T_c)\theta (T-T_c),$$

$$v_n\{2\}=\frac{1}{⟨N(N-1){⟩}_{\mathrm{ev}}}(⟨\mathrm{Re}\{Q_n{Q}_n^{\ast }\}-N{⟩}_{\mathrm{ev}}),$$

$$C(\mathbf{q})$$

Full text

Introduction In heavy ion collisions strong final state interactions that are well described within the framework of relativistic hydrodynamics produce azimuthally anisotropic particle distributions that are strongly correlated with the initial event geometry Heinz and Snellings (2013); Gale et al. (2013a); de Souza et al. (2015); Song et al. (2017). Measurements of produced hadrons in collisions of protons and other small nuclei with heavy ions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) have observed very similar azimuthal anisotropies as in heavy ion collisions Adare et al. (2013); Chatrchyan et al. (2013a); Abelev et al. (2013a); Aad et al. (2012) (also see the review Dusling et al. (2016)). Naturally, it was suggested that these anisotropies in small systems could have the same origin as in heavy ion collisions, namely the hydrodynamic response of the produced medium to the initial shape of the interaction region in the transverse plane Bozek (2012); Bozek and Broniowski (2013a, b, c); Bozek et al. (2013); Qin and Muller (2014); Werner et al. (2014); Kozlov et al. (2014); Romatschke (2015); Shen et al. (2017); Weller and Romatschke (2017).

It has also been shown that in high energy collisions of hadrons and nuclei, gluons are produced with intrinsic anisotropic azimuthal momentum distributions (for recent reviews see Dusling et al. (2016); Schlichting and Tribedy (2016)). These distributions have some features often attributed to hydrodynamic behavior Schenke et al. (2016); Dusling et al. (2017) but are uncorrelated with the global geometry of the interaction region Schenke et al. (2015); McLerran and Skokov (2017). Furthermore, it is not clear whether they will survive potential final state interactions. In this work we focus on high multiplicity events and assume that the initial state momentum correlations are fully erased by strong final state interactions. Thus anisotropies in final observables are enitrely due to the system’s response to the initial geometry. We further assume that relativistic fluid dynamics is applicable in the presence of fairly large gradients (as present in small collision sytems such as p+Pb) and at times $\lesssim 1\,{\rm fm}$, something which has been extensively argued to be the case in the literature Chesler and Yaffe (2010); Wu and Romatschke (2011); Heller et al. (2012); van der Schee (2013); Casalderrey-Solana et al. (2013); Kurkela and Zhu (2015); Keegan et al. (2016a); Chesler (2016); Attems et al. (2017); Romatschke (2017a, b).

The initial state and early time evolution in this work is described by the IP-Glasma model Schenke et al. (2012a, b). This model has been shown to lead to very good agreement with a wide range of experimental data Gale et al. (2013b); Schenke and Venugopalan (2014), and, in contrast to some other models, was shown to be compatible with data in an extensive Bayesian parameter estimation study Bernhard et al. (2016). The only other models that are in good agreement with a similar amount of experimental data, including event-by-event distributions of flow harmonics Gale et al. (2013b); Niemi et al. (2013); Schenke and Venugopalan (2014); Aad et al. (2013), are the EKRT model Niemi et al. (2015), and a particular version of TRENTo Moreland et al. (2015).

The field energy-momentum tensor from the IP-Glasma model provides the initialization for the relativistic viscous hydrodynamic simulation Music Schenke et al. (2010, 2011, 2012c). In the last step of the calculation we determine spatial and momentum distributions of all particles and feed them into the hadronic cascade UrQMD Bass et al. (1998); Bleicher et al. (1999) to obtain our final results.

Hydrodynamic evolution of IP-Glasma initial states in small systems has been explored before Bzdak et al. (2013); Schenke and Venugopalan (2014). The most notable result from the analysis in Schenke and Venugopalan (2014) was that the measured harmonic flow coefficients in p+Pb collisions were underestimated significantly. This stood in sharp contrast with calculations using the Monte-Carlo (MC) Glauber model for the initial state, which lead to much better agreement with experimental flow harmonics in small systems Bozek (2012); Bozek and Broniowski (2013a, b, c); Bozek et al. (2013); Qin and Muller (2014); Werner et al. (2014); Kozlov et al. (2014); Romatschke (2015); Shen et al. (2017); Weller and Romatschke (2017). As anticipated in Schenke and Venugopalan (2014) the reason for the disagreement with experimental data is the underestimation of subnucleonic fluctuations in the IP-Glasma model - the profile of a nucleon’s color charge density was assumed to be Gaussian in the transverse plane (called ’round proton’ in the following). In the Yang-Mills framework, the interaction region follows much more closely the shape of the smaller projectile than in MC Glauber implementations. This is why in p+A collisions one becomes sensitive to the assumed proton shape and its fluctuations, while in d/3He/A+A collisions spatial eccentricities are dominated by nucleon size scale fluctuations.

Because fluctuations of the proton shape constitute an additional degree of freedom, their inclusion requires an additional external constraint. In the spirit of the IP-Glasma model, whose parameters are determined by HERA deeply inelastic scattering data Rezaeian et al. (2013), we employ additional data from HERA on diffractive $J/\Psi$ production Breitweg et al. (2000); Chekanov et al. (2002, 2003); Aktas et al. (2003, 2006); Alexa et al. (2013). In a recent work by two of the authors the incoherent diffractive cross section for exclusive $J/\Psi$ production was used to constrain the degree of fluctuations of the gluon distribution in the proton within the IP-Glasma model Mäntysaari and Schenke (2016a, b). In this work we employ the model using three gluonic hot spots Schlichting and Schenke (2014) with parameters determined in Mäntysaari and Schenke (2016b). Similar models were discussed recently in Albacete and Soto-Ontoso (2016); Welsh et al. (2016); Cepila et al. (2017); Weller and Romatschke (2017); Mäntysaari and Schenke (2017); Moreland et al. (2017).

We show in the following that strong final state interactions turn initial shape fluctuations of the proton’s gluon distribution into observable final state particle correlations. The values of the resulting flow harmonics differ from those obtained for round protons by up to a factor of 5. This means that high energy p+A collisions provide unprecedented access to analyze the shape of the proton and its fluctuations. The constrained proton shape parameters from Mäntysaari and Schenke (2016a, b) together with temperature dependent medium properties similar to those employed in a previous work Denicol et al. (2016) allow for a quantitative description of experimentally measured particle spectra, flow harmonics, and Hanbury-Brown-Twiss (HBT) radii in p+Pb collisions at $5.02\,{\rm TeV}$.

Fluctuating proton in the IP-Glasma initial state As discussed above, the original IP-Glasma model assumed a 2D-Gaussian spatial shape for all nucleons. This assumption lead to a dramatic underestimation of flow harmonics in p+A collisions Schenke and Venugopalan (2014). Here we introduce an additional substructure to all nucleons, using three hot spots whose positions in the transverse plane are Gaussian distributed with width $B_{qc}$. The density profile of each hot spot in the transverse plane is also assumed to be Gaussian

[TABLE]

with width parameter $B_{q}$. Other than this substitution for the round nucleon, the implementation of the IP-Glasma model is exactly as described in Schenke et al. (2014). The infrared regulator $m$ is determined together with $B_{qc}$ and $B_{q}$ from the analysis of HERA data in Mäntysaari and Schenke (2016b). The values are $m=0.4\,{\rm GeV}$, $B_{qc}=3\,{\rm GeV}^{-2}$, and $B_{q}=0.3\,{\rm GeV}^{-2}$.

Because of the shorter lifetime in small systems, the initial viscous stress tensor from the classical Yang-Mills (CYM) simulation plays a more important role compared to heavy ion collisions. We thus include the full $T^{\mu\nu}_{\rm CYM}$ when initializing the hydrodynamic simulation as opposed to only extracting the energy density $\varepsilon$ and transverse flow velocities $u^{\mu}$ as done in previous works. Using the ideal equation of state $\varepsilon=3P$ of the classical Yang-Mills system we have

[TABLE]

We find the spatial dependence of components of $\pi^{\mu\nu}$ to be similar to the initial Navier-Stokes value when using a shear viscosity of $\eta/s=0.1$. We further include a correction for the initial value for the bulk component of the stress energy tensor to account for the different equation of state on the hydrodynamic side of the simulation, which is taken from lattice QCD Huovinen and Petreczky (2010). This matching is performed at the initial switching time $\tau_{0}$, which we vary from $\tau_{0}=0.2\,{\rm fm}$ to $\tau_{0}=0.4\,{\rm fm}$.

Hydrodynamics and hadronic cascade In the hydrodynamic simulation Music we use the same second order transport parameters as in Ryu et al. (2015). We employ a temperature dependent shear viscosity to entropy density ratio with similar features as those found to describe rapidity dependent flow harmonics at RHIC energies in Denicol et al. (2016). These features are a strong increase with decreasing temperature in the hadronic phase and no to little temperature dependence in the QGP phase. To be precise, we use the following parametrization

[TABLE]

where $a=15$ GeV*-1*, $b=1$ GeV*-1*, and $(\eta/s)_{\rm min}=0.08$ at a temperature of $166\,{\rm MeV}$. We note that both the minimal value and the slope in the hadronic phase ($a$) are larger than what was used in Denicol et al. (2016). The reason for needing larger viscosities to describe the p+Pb data could be that i) we are using a different initial state model and ii) we use parameters for the initial proton shape determined at $x\approx 10^{-3}$, while typical $x$ values in 5.02 TeV p+Pb collisions at midrapidity are $x\approx 2\cdot 10^{-4}$. Explicit small $x$ evolution via solution of the JIMWLK equation Jalilian-Marian et al. (1997, 1999); Iancu et al. (2001) that we do not include here would lead to the reduction of eccentricities towards this lower $x$ Schlichting and Schenke (2014), leading in turn to an extraction of smaller viscosity values. The inclusion of this energy dependence is numerically intensive but will be a very interesting project for the future.

We also employ an effective shear viscosity of $\eta/s=0.2$ for comparison. The bulk viscosity to entropy density ratio $\zeta/s$ used is similar to that in Denicol et al. (2009, 2016); Ryu et al. (2015), the difference being an exponential reduction of $\zeta/s$ at low temperatures to keep non-equilibrium corrections to the thermal distribution functions Rose et al. (2014) under control. We note that in small systems the effect of bulk viscosity on the evolution is essential to reproduce the $\langle p_{T}\rangle$ of identified hadrons, even more so than in heavy ion collisions Ryu et al. (2015). A potential problem is that bulk viscosity in combination with the large expansion rates in p+Pb collisions can lead to negative effective pressures. We do not stop the hydrodynamic evolution when such negative values appear, whose existence has lead to the discussion of the potential phenomenon of cavitation Rajagopal and Tripuraneni (2010); Bhatt et al. (2011); Habich and Romatschke (2014); Denicol et al. (2015); Sanches et al. (2015). We note that anisotropic hydroynamics Nopoush et al. (2014); Alqahtani et al. (2015) avoids (or at least reduces in a more general implementation) negative pressures by resumming certain viscous corrections.

At a temperature of $T_{\rm switch}=155\,{\rm MeV}$ we determine the switching surface and sample particles that then propagate in the hadronic cascade model UrQMD. Particle momentum distributions are sampled first using the probability distribution function obtained from integrating the Cooper-Frye formula over the whole hyper-surface. Then the spatial positions of the sampled particles are determined from the differential Cooper-Frye formula. This sampling procedure introduces the least distortion in the momentum distribution of particle samples when the Cooper-Frye formula takes on negative values in certain regions of the hypersurface Shen et al. (2016).

We have checked by explicit calculation that for charged hadrons the effect of rescattering in UrQMD is negligible. Proton spectra and $v_{n}(p_{T})$ are slighty blue shifted because of the additional hadronic scatterings Shen et al. (2017).

Results We begin by presenting results for the average transverse momentum $\langle p_{T}\rangle$ of identified particles as a function of charged particle multiplicity in Fig. 1. Using both temperature dependent $\eta/s$ and $\zeta/s$ and a switching time of $\tau_{0}=0.4\,{\rm fm}$ we find good agreement with experimental data from the ALICE collaboration for charged pions, protons, and $\Lambda$’s. The $\langle p_{T}\rangle$ of charged kaons is underestimated. We note that as discussed for heavy ion collisions in Ryu et al. (2015), the inclusion of bulk viscosity is essential in order not to overestimate $\langle p_{T}\rangle$. Without bulk viscosity, the pion $\langle p_{T}\rangle$ is overestimated the most, by approximately $50\%$. The effect of using the constant effective $\eta/s$ is weak as is the effect of a smaller switching time $\tau_{0}=0.2\,{\rm fm}$, which is not shown here.

Having established the agreement with measured transverse momentum spectra, which is almost entirely determined by $\langle p_{T}\rangle$, we now present results for $v_{n}$ from two-particle correlations. To compute $v_{n}\{2\}$ using particle samples from UrQMD, we first construct the flow vector $Q_{n}=\sum_{i}w_{i}e^{in\phi_{i}}$, where the sum $i$ runs over all particles of interest with $0.3\,{\rm GeV}<p_{T}<3\,{\rm GeV}$ (when comparing to CMS results), and the weights are set to $w_{i}=1$. The two particle cumulant $v_{n}\{2\}$ is then computed as

[TABLE]

where $N$ is the number of particles included in the calculation of $Q_{n}$ and $\langle\cdot\rangle_{\rm ev}$ is the average over events. In practice, we sample the hypersurface from each hydrodynamic event 5000 times and run UrQMD for each of these particle configurations. For the evaluation of $v_{n}\{2\}$ we combine the UrQMD output of all 5000 runs to collect enough statistics and suppress short range correlations from e.g. resonance decays. The latter effect is desired because the measurement uses a large pseudo-rapidity gap of $|\Delta\eta|>2$ between the two particles, also eliminating short range correlations.

In Fig. 2 we see that above a multiplicity of $N_{\rm trk}^{\rm offline}\simeq 80=2\langle N_{\rm trk}^{\rm offline}\rangle$, the experimental $v_{n}\{2\}$ are well reproduced by our calculation, considering the uncertainty from the initial switching time, which we vary from $\tau_{0}=0.2\,{\rm fm}$ to $\tau_{0}=0.4\,{\rm fm}$.111To represent our results as a function of the experimentally determined $N_{\rm trk}^{\rm offline}$, we rescaled our computed multiplicities by a factor $\langle N_{\rm trk}^{\rm offline}\rangle/\langle dN_{\rm ch}/d\eta\rangle$. Below $N_{\rm trk}^{\rm offline}\simeq 80=2\langle N_{\rm trk}^{\rm offline}\rangle$ the $v_{n}\{2\}$ are overestimated. This could indicate the failure of the hydrodynamic framework at low multiplicity. However, the experimental procedure to subtract correlations from peripheral events may cause a too fast decrease of $v_{n}\{2\}$ with decreasing multiplicity.

We repeat the same comparison in Fig. 3 using the constant effective $\eta/s=0.2$. We do not find a large difference to results with temperature dependent $\eta/s$ as long as the temperature does not rise too rapidly in the QGP phase. A more rapid rise ($b>2\,{\rm GeV}^{-1}$) was already ruled out in Denicol et al. (2016) and is not considered here. As was already visible in Fig. 2, the earlier switching to hydrodynamics leads to an increased anisotropic flow. This dependence on the switching time indicates that a more careful treatment of the non-equilibrium early time evolution is required.

In Fig. 4 we present our results for the transverse momentum dependent $v_{2}\{2\}(p_{T})$ and $v_{3}\{2\}(p_{T})$ in one multiplicity class and compare to experimental data from the ATLAS collaboration Aad et al. (2014). Here, the reference bin for the second particle is $1\,{\rm GeV}<p_{T}<3\,{\rm GeV}$. We find excellent agreement with the experimental data for the temperature dependent $\eta/s$. Results for the effective $\eta/s=0.2$ are slightly lower (not shown) but agree within statistical errors.

In Fig. 5 we present the transverse momentum dependent $v_{2}$ of identified pions and protons in $0-20\%$ central events. We compare with experimental data from the ALICE experiment Abelev et al. (2013b), where the reference $p_{T}$ bin is $0.3\,{\rm GeV}<p_{T}<4\,{\rm GeV}$. The characteristic mass dependence of $v_{2}(p_{T})$ is well reproduced for $p_{T}<2\,{\rm GeV}$ and when comparing to results with correlations from 60-100% central events subtracted.

In Fig. 6 we present the $v_{2}(p_{T})$ of $\Lambda$ baryons and $K^{0}_{S}$. Again, a clear mass dependence is visible. Agreement with experimental data from the CMS collaboration Khachatryan et al. (2015) is very good. The results in Figs. 5 and 6 demonstrate that the mass ordering of $v_{2}$, while also possible to obtain from initial state momentum correlations and string fragmentation Schenke et al. (2016), is quantitatively reproduced in our hydrodynamic simulations of p+Pb collisions.

Finally we present results for the HBT radii $R_{\rm long}$, $R_{\rm side}$, and $R_{\rm out}$, using both UrQMD results that include hadronic rescattering, and a calculation, independent of UrQMD, that includes only resonance decays, similar to Therminator Chojnacki et al. (2012). For a given range in particle pair momentum $\mathbf{k}_{T}=(\mathbf{p}_{T}^{1}+\mathbf{p}_{T}^{2})/2$, and under the assumption of a plane wave superposition for the pion wave function, the two-particle correlation function can be written as

[TABLE]

where $ij$ runs over all particle pairs $(i,j)$ within the $k_{T}=|\mathbf{k}_{T}|$ bin, and $q^{\mu}_{ij}=p^{\mu}_{i}-p^{\mu}_{j}$ and $x^{\mu}_{ij}=x_{i}^{\mu}-x_{j}^{\mu}$, with $x^{\mu}_{i}$ the space-time four vector of the last interaction point of particle $i$. $\langle N_{\rm pair}\rangle$ is the average total number of all pairs from single events, $\langle N_{\mathrm{mixpair}}\rangle$ the average total number of pairs from different (hydrodynamic) events, and $\langle N_{\mathrm{mixpair}}(q)\rangle$ is the average number of pairs from different events in the considered $q$-range. For this analysis we also perform many UrQMD runs for a single hydrodynamic event. This results in values of $\langle N_{\rm pair}\rangle$ and $\langle N_{\rm mixpair}\rangle$ of the order of $10^{8}$ per $k_{T}$ bin.

The constructed two-pion HBT correlation function is fitted to the Pratt-Bertsch parameterization in the longitudinally co-moving system (LCMS) Pratt (1986); Bertsch et al. (1988). We note that the shape of $C(\mathbf{q})$ is rather different from a Gaussian once resonance decay contributions are included. This non-Gaussianity results in different HBT radii depending on the fit range Plumberg and Heinz (2016). This is illustrated by the shaded bands in Fig. 7, where we show $R_{\rm long}$, $R_{\rm side}$, and $R_{\rm out}$ along with the ratio $R_{\rm out}/R_{\rm side}$ as functions of $k_{T}$ for temperature dependent $\eta/s$. For the radii, the upper end of the bands results from fitting to the smallest momentum range $0.05\,{\rm GeV}<|\mathbf{q}|<0.08\,{\rm GeV}$, the lower end to the largest, $0.05\,{\rm GeV}<|\mathbf{q}|<0.2\,{\rm GeV}$. The systematic error quoted by the ALICE collaboration includes the variation of the upper limit of the fit range from $0.3\,{\rm GeV}$ to $1.1\,{\rm GeV}$. The limited statistics at large $|\mathbf{q}|$ in our calculation constrains the possible fit range to smaller values.

We find that for $R_{\rm out}$ the lower end of both bands (dashed lines for the case of UrQMD including scattering in the later hadronic phase, solid lines for the case of only including decays) are compatible with the experimental data from the ALICE Collaboration Adam et al. (2015). The other two radii are overestimated by both calculations, resulting also in a slight underestimation of the ratio $R_{\rm out}/R_{\rm side}$, particularly in the case of using UrQMD. A smoother initial state (e.g. from small $x$ evolution) and reduced bulk viscosity could overcome this tension. Results with the constant effective $\eta/s$ are almost identical to the ones shown.

These results are comparable to previously published HBT radii from calculations using various other initial state models Bozek and Broniowski (2013b); Shapoval et al. (2013); Romatschke (2015). In our calculation, the smaller (compared to e.g. the MC-Glauber model Bozek and Broniowski (2013b)) initial size of the system is overcome by a more explosive expansion and the accompanying extended lifetime due to entropy production from bulk viscous effects, leading to HBT radii being slightly larger than the experimental measurement.

Discussion Early-time non-equilibrium phase: Given the shorter life-time compared to heavy ion collisions, small systems such as p+Pb collisions are more sensitive to the early time non-equilibrium stage of the evolution. We have performed calculations that initialized the hydrodynamic simulation directly with the energy momentum tensor from the classical Yang-Mills IP-Glasma calculation at times $\tau_{0}=0.2\,{\rm fm}$ and $0.4\,{\rm fm}$. All deviations from equilibrium and isotropy are thereby absorbed into viscous corrections on the hydrodynamic side of the simulation. When not including the initial viscous stress tensor, values of integrated $v_{n}$ increase by up to 50% (for $\tau_{0}=0.4\,{\rm fm}$, $\eta/s=0.2$ for $3-3.75\langle N_{\rm trk}^{\rm offline}\rangle$) underlining the importance of early times in p+Pb collisions. So it is desirable to improve on this prescription in the future by including a realistic dynamic description of the non-equilibrium stage, which will smoothly interpolate between the IP-Glasma initial state and the viscous hydrodynamic description, without requiring very large initial viscous shear stress tensors in the hydrodynamic stage. A first attempt at this, using an effective kinetic theory, was presented in Keegan et al. (2016b).

Off-equilibrium corrections: Surprisingly, effects of off-equilibrium corrections to the thermal distribution functions largely cancel in final results for $p_{T}$-spectra or $v_{n}(p_{T})$, however, studying the distribtution of corrections relative to the thermal distribution from all freeze-out surface cells, we find a significant share (10% for $p_{T}\sim 0.45\,{\rm GeV}$, 25% for $p_{T}\sim 1\,{\rm GeV}$, 45% for $p_{T}\sim 1.5\,{\rm GeV}$ for pions) of large (shear) corrections ($|\delta f|/f\gtrsim 100\%$). This demonstrates that our results are plagued by large viscous corrections, in particular for $p_{T}\gtrsim 1\,{\rm GeV}$. Nevertheless, $p_{T}$-integrated quantities are dominated by low $p_{T}$ contributions and less sensitive to these problems.

Boost-invariance: The IP-Glasma model is by construction boost invariant. It will be very interesting to use an extension of this model discussed in Schenke and Schlichting (2016) to study p+Pb collisions within 3+1 dimensional hydrodynamic simulations. We note, however, that earlier work using a MC-Glauber initial state has found that at $\sqrt{s}=5.02\,{\rm TeV}$ differences of midrapidity $v_{n}$ between full 3+1 dimensional and boost invariant 2+1D simulations are negligible Shen et al. (2017).

Conclusions We have presented results of particle production, anisotropic flow, and HBT radii in $\sqrt{s}=5.02\,{\rm TeV}$ proton lead collisions at the LHC from a hydrodynamic framework. We use an initial state model that is based on an effective theory of QCD and compatible with a wide range of heavy ion data (IP-Glasma). This is coupled to a relativistic viscous hydrodynamic simulation (Music), which finally is connected to a microscopic hadronic afterburner (UrQMD).

The main new ingredient in this work is the inclusion of subnucleonic fluctuations. These were independently shown to be highly relevant by comparison of the IP-Glasma calculation of diffractive $J/\Psi$ production with HERA data, which also allowed us to constrain the degree of subnucleonic fluctuations quantitatively.

We have demonstrated that including subnucleonic fluctuations increases flow harmonics in p+Pb collisions by approximately a factor of 5 compared to previous calculations using round nucleons Schenke and Venugopalan (2014). Furthermore, we have shown that the multiplicity dependence of transverse momentum spectra, $v_{2}$ and $v_{3}$, as well as the transverse momentum dependence of the charged hadron and identified particle flow harmonics is well described within our framework, at least for multiplicities above $\sim 2$ times the average multiplicity. For lower multiplicities, we expect initial state momentum correlations (see Dusling et al. (2016); Schlichting and Tribedy (2016)) to become more dominant, and more sensitivity to the experimental method for removing jet-like correlations (new methods show little decrease of e.g. $v_{2}\{4\}$ with decreasing multiplicity Jia et al. (2017)).

With the exception of $R_{\rm out}$, HBT radii are overestimated. However, a direct comparison with experimental data is difficult because of the non-Gaussian shape of the experimental and theoretical correlation functions, and limited statistics of the calculation at large momentum differences.

Future analyses with much higher statistics could be used to constrain the fine details of the subnucleonic structure in the initial state, in particular when higher order cumulants and other more complex multi-particle correlators can be studied. With that we could get unprecedented access to the fluctuating shape of protons and its energy dependence via proton-nucleus collisions.

Acknowledgments We thank Christopher Plumberg for useful discussions about HBT correlations. The authors are supported under DOE Contract No. DE-SC0012704. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. BPS acknowledges a DOE Office of Science Early Career Award.