Hanbury Brown--Twiss Interferometry and Collectivity in Small Systems

Hanbury Brown--Twiss interferometry (HBT) provides crucial insights into both the space-time structure and the momentum-space evolution of ultrarelativistic nuclear collisions at freeze-out. In particular, the dependence of the HBT radii on the transverse pair momentum $K_T$ and the system charged multiplicity $dN_{\mathrm{ch}}/d\eta$ may reflect the mechanisms driving collective behavior in small systems. This paper argues that certain features observed in the multiplicity dependence of the HBT radii can be naturally understood if small systems evolve hydrodynamically at high-multiplicity. This study thus establishes a baseline for the multiplicity dependence of HBT in hydrodynamics which may prove useful in discriminating between competing models of collectivity in nuclear collisions.


I. INTRODUCTION
The existence of collective, fluid-like behavior in relativistic nuclear collisions, from pp to A-A , is by now well established [1,2]. Understanding the precise origins of this collective behavior, however, remains one of the foremost outstanding challenges in the field. To date, a number of explanations of this phenomenon have been proposed, including CGC-type models with collectivity built into the initial state [3,4], "escape mechanism" models which effectively generate collectivity kinematically [5], approaches based on string hadronization models [6,7], "one-hit" dynamical models [8], and relativistic hydrodynamics [9][10][11]. 1 The ability to discriminate between competing models of collectivity is therefore urgently needed and requires both quantitative predictions and comparison with experiment. In this context, femtoscopic observables, such as those derived from Hanbury Brown-Twiss (HBT) interferometry, offer a powerful and complementary glimpse into the space-time structure and dynamical evolution of nuclear collisions at freeze-out [13]. The most widely used of these observables, the "HBT radii," reflect collective effects in a number of ways, particularly in their dependence on the transverse pair momentum and on the system's charged multiplicity [14][15][16][17]. For instance, when comparing large and small systems at fixed multiplicity, initially more compact systems (pp , p-A ) need to develop stronger transverse flow in order to reach the same final freeze-out volume as attained in larger systems [17]. This enhanced flow, which is driven by the larger initial density gradients present in small systems [18], is a direct prediction of hydrodynamics and leads to a measurable ordering pp < p-A < A-A in the radii extracted from different systems at the same multiplicity. This study will explore the implications of the enhanced radial flow produced by hydrodynamics in small systems for the multiplicity dependence of the HBT radii.
The multiplicity dependence of HBT in small and large systems has already been studied in a fair number of experimental analyses [19][20][21][22]. One notable feature of these measurements is that the collision system's volume, when estimated from the HBT radii, appears to scale linearly with the charged particle multiplicity dN ch /dη, while the individual radii each scale linearly with (dN ch /dη) 1/3 . Moreover, the dependence of the individual radii in large systems is seen to fall roughly onto a single universal, approximately linear trajectory in (dN ch /dη) 1/3 , regardless of collision species [13]. This is exactly what one would expect to find if pp , p-A , and A-A collisions are all driven by hydrodynamics.
What is initially surprising, then, is to find in the data that the individual radii across collision systems should differ not only in their magnitudes, as implied by the enhanced flow in small systems, but also in the slopes of their respective (dN ch /dη) 1/3 dependences.
There are two specific features to be noted. First, in large systems, as noted above, each radius has a slope which is approximately independent of the collision species, whereas in pp and p-A , the corresponding slopes tend to be considerably smaller (the exception is the 'side' radius, as we will see below). I will refer to this feature of the data as the slope non-universality exhibited by the radii in small systems. In addition, not only do the slopes in small systems tend to deviate from the universal slope of large systems, but they also disagree more significantly amongst themselves: the rough ordering of the slopes within a fixed system is There is therefore also a slope hierarchy exhibited by the different radii which depends on which collision system is being considered. This hierarchy is clearly present in both large and small systems, but varies in strength between them. The presence of these features in the multiplicity dependence of HBT in pp , p-A , and A-A implies radical differences in the space-time evolution of large and small systems.
The preceding observations seem to introduce some unwelcome complexity into an otherwise simple situation. Taken together, the two features just identified -the non-universality and hierarchy of the slopes -appear to stand in tension with the intuitive expectation that hydrodynamics should lead systems of all sizes to evolve in similar ways. This in turn raises the crucial question of whether the features appearing in small systems are indeed signatures of genuine, hydrodynamic collectivity or of something else and, more generally, whether the mechanisms driving the dynamical evolution of small systems are the same as those at play in larger systems [23,24].
The goal of this paper is to show how the features just identified in the multiplicity dependence of the HBT radii arises naturally within the context of hydrodynamics. More precisely, I will show using a simplified hydrodynamic model that both the non-universality and the hierarchy exhibited by the slopes of the (dN ch /dη) 1/3 dependence of the HBT radii across different collision systems emerges naturally within the hydrodynamic paradigm, at least at sufficiently large multiplicities. While the simplifications used in my hydrodynamic modeling limit my discussion here to a somewhat qualitative level, they can (and should) be removed for a more quantitative interpretation of the experimental data in future work.
The outline of this paper is as follows. In Sec. II I review the basic elements of HBT interferometry and show how to justify the interpretation of the HBT radii in terms of the space-time structure of the freeze-out surface. In Sec. III, I take a closer look at the data which most clearly illustrate the significant differences in the radii and their multiplicity scaling when compared across various collision systems. In Sec. IV I will show using a highly simplified hydrodynamic model how the discrepancies observed in Sec. III might reflect fluid dynamical behavior in small collision systems. Finally, Sec. V will summarize the main results and offer some suggestions for future work which will flesh out these ideas in a more quantitative fashion.

II. FORMALISM
A. The correlation function HBT interferometry relies on the existence of Bose-Einstein or Fermi-Dirac correlations between pairs of identical particles. The techniques underlying HBT have been developed and reviewed extensively elsewhere [13,25]. Here, I will briefly present the essential elements which are necessary to establish my notation and to show how the space-time structure of the source may be inferred.
The basic observable of HBT interferometry is the two-particle correlation function, defined by (1) Ideally, it is constructed so as to reduce to unity in the absence of actual Bose-Einstein correlations between identical particle pairs -in our case, pairs of π + bosons -produced by the collision event. Theoretically, it is usually convenient to consider instead of (1) the equivalent correlation function evaluated in terms of the relative momentum q = p 1 − p 2 and the pair momentum K = (p 1 + p 2 )/2 and relate it directly to an "emission function" (or "source function") S(x, K): S can be thought of in essence as a quantum-mechanical phase space (or 'Wigner') distribution [26] which roughly characterizes the probability to emit a particle from position x with momentum K. The step from (2) to (3) makes use of the so-called "smoothness assumption" [27] which is well-justified for large sources such as A-A collisions, but becomes questionable in pp and p-A collisions. This assumption will be relaxed in the full analysis which follows, although its effects turn out to be mostly negligible when evaluated quantitatively.
In any event, the width of the correlation function in q reflects the space-time structure of the underlying source at a fixed K. This structure can be inferred by suitably parameterizing the correlation function with a functional form such as The R 2 ij ( K) and λ( K) are extracted as free parameters, obtained by fitting (4) to one of the theoretical correlation functions (2) or (3). The quantities R 2 ij ( K) are known as the HBT radii and quantify the space-time structure of the emitting source, and λ( K) is an ad hoc factor which typically deviates from unity when effects due to resonance decays [28] or coherent pion production [29][30][31][32] are important. These effects will be neglected in this study, meaning that λ( K) = 1. 2 In the special case of a perfectly Gaussian source (and making use of the smoothness assumption) [26], one can perform the Fourier integrals in (3) analytically, yielding an exact relation between the R 2 ij ( K) and space-time variances of the underlying source [33,34]: Here the averages are taken with respect to the emission function: Additionally, the pair velocity β is given by and the separate terms ( x ixj , . . .) comprising the righthand side of Eq. (5) are known as the "source variances" [35]. Although the "Gaussian source approximation" is not used in this analysis, it will be useful here in interpreting and developing intuition for the results presented below.
It is worth emphasizing that the relations (5) -(7) provide the essential connection between the R 2 ij and the spatial and temporal characteristics of the underlying source function, which justifies the usual interpretation of the HBT radii in terms of the space-time geometry of the collision system at freeze-out. Nevertheless, the radii do not reflect only spatial lengthscales in the system, but necessarily represent a mixture of spatial and temporal information together.

B. The emission function
To proceed further, we need to specify the emission function S which governs the particle production process in a nuclear collision. For the systems studied here using hydrodynamics, the emission function can be defined straightforwardly according to the standard Cooper-Frye prescription [36]: where u(x) is the local flow velocity profile and Σ(y) signifies the freeze-out surface over which the integral is evaluated. Viscous corrections to the distribution function have been neglected here for simplicity. This is reasonable, as the precise form and magnitude of these corrections are still not extremely well-constrained theoretically [37,38], and in any event have little effect on the qualitative behavior of the R 2 ij obtained from hydrodynamic simulations with smooth or event-averaged initial states [35]. 3

C. Initial Conditions and Hydrodynamics
Hydrodynamics of course requires initial conditions. The initial conditions for this analysis were generated using the MC-Glauber model [40], including fluctuations of both the nucleon positions and collision-by-collision multiplicity fluctuations [10]. Event-averaged initial conditions were generated for each system (pp , p-Pb , Pb-Pb ) in 10% centrality-class intervals (0 − 10%, . . . 90 − 100%) by cutting on the total initial entropy at mid-rapidity, as described also in [10]. The overall normalization for each system was adjusted so that the system yielded a benchmark value of the charged particle pseudorapidity density dN ch /dη, obtained from experimental measurements for that system in a given reference centrality class. The benchmark value of dN ch /dη in each system's respective reference class is given in Table I. Several higher centrality classes (0-1%, 0-0.1%, 0-0.01%, and 0-0.001%) were also generated for pp and p-Pb .
The hydrodynamic evolution was performed using the 2+1D iEBE-VISHNU package [10] with specific viscosities η/s = 0.08 and ζ/s = 0 and the s95p-v1 equation of state 3 The same is not necessarily true for the hydrodynamic simulations themselves, as viscous effects influence not only the amount of particle production (and therefore the final charged multiplicity) but also contribute to the transverse flow [39] and consequently affect the shape of the freeze-out surface as well. For these reasons, some viscous effects have been retained in the hydrodynamic simulations presented below, despite being excluded from the calculation of the HBT radii. [44]. Hydrodynamics was initialized at a proper time τ 0 = 0.6 fm/c without including any preequilibrium effects. This is already a significant assumption, especially for small systems: hydrodynamics is typically valid only after preequilibrium dynamics (e.g., [45]) have enabled the system to 'hydrodynamize' on a timescale τ hydro ∼ O(1/T ) set by the temperature, which varies with the size of the collision system [46]. However, the fact that this simplification omits an important source of transverse flow in nuclear collisions means it is likely to underestimate the effects on the HBT radii [47]. Since including preequilibrium flow would likely only strengthen the conclusions drawn in this work, it will be neglected here for simplicity.
Once it is initialized, the hydrodynamic phase evolves in the usual way and is terminated when the system has cooled to a freeze-out temperature of T fo = 120 MeV, which is a typical value (cf., e.g., [48]). This is done in lieu of terminating at a higher temperature and evolving subsequently with a hadron cascade. After freeze-out, particle yields are obtained by evaluating Eq. (9) numerically as an integral over the freeze-out surface.

The correlation functions (2)-(3) can be similarly evaluated in terms of Cooper-Frye-like
integrals over the freeze-out surface [49]. In this case, the correlation function (2) is first evaluated on a fixed grid of points in q, K T , and the transverse pair momentum angle Φ K .
For each K T and Φ K , it is then fit to Eq. (4), which gives a set of R 2 ij as functions of K T and Φ K . More details of the fitting procedure are described in [49]. No systematic (e.g., fit-range [50]) uncertainties have been assessed for the fits in this study.
Once the fit radii are obtained, they are averaged separately over their angular dependence, finally yielding them as functions of K T only. Since the radii are azimuthally averaged, they are basically insensitive to differences between the x and y directions in large systems which tend to change with centrality [51].
This highly simplified hydrodynamic model allows us to concentrate on the essential features of interest in this study, namely, the connection between the presence of enhanced flow in small systems and the resultant scaling of the R 2 ij with multiplicity.

III. EXPERIMENTAL CONTEXT
The experimental motivation for this study originated from two analyses alluded to previously which explored the multiplicity dependence of HBT in pp and p-Pb collisions [20,22]  In the slopes of the HBT radii vs. (dN ch /dη) 1/3 , one observes both non-universality (radii have different slopes in large vs. small systems) and hierarchy (different radii possess differing slopes in a given system). The STAR results were published in [52,53]. The ALICE results for pp , p-Pb , and Pb-Pb were given in Refs. [20,22,54]. The fit lines were added by the author to guide the eye.
See the text for further discussion. and compared it with similar measurements for larger systems (such as Au+Au [52,53] and Pb-Pb [54]). Several of these measurements are shown in Fig.1.
There are two noteworthy features in this data which have been already discussed extensively in the literature [19-22, 55, 56] and which were briefly described in the Introduction.
First, the radii across different systems exhibit discontinuities in their magnitudes: that is, the values of the radii appear to differ across collision systems at fixed multiplicity. Second, and somewhat related to the first point, the radii exhibit different slopes with (dN ch /dη) 1/3 in the various directions and collision systems. Fit lines (dashed) have been included in A more thorough inspection of the complete datasets in Fig. 1 shows that these features also depend strongly on the K T value for which they are plotted (this will be seen clearly below in Sec. IV).
It is important to emphasize here that the hierarchy and non-universality of the slopes in Fig. 1 are entirely independent concepts: one could have had completely universal slopes which exhibited a hierarchy (i.e., were independent of collision system, but differed for each radius), or one could as easily have had no hierarchy between the various radii, but a nonuniversal slope for each radius whose value depended on the collision system. In the present case, of course, a mixture of both is found in the data of Fig. 1. One finds in A-A , for instance, that the slope of R o is somewhat larger (∼ 0.7) than that of R s or R l (∼ 0.6), whereas for pp , the trend is reversed: the R o fit has a slope comparable to R l (∼ 0.3), while the R s slope is considerably steeper (∼ 0.5). This paper is an attempt to organize these various observations within a single, coherent framework.
As already noted, many other works have already observed the discrepant behavior in R o when compared with R s and R l [19], although these observations have sometimes been made only for larger collision systems [13]. Previous theoretical work has explored the implications of the K T -dependence of the radii [15] but has not specifically considered the role of the (dN ch /dη) 1/3 dependence in the radii. Refs. [16,57] have further emphasized the importance of flow and the space-time structure of the source for understanding the (dN ch /dη) 1/3 dependence of the pp and Pb-Pb radii as modeled by UrQMD, but do not appear to have analyzed the same dependence in detail from the perspective of hydrodynamics.
For the present study, the goal is to explore specifically whether the differences in the (dN ch /dη) 1/3 dependence between the various radii and collision systems can be naturally understood in terms of the space-time picture provided by hydrodynamics. Since the hydrodynamic formalism used here is highly simplified in the interest of clarity, the focus will be placed on obtaining a qualitatively plausible understanding of how hydrodynamics describes the space-time evolution in different collision systems, rather than attempting to quantitatively reproduce the data in detail. This is the subject to which we turn next.

IV. RESULTS
In this work, I have applied the formalism covered in Sec. II to the systems studied in the experimental analyses described in Sec. III. The results are presented in this section and are organized into three areas. In order to understand the multiplicity scaling of the HBT radii, we must first appreciate the ways in which changing the multiplicity in different systems affects the shape of the freeze-out surfaces themselves whose structure the HBT radii are  Pb-Pb collisions. While increasing the multiplicity proportionately increases the size of the system, it also tends to distort the shape of the freeze-out surface. Several additional centrality classes are shown for pp , in order to illustrate the changes in shape which occur at high multiplicity.
supposed to characterize. We begin by examining and comparing these surfaces directly in Fig. 2 for different systems and centrality classes. This will lead us to consider how this space-time structure can be manifested in the HBT radii, and it is at this point that the Gaussian source approximation will prove useful for guiding intuition, namely, by relating the HBT radii directly to space-time variances of the source function. Finally, having an intuitive feeling for how the multiplicity influences HBT on the basis of hydrodynamics, we finally consider the radii themselves which are extracted according to the above formalism.
A. Freeze-out surfaces First, since HBT interferometry reflects the space-time structure of the emitting source in nuclear collisions, it is crucial to examine how the freeze-out surfaces themselves evolve with the system's multiplicity. This is shown in Fig. 2 for the centrality classes under consideration; similar plots were also studied in Ref. [17]. Several additional classes for pp have also been added at large multiplicities for illustrative clarity. One notices immediately a conspicuous difference in the behavior of Pb-Pb collisions as compared with p-Pb and pp collisions, especially at large centralities. In the former, the scaling with multiplicity primarily affects only the enclosed space-time volume of the system, without dramatically altering the shape of the freeze-out surface itself. In small systems, however, the growth of the enclosed volume with multiplicity is less important than changes to the shape of the freeze-out surface, especially at high multiplicity. Remarkably, in extreme pp collisions, the system freezes out in the center first, followed by freeze-out at the edges. Viewed as an animation, one would see such a system as a ring of quark-gluon plasma in the transverse plane, expanding and narrowing until final freeze out at a time τ − τ 0 ∼ 5 fm/c and radius r ∼ 4.5 fm. One should therefore expect significant differences in the scaling of the radii in Pb-Pb collisions when compared to that in pp or p-Pb collisions at similar multiplicities.
One also notices, by comparing the slices along the x and y axes, that pp and p-Pb collisions exhibit greater rotational symmetry than Pb-Pb : this reflects the fact that the increasing importance of event-by-event fluctuations in the location and violence of collisions between subconstituents destroys the rather tight correlation of collision centrality and impact parameter observed in collisions between large nuclei when going to small collision systems [58,59]. It is also worth underscoring that the contours in Fig. 2 correspond to fixed centrality classes, not necessarily fixed multiplicities. The comparison at fixed multiplicity will be shown below.
From these reflections we may already draw a very important preliminary conclusion: hydrodynamics does not in general predict a universal scaling of the R 2 ij with dN ch /dη which is irrespective of the system size. Conversely, even highly simplified hydrodynamic models (like the one considered here) predict a non-universal scaling for high-multiplicity pp collisions. Although this observation is focused on HBT and a particular definition of the multiplicity, it presumably applies to other space-time observables and definitions of the multiplicity as well.

B. The emission function S(x, K)
The fact that the freeze-out structure scales differently in large than small collisions should be reflected in the radii as well. This can be justified quantitatively by considering the scaling of the source variances entering the R 2 ij that are obtained using the Gaussian source approximation. Recall that, using this approximation, the HBT radii can be related directly to the space-time structure of the underlying emission function. The multiplicity dependence of the HBT radii is therefore approximately reflected in the corresponding behavior of the brightest emission regions at a fixed value of K T .
To see this more clearly, we write out explicitly the R 2 ij of interest on the basis of Eq. (5): These relations depend on a total of six source variances: three geometric terms ( x 2 o , x 2 s , x 2 l ), two cross terms ( x ot , x lt ), and a purely temporal term t 2 . Each term probes a different spatiotemporal dimension of the effective emission region which dominates particle production for a given K. Thus, for instance, x 2 o represents the size of a given emission region in the out direction (along the direction of K in the transverse plane).
Similarly, t 2 represents the spread in times over which particles at a given K were typically emitted, while the x ot represents the degree of correlation between the out and time coordinates of particle emission.
Now we wish to see how the scaling of these source variances is reflected in the emission regions shown in Fig. 3, which compares pp , p-Pb , and Pb-Pb collisions at a fixed dN ch /dη = 100 as was done in Ref. [17]. Note that, as was observed in [17], the condition of equal multiplicity requires each freeze-out surface portrayed in Fig. 2 to have the same co-moving volume; visual inspection shows that this condition is very closely satisfied. In addition, however, the color map in Fig. 3 projects the emission function S(x, K) for each system and 0-0.00025%, respectively, with the normalizations determined from Table I. directly onto its respective freeze-out surface, in order to illustrate how the space-time structure of S is influenced by K T and the collision system under consideration. 5 In each panel of Fig. 3, the outlines of the freeze-out surface along the x and y axes are shown as thin, dashed white lines to guide the eye. Bright yellow regions have the highest pion emissivity and will tend to dominate the corresponding source variances as well; dark purple regions contain fluid cells producing the fewest pions and will accordingly have little effect on the source variances. The color scale is normalized between 0 (minimum emissivity) and 1 (maximum emissivity).
Thus, one observes that at small K T , particle emission happens mostly at late times τ ∼ 5 − 8 fm/c within 2 − 4 fm of the origin (r = 0). In Pb-Pb , emission at small K T occurs later and faster than emission at large K T , which originates mainly from the edge of the system over a larger and earlier spread of times. In pp and p-Pb , by contrast, small K T emission happens earlier and more rapidly than large K T emission, as a consequence of the systems' freezing out sooner in the centers than at the edges. Adjusting the K T window thus influences where the emission function is brightest, and thereby provides a tunable filter with which to probe different portions of the freeze-out surface in a controlled way.
The shifting of the highest emissivity regions with K T leads to the well-known K T scaling of the radii [13,14] which can be identified directly in the reduced sizes of the bright yellow regions at large vs. small K T in Fig. 3. This effect is present in all three systems and is a consequence of collective flow: particles emitted from the system's center tend to belong to pairs with relatively small K T values, since the collective motion is comparatively weak there. Particles emitted from the system's edge, on the other hand, are produced by fluid elements which already possess a strong transverse velocity component, and consequently emit particles preferentially with large momenta moving in the same direction. Large-K T 5 Note that, in the case of pp , the normalization of the minimum-bias initial conditions was retuned in order to reach dN ch /dη = 100, since these events are too rare to be conveniently reproduced with the normalization fixed by Table I. This is only a justifiable trick in pp where the multiplicity is not correlated with the impact parameter b the way it is in larger systems, so that the rescaling does not alter the subsequent evolution significantly. In any event, the retuning to dN ch /dη = 100 is used here for purely illustrative purposes, in order to indicate how different systems may possess dramatically different spacetime structures, even at the same multiplicity. emission is thus dominated by the fluid cells at the edge of the system, leading to more compact emission regions, and causing the HBT radii to decrease accordingly [34].
Moreover, the highly elongated 'wing-like' structure of the small systems' freeze-out surfaces is also a result of enhanced collective flow, in which the center of the source freezes out well before the edges do. In this sense, small systems at high multiplicity are quite literally exploding 'rings of fire,' from the perspective of hydrodynamics. This enhanced collective flow in small systems originates from a combination of their reduced sizes (generating larger initial density gradients) and the higher temperatures produced in their interiors [60] which generate a more violent response due to a larger speed of sound (e.g., [61,62]).
Notably, this elongated freeze-out structure leads visibly at large K T to a strong, positive correlation betweenx o (∼ r) andt, implying that x ot > 0 in small systems at high multiplicity. This is opposite to the behavior of large systems at the same multiplicity in hydrodynamics, which clearly tend to havex o andt negatively correlated with one another, implying that x ot < 0 in these systems.
The freeze-out geometry also has implications for the scaling of the other source variances with multiplicity, as shown in Fig. 4 for pp , p-Pb , and Pb-Pb at K T = 450 MeV. Since we work here in the LCMS frame [63,64], in which the longitudinal component of the pair momentum vanishes (K L ≡ 0) for each pair by definition, β L = 0 as well, so that both R 2 s and R 2 l are dominated completely by the system's spatial geometry. The geometric variances ( x 2 o , x 2 s , x 2 l ) shown in Fig. 4 grow approximately monotonically with multiplicity in both large and small systems, reflecting the steady scaling which is already visible in Fig. 2. x 2 s (red up-triangles) is found to be consistently larger than x 2 o (green circles) in all systems, but both scale with dN ch /dη in essentially the same way. x 2 l (orange diamonds) grows more rapidly than x 2 o or x 2 s , reflecting the extended shape of the system in the longitudinal direction.
The behavior of the temporal variances ( x ot , t 2 ) is somewhat more interesting. In Pb-Pb collisions one observes that the emission duration t 2 (purple squares) grows monotonically with multiplicity, while the correlation term x ot < 0 everywhere and decreases monotonically with multiplicity (blue down-triangles). In pp and p-Pb , on the other hand, this monotonic behavior is lost: the emission duration eventually "levels off" and the correlation term actually turns positive x ot > 0 at a critical value of the multiplicity, owing to the wing-like structure shown already in Fig. 3. The panels (a-c) are compared side-by-side in panel (d) on the same scale, in order to make the differences between systems more apparent. Here, K T = 450 MeV. One observes that most of the splitting between systems emerges in the temporal or longitudinal variances ( x 2 l , x ot , t 2 ), whereas less splitting is visible in the transverse geometry ( x 2 o , x 2 s ). All source variances have been averaged azimuthally over Φ K .
In Fig. 4d, we show the same source variances as in panels (a-c), but overlaid on the same set of axes to facilitate direct comparison. We observe a number of critical features.
First, the longitudinal variance scales strongly with the size of the collision system, with noticeable splitting occurring above dN ch /dη > ∼ 12. Similar splitting is seen in x ot and t 2 , for which the freeze-out geometry dictates dramatically different behavior in the different systems, as we have seen previously. Interestingly, the transverse source variances x 2 s and, in particular, x 2 o change surprisingly little between large and small systems when holding the multiplicity fixed. At a superficial level this reflects the observation made in [17] that, if freeze-out occurs at constant density, for fixed multiplicity the co-moving volume must be the same in all collision systems. However, HBT radii are known not to measure the entire freeze-out volume, but only some fraction of it, known as the 'homogeneity volume' [65], whose size is affected by the collective expansion rate at freeze-out. As already discussed and explicitly seen in Figs. 2 and 3, this expansion rate increases from Pb-Pb to pp collisions; Fig. 4d shows that, at fixed multiplicity this increase in radial flow leads to a decrease of x 2 s from Pb-Pb (solid) to p-Pb (dash-dotted) to pp (short-dashed), as anticipated in [17]. Contrary to the expectations in [17], however, this effect is much weaker for x 2 o and can therefore not explain the experimentally observed significantly larger variation of R 2 o than R 2 s when going from Pb-Pb to pp at fixed multiplicity. Instead, as we will discuss next, this last feature can be understood by studying the system size dependence of the other contributions to R 2 o in Eq. (10), caused by the qualitative change in the shape of the freeze-out surface exhibited in Fig. 3.
We are finally in a position to consider the actual multiplicity dependence of the HBT radii in hydrodynamics. This is shown in Fig. 5. As expected, R s and R l follow a nearly universal, approximately linear scaling with (dN ch /dη) 1/3 in all three systems. Since they are dominated by the spatial geometry of the system, their scaling reflects the extensive nature of dN ch /dη, which should be proportional to the system's volume.
R o , on the other hand, is sensitive to both the spatial and temporal sizes of the source, as well as the correlation between the two. As we have just seen, the latter behaves very differently in large and small collision systems at a fixed multiplicity: the wing-like geometry induced by strong collective flow forces a change of sign in the correlation term x ot and a concomitant leveling off of the emission duration t 2 . The combination of these effects is that, in high-multiplicity pp and p-Pb collisions, R o exhibits a much weaker scaling with multiplicity than either R s or R l , whereas in Pb-Pb collisions, the comparatively weaker flow allows R o to grow at a rate similar to that seen in R s and R l . The shallow scaling of R o with (dN ch /dη) 1/3 in pp and p-Pb becomes especially pronounced at higher multiplicities, leading to a slower overall growth with multiplicity. This is how hydrodynamics explains the slope hierarchy observed in the data.
Hydrodynamics may also allow an understanding of the slope non-universality visible across collision systems. Hydrodynamics predicts very similar slopes for R s in all systems, a feature which seems to be fairly well borne out by the data (cf. panel (b) of Fig. 1). For R o , the slope in p-Pb falls squarely in between those of pp and Pb-Pb at large multiplicities, again in surprisingly good agreement with data.
The R l data initially seem to violate the qualitative hydrodynamic tendencies, showing similar slopes between p-Pb and Pb-Pb , but a significantly smaller slope in pp (cf. panel (c) of Fig. 1). On closer inspection, however, the discrepancies may not be as bad as they first appear: at large multiplicities, there is a small but detectable splitting in the model slopes of R l which mirrors that seen in R o and qualitatively, if not quantitatively, resembles the splitting seen in the data. Although a thorough resolution of this tension is beyond the scope of the current work, it is worth speculating as to how the tension originates and how it might be alleviated. The qualitative similarity may originate from the fact that both and might explain the larger longitudinal slope hierarchy seen in the data. Hydrodynamics will therefore tend to underestimate the splitting of R l until it is supplemented with more realistic features, such as a hadronic rescattering phase, which naturally produces a mixing of longitudinal and temporal source properties as is inherent to experimental analyses with finite statistics [64]. Of course, this proposed explanation is highly speculative and should be considered with appropriate caution. For present purposes, it is sufficient to note that the model-to-data discrepancies seen between Figs. 1(c) and 5(c, f) may be plausibly attributed to the simplicity of the model used here, and need not imply any intrinsic limitations of hydrodynamics itself.
There is another respect in which the model used here fails to completely represent the data. It is clear that the multiplicity dependence of R o predicted by hydrodynamics is not linear in small systems and that there is even some slight curvature visible in the Pb-Pb curves shown here. The reason for this is the competition between spatial and temporal information which influences R o : because the freeze-out structure does change with multiplicity (due to changes in the amount of flow) in both large and small systems, one expects these changes in shape to produce deviations from the otherwise linear dependence which would result from a pure rescaling of the system size. Because all hydrodynamic systems change both size and shape with multiplicity (cf. Fig. 2), the scaling of the radii is not in general expected to be perfectly linear.
However, the fact remains that the current hydrodynamic model predicts non-linear (dN ch /dη) 1/3 dependence of the radii which is not obviously reflected in the data, and only reproduces the observed slope hierarchy and non-universality at multiplicities large enough to generate the wing-like structure of Fig. 3. Thus, one might worry that the R o scaling in small systems is too similar to the steeper R s and R l scaling at low multiplicities for the connection with hydrodynamics to be justifiably drawn in this regime. While it remains to be seen whether a more sophisticated model would reproduce the trends observed in the data, it is worth pointing out that nothing in principle prevents hydrodynamics from being applicable even to systems with very low multiplicity [17], and the more crucial question is whether sufficient flow can be generated to weaken the R o scaling also at smaller dN ch /dη once effects like preequilibrium flow [66] have been included in the analysis. This question will have to be answered with a more advanced model than the one employed here.
In the interest of clarity, it is helpful to present the data alongside the model results in a way which isolates the behavior of the slopes in the radii as themselves functions of K T . This can be done by estimating the slope (using simple linear fits) for a given radius in each system separately. Fig. 5 shows the (dN ch /dη) 1/3 dependence of the model results explicitly for only two values of K T , and as already noted, in both cases it is clear that the differences in slopes only begin to emerge above a certain value of the multiplicity which depends on the system in question. For pp and p-Pb shown in Fig. 5, this seems to happen for (roughly) the five largest centrality bins considered, corresponding to dN ch /dη > ∼ 13 in pp , and dN ch /dη > ∼ 42 in p-Pb . Pb-Pb has a nearly constant slope for all multiplicities shown, although R o shows a slight curvature. To isolate the slopes in these high-multiplicity regimes, we fit each radius against (dN ch /dη) 1/3 to a straight line over the five largest centralities and extract the corresponding slope; cf. the illustrative fits included in Fig. 5a. We perform a similar exercise for the experimental datasets presented in Fig. 1, using all published centrality or multiplicity datasets in the pp , p-Pb , and the combined STAR datasets. For each radius and collision system in the model or data, we extract the approximate linear slope as just described and plot it as a function of K T .
The result of this (somewhat heuristic) fitting exercise is plotted in Fig. 6. It needs to be reiterated that this should not be taken as a serious, quantitative comparison between the model results computed in this work and the experimental data. Rather, this is an efficient way of assessing to what extent the high-multiplicity, flow-driven behavior exhibited by small systems is capable, at least in principle, of explaining the non-trivial features observed in the (dN ch /dη) 1/3 -dependence of the HBT radii. The slopes dR ij /dM , with M ≡ (dN ch /dη) 1/3 , are shown for each radius and collision system used in this study. Note that the experimental slopes for A-A have actually been taken from the STAR datasets for Au+Au and Cu+Cu collisions, whereas the model A-A slopes were taken from the Pb-Pb radii plotted in Fig. 5.
Since the slopes of the radii are approximately universal in A-A , this approximation is a reasonable one [13].
The features of non-universality and hierarchy in the slopes can now be clearly seen.
Non-universality is manifested by comparing curves in the same panel, which reveals how the slope changes across collision systems. For instance, panel (a) show that in both the data and the model there is a dramatic splitting of the R o slope in pp , p-A , and A-A , with the magnitude of the slopes ordered by system size (and the pp scaling even turning negative at a certain K T ). By constrast, the slope of R s is approximately universal in all systems, in both the data and the model. Perhaps the largest tension between the model and data on this count occurs in R l , for which the slopes in the pp data are considerably lower than those in p-A and A-A , while the corresponding model calculations show no such dramatic splitting. Even so, the problem could still be more quantitative than qualitative, and a model which relaxed the assumption of boost invariance would likely see an improvement of the agreement with R l .
What I have here called the 'slope hierarchy' is also easily visible in Fig. 6 by comparing same color curves in different panels, i.e., by comparing the behavior of different radii in the same system. Here the verdict is generally encouraging: the R o slope is generically far smaller than that of R s or R l , in both model and data. In general, we also find that the slopes mostly decrease with K T , in both model and data (with a few exceptions) . This is what one should expect to find: larger multiplicities lead to more distorted geometries (cf. Fig. 2) which tend to exacerbate the associated K T scaling which probes close to the edge of the system (cf. Fig. 3).
Finally, it has to be emphasized yet again that the failure of the non-universality and hierarchy of slopes to persist to small dN ch /dη is liable to change with the use of a more realistic model. It is also certainly interesting to consider the possibility that hydrodynamics is valid in pp only at sufficiently high multiplicities, but merges smoothly to some alternative formulation at smaller multiplicities. Nevertheless, regardless of whether hydrodynamics turns out to be valid at low multiplicities, it is still the case that the scaling at sufficiently large multiplicities reproduces much of the observed slope hierarchy and non-universality of the radii in a natural and automatic way. Whether the situation at low multiplicities improves once it is coupled with state-of-the-art hydrodynamic simulations will be investigated in a future study.

V. CONCLUSIONS
In this paper, I have considered the multiplicity dependence of the HBT radii in detail from the perspective of a simplified, hydrodynamic model. The advantage of using such a simplified model is that the most important conclusions are easy to draw. In this case, by comparing the radii measured in large and small systems at a fixed multiplicity, one sees that, according to hydrodynamics, the strength of collective flow in the latter systems is disproportionately stronger than in the former systems, leading to measurable effects on the structure of the freeze-out surface itself.
The fact that R 2 o contains a mixture of space and time information, whereas R 2 s and R 2 l are dominated mainly by spatial geometry, implies that the multiplicity scaling in the out direction in small systems should be substantially weaker than in the side or long directions in systems with strong collective flow. The slope hierarchy observed in the data is therefore an automatic consequence of the hydrodynamic paradigm presented here, suggesting that the evolution of high-multiplicity pp collisions may be driven by a violent, hydrodynamic response to initial density gradients in these systems. More generally, it is clear that the multiplicity dependence of HBT can be used to place non-trivial constraints on any model of collectivity in nuclear collisions, and may provide additional discriminating power for adjudicating between the various mechanisms which have been proposed to explain it. Although a number of features in the calculations outlined above will certainly change quantitatively with additional theoretical improvements, the primary connections between strong collective flow in small systems, the non-trivial evolution of the freeze-out surfaces in large and small systems with multiplicity, and the resulting effects on the HBT radii, are all expected to survive a more rigorous analysis.
Nevertheless, the primary drawback of employing a highly simplified model like the one used here is that one needs to relax a number of strong assumptions and approximations before the model's implications can be taken seriously in a quantitative way. In this vein, there remain several important directions in which this work will be extended, some of which are already underway.
First, a more sophisticated treatment of initial state fluctuations, viscous corrections, and the inclusion of resonance decay effects and a hadronic rescattering phase are features which must be incorporated before a quantitative comparison with data may be legitimately performed. Moreover, as pointed out in [23], systematic uncertainties arising from the parameterization and construction of the correlation function have not been treated in this work and must in general be handled with great care [49,50].
Second, one should consider alternatives to hydrodynamics which might also be capable of reproducing the essential behavior seen in the data. One example is Pythia [67] and its recent extension, Angantyr [68], to the description of collective effects in both small and large collision systems [69]. Alternatively, one could consider to what extent K T scaling and (dN ch /dη) 1/3 scaling of the femtoscopic radii could reflect collectivity which is generated by models with initial-state correlations, such as CGC/IP-Glasma [24,70,71].
In the latter case, initial-state models, which predict some flow already at early times, suggest different kinematic initial conditions in nuclear collisions than those which generate the same flow dynamically throughout the collision history, i.e., in response to initial density gradients. This difference would be undetectable using momentum-space (flow) information alone, but might be identifiable by examining suitable momentum-space observables (e.g., p T , v 2 ) simultaneously with complementary coordinate-space observables (e.g., azimuthally sensitive R 2 ij (K T , Φ K )). For this reason, for a fixed amount of flow and multiplicity, initial-state models will generically predict smaller source radii than hydrodynamic models. Conditioning on both momentum-space and space-time information simultaneously therefore offers constraining power which the use of momentum-space information alone does not.
Third, one should in principle consider the effects of replacing the mid-rapidity charge particle density with forward/backward multiplicity estimators, for which the scaling of the radii could change significantly and might also reveal distinct insights into the systems' dynamics. Doing so quantitatively would of course require relaxing the assumption of boost invariance used here, and would involve solving the full baryon density equation of motion in conjunction with the usual hydrodynamics in 3+1 dimensions.
Fourth and finally, one can also explore whether additional constraining power can be provided by studying higher moments of the distributions of femtoscopic radii which result from the incorporation of event-by-event fluctuations [49,51,71]. By combining these together with other collections of observables, such as the wide-ranging set of radial and anisotropic flow measurements, one could thereby place non-trivial constraints simultaneously on both the space-time and momentum-space evolution of a wide range of high-energy nuclear collisions. A major effort along these lines is deferred to future work.