Do nuclear collisions create a locally equilibrated quark-gluon plasma?

Experimental results on azimuthal correlations in high energy nuclear collisions (nucleus-nucleus, proton-nucleus and proton-proton) seem to be well described by viscous hydrodynamics. It is often argued that this agreement implies either local thermal equilibrium or at least local isotropy. In this note, I present arguments why this is not the case. Neither local near-equilibrium nor near-isotropy are required in order for hydrodynamics to offer a successful and accurate description of experimental results. However, I predict the breakdown of hydrodynamics at momenta of order seven times the temperature, corresponding to a smallest possible QCD liquid drop size of 0.15 fm.


II. INTRODUCTION
If nuclear collision experiments do not probe near-equilibrium matter, then this may have a number of consequences which to my knowledge have not been appreciated before, providing the motivation for this note. Firstly, it would imply that the nuclear experiments do not (directly) probe equilibrium QCD properties as those calculated in firstprinciple lattice QCD calculations. Depending on the degree of non-equilibrium, experiments may be closer or farther away from the QCD phase diagram plane spanned by temperature and baryon chemical potential. While for hydrodynamics, a projection from non-equilibrium space to the equilibrium plane is provided by e.g. the Landau matching condition, for other observables such a projection is not explicitly known. For instance, it is possible that the phenomenon of critical fluctuations associated with the experimental search for a QCD critical point would get modified when experiments probe QCD away from equilibrium.
The understanding that the matter created in high energy nuclear collisions does not need to equilibrate or isotropize locally in order for hydrodynamics to be quantitatively applicable would imply that the "early thermalization puzzle" is to some extent not a genuine puzzle (see more details about this below). On the other hand, having experimental access to a non-equilibrium quantum system could lead to new directions in the field such as e.g. observing nonequilibrium entropy production, properties of non-thermal fixed points or off-equilibrium photo-production.
Finally, if systems created in nuclear collisions do not equilibrate this could naturally explain why proton-proton collision data on azimuthal correlations appears to be so similar to data obtained in nucleus-nucleus collisions. Once small gradients or near-equilibrium is no longer a requirement, hydrodynamics will generically convert initial state geometry and fluctuations into correlations, thus making large and small systems look alike in their azimuthal correlation signals. Pushing this idea even further would imply that any lump of sufficiently high energy density could expand according to the laws of hydrodynamics (with one important caveat which will be discussed below). A natural consequence of this would be the presence of exponentially falling (thermal) spectra as well as potential azimuthal correlations in e + +e − collisions.

III. A HISTORICAL PERSPECTIVE
In recent years, it has been demonstrated that experimental results obtained in relativistic nuclear collisions are well described by hydrodynamic simulations. Based on the paradigm that hydrodynamics requires near-equilibrium in order to be applicable, this successful hydrodynamic description has been interpreted as evidence for a locally equilibrated state of matter (dubbed the quark-gluon plasma) in high energy nuclear collisions.
The question on how the system created in high energy nuclear collisions reaches or at least comes close to equilibrium subsequently has led to a number of developments. In particular, it was realized that because of the expansion of the matter into the vacuum following the nuclear collision, the system would cool and thus freeze into a hadronic gas quickly, at which point the type of correlations observed in experiment could not longer be built up. Thus, it became apparent that a fluid dynamic approximation to the system dynamics had to start early, on a time-scale of τ ∼ 1 fm/c or less [1][2][3][4].
In a seminal article entitled Bottom-up Thermalization [5], Baier, Mueller, Schiff and Son calculated the time when gluons at weak coupling α s ≪ 1 reached thermal equilibrium in a heavy-ion collision, finding [6] τ ≥ 1.5α −13/5 which leads to τ > ∼ 6.9 fm/c for Q s ∼ 1 GeV and α s ∼ 0.3. Clearly, this result seemed to be in tension with the starting time required from the agreement between hydrodynamics and experimental data.
Arnold, Lenaghan, Moore and Yaffe [7] pointed out that a possible way out of the dilemma was that full thermalization was actually not required for a hydrodynamic description, and that local (near-) isotropy of the pressure tensor was sufficient. Thus, the attention of the field shifted towards finding a mechanism to quickly achieve local isotropy (isotropization) rather than full thermalization in high energy nuclear collisions.
One such possible mechanism was that of non-abelian plasma instabilities, specifically the non-abelian Weibel instability, which had been extensively studied by Mrowczynski since the 1980s [8,9]. In a series of numerical studies by a number of groups the growth and saturation of these plasma instabilities was determined for non-expanding systems [10][11][12][13][14], see Refs. [15,16] for a review. While corroborating the initial exponential approach towards isotropy, these numerical studies suggested the system to stall at large pressure anisotropies once the plasma instabilities reached the non-perturbative non-abelian scale and could no longer grow exponentially. Even worse, later studies in expanding systems aiming at more realistically describing experimental nuclear collisions indicated that the effect of plasma instabilities was delayed/diminished to an extent that they could not lead to local pressure isotropy in a time-scale relevant for nuclear collisions at RHIC and the LHC [17][18][19][20].
While full isotropization seemed difficult to achieve within a weak-coupling QCD based framework, it appeared to be reachable much faster in gravitational duals of gauge theories in the limit of infinite coupling. For instance, Chesler and Yaffe [21] report isotropization to occur at τ ∼ 0.7/T for a non-expanding system, roughly translating to τ ≃ 0.35 fm/c when assuming T ∼ 0.4 GeV. However, similar to the case of plasma instabilities, isotropization does take more time when considering the case of expanding systems such as those in nuclear collisions, because expansion constantly tries to drive the system away from local isotropy. This is the reason why in newer studies including expansion [22,23], the isotropization time gets delayed. In particular, it eventually became clear that even for systems at infinite coupling strength the system does not isotropize early. Rather, even at infinite coupling, the pressure anisotropy exceeds 10 percent for all times τ < ∼ 10 fm/c [24]. For completeness, it should be noted that when including inelastic processes in a weak-coupling based descriptions, recent studies [25,26] have demonstrated the approach to isotropy, albeit at times later than those found for infinitely strongly coupled gauge theories. (This better be the case). Thus, the approach to isotropization in expanding gauge theories is now understood both at weak and strong coupling, and indicates long times.
Despite the impressive progress made, I believe it is a correct statement to say that at phenomenologically relevant times of τ ∼ 1 fm/c following a nuclear collision, no theoretical approach (be it weakly coupled or strongly coupled) finds the longitudinal and transverse pressure to agree with each other to better than a factor of two. Obviously, a pressure anisotropy of 50 percent is not close to an isotropic system, let alone a system in thermal equilibrium. By the criterion of Arnold, Lenaghan, Moore and Yaffe, hydrodynamics should not apply.
But it does.

IV. HYDRODYNAMIZATION OR THE ONSET OF HYDRODYNAMIC APPLICABILITY
• Q: How do you people know hydrodynamics applies for pressure anisotropies of 50 percent or larger?
• A: We checked.
Let us consider the following numerical experiment. Take matter described by gauge/gravity duality at infinite coupling or alternatively described by kinetic theory at some finite (constant) value of the coupling. Let the matter be initially at rest in equilibrium with some temperature T i in flat Minkowski space. Then, at a time t ∼ 0, the spacetime suddenly starts to expand in one dimension so that it effectively mimics the effects of so-called Bjorken flow [27]. The symbols in Fig. 1 show the response of the matter (at various values of the coupling λ) in terms of the ratio of longitudinal to transverse pressure as a function of time. The matter is initially in equilibrium so P T = P L (zero pressure anisotropy) and also tends to equilibrium at late times when the expansion becomes very slow. However, for t ≃ 0 when the expansion is most rapid, the matter is clearly not locally isotropic, and deviations from local isotropy become larger as the coupling is decreased.
Also shown in Fig. 1 are results from a hydrodynamic gradient expansion to first and second order in gradients, respectively. One observes that hydrodynamics quantitatively matches the exact results whenever the pressure anisotropy is 50 percent or smaller.
This 'unreasonable success' of hydrodynamics in describing systems with pressure anisotropies of order unity is neither limited to this one example nor to AdS/CFT dynamics, nor exclusively to previous work by the present author, cf. Ref. [23,[28][29][30][31][32].
While of course no general proof, the above numerical experiment indicates that hydrodynamics is able to give quantitatively accurate descriptions even when the matter not locally isotropic. The timescale at which hydrodynamics first is able to closely approximate the subsequent dynamics of the exact underlying microscopic theory has been dubbed hydrodynamization time [33]. At the hydrodynamization time, the matter is typically not locally isotropic. So what sets the timescale for the onset of the applicability of hydrodynamics?

V. HYDRODYNAMIC VERSUS NON-HYDRODYNAMIC MODES
What is hydrodynamics? The equations of hydrodynamics can be derived using a multitude of approaches. Some assume the system to be close to thermal equilibrium, others assume a weakly coupled microscopic particle description (kinetic theory). In my opinion, the most general derivation of hydrodynamics follows the approach of effective field theory (EFT).
According to this viewpoint, hydrodynamics is the EFT of long-lived, long-wavelength excitations consistent with the basic symmetries of the underlying system. The fundamental variables of the EFT are that of a fluid: pressure P, (energy) density ǫ and fluid velocity u a . To lowest (leading) order in the EFT, only algebraic combinations of these quantities will enter the description 1 . Corrections can then be systematically obtained by considering gradients of the fundamental variables.
Applying this EFT approach to the energy-momentum tensor for a relativistic system in three dimensions leads to the well-known expansion where g ab denotes the space-time metric, η is the shear viscosity coefficient and the symbols <> denote a particular symmetric projector that dedicated readers can look up e.g. in Ref. [35]. With Eq. (2), conservation of energy and momentum ∇ a T ab = 0 then are the relativistic Navier-Stokes equations. Many articles have been written about non-causality of the relativistic Navier-Stokes equations; I will simply ignore this issue here because it is somewhat tangential to the following discussion. The above EFT derivation does at no point invoke the presence of an underlying particle-based, kinetic description of the matter. However, given the requirement of the small gradients, it does seem to require the system to be close to isotropy. So what if gradients were not small in a particular situation of interest? Obviously, stopping at first order in a gradient expansion would not be a good approximation. However, one could try to include higher order gradient corrections to obtain a good approximation. I will try to elucidate what happens in this case through a particular example.
For pedagogical purposes, let me pick the example of N = 4 SYM at infinite coupling undergoing Bjorken expansion that has been worked through in a tour-de-force paper by Heller, Janik and Witaszczyk [36]. In this case one has a high degree of symmetry, and the only relevant gradient is ∇ · u = 1 τ . The equations of motion then lead to a solution for the temperature T as a function of τ which may be systematically calculated for small gradients (or equivalently large τ ). (Note that the actual dimensionless expansion parameter is 1 τ T which scales as τ −2/3 in the hydrodynamic limit). Calculating the temperature T (τ ) in a hydrodynamic gradient expansion to order 240 leads to a series of the form [36] T whereτ = τ τ0 , α n constant, and τ 0 setting the initial time (or equivalently temperature) scale. If the gradient expansion was convergent, then we would have succeeded in a (high order) theory of hydrodynamics that was unconditionally applicable also when the gradients are large. Given that for this theory a very large number of coefficients α n had to be calculated, it would be cumbersome if not impossible to generalize this approach to situations with a much lower degree of symmetry (e.g. nuclear collisions), but at least in principle, it would work! Unfortunately, there is mounting evidence that the hydrodynamic gradient expansion generally is not a convergent series. In the cases that have been examined in detail (N = 4 and N = 2 * SYM at infinite coupling, weakly coupled kinetic theory in the relaxation time approximation and Müller-Israel-Stewart (MIS) theory) it was found that α n ∝ n! for large n, thus making the gradient expansion a divergent series [36][37][38][39][40].
However, not all is lost. It turns out that when inspecting the analytic structure of gradient expansions at high orders, it is possible to use a generalized Borel resummation to rewrite the above series for the case of N = 4 SYM as where T hydro (τ ) is well approximated by the first few orders in (3) as long asτ is not too small. In the above expression, ω Borel ≃ ±3.1193 − 2.7471i, and both γ, ω 1 have been calculated in Ref. [36].
There are two things to note about the resummed result (4). First, the exponential multiplying the coefficient γ in (4) can not be recast in terms of the hydrodynamic gradient expansion. It is a truly non-hydrodynamic mode, and its presence explains why the naive hydrodynamic gradient series is divergent.
Second, the numerical value ofω Borel is not an arbitrary number. It happens to be consistent with the first non-hydrodynamic quasi-normal mode frequency of a near-equilibrium 5d Schwarzschild-AdS black holê calculated by Starinets in Ref. [41]. A quasi-normal mode corresponds to a pole of a two-point function in the complex frequency plane located at ω = 2πTω QNM . Linear response then implies the presence of a contribution of the form e iωτ = e i2πωT τ to the one-point function which upon using the leading order expression in (3) for T (τ ) then leads to a result e ∼iωτ 2/3 consistent with (4), cf. [42]. While only the first non-hydrodynamic quasi-normal modeω (1) QNM has been identified in the Borel transform of the gradient expansion, it is likely that all higher non-hydrodynamic modes also will contribute likewise to T (τ ), which has been anticipated through the ellipses in (4). In fact, Buchel, Heller and Noronha were able to show that for the case of N = 2 * SYM the first 10 quasi-normal modes could be obtained from the relevant Borel transform [38].
Thus, the following picture emerges: a naive hydrodynamic gradient expansion of the energy-momentum tensor is divergent because of the presence of other, non-hydrodynamic degrees of freedom. However, the contribution from these non-hydrodynamic modes may be either resummed via a generalized Borel transform, or anticipated through explicitly calculating the non-hydrodynamic mode structure of the energy-momentum tensor for the theory under consideration: The term T hydro (τ ) in (4) is a "generalized" hydrodynamic piece has been dubbed "hydrodynamic attractor" or "all order hydrodynamics" by various authors [37,43]. As remarked above, it is well approximated by a low-order gradient expansion approximation even in regime when gradients are moderately strong (see e.g. Ref. [37]). Thus, even though different in principle, it is entirely conceivable that for many applications involving gradient terms of order unity, T ab hydro will quantitatively be well approximated by the naive, low-order gradient expansion (2) because the divergence of the gradient approximation only becomes apparent when including higher orders. High-temperature perturbation theory in QCD exhibits a similar feature where the leading order (g 2 ) correction to the free energy offers a quantitatively reasonable description of lattice QCD data even for g ≃ 1 [44], while the inclusion of higher order corrections clearly exhibits the divergent series nature of a naive perturbative expansion.
Contrary to the hydrodynamic contribution, the non-hydrodynamic contribution T ab non−hydro will in general not have a universal form, but rather be dependent on the particular underlying microscopic description under consideration ("microscopic" in the sense of QCD, not in the sense of quasi-particles) . This is most easily elucidated when considering the small amplitude (but arbitrary gradient) linear response of the energy-momentum tensor to an initial source S cd (x), where G ab,cd R (ω, k) is the retarded two-point function of the energy-momentum tensor (cf. [45]). For instance, in the case of N = 4 SYM at infinite coupling, G 00,00 (ω, k) possesses only poles in the complex plane. Two of these poles may be uniquely identified as hydrodynamic sound poles, located at ω h = ±c s |k| − 2iη|k| 2 3s when |k| ≪ 1. In addition to the hydrodynamic poles, there is an infinite number of (pairs) of non-hydrodynamic quasinormal modes located at ω = ω nh , . . . [47]. Performing the frequency integration in (7) will pick up contributions at all of these poles, immediately leading to nh t a n (k) = δT 00 hydro + δT 00 non−hydro , where the coefficient functions a h (k), a n (k) depend on the residues from the integration as well as the source function S(k). The integral in (8) cleanly separates into a hydrodynamic piece governed by the sound mode dispersion ω h (k) and an infinite sum over non-hydrodynamic contributions with dispersion ω (n) nh (k). More important than the realization that the energy-momentum tensor can be split into a hydrodynamic and a nonhydrodynamic piece is the fact that the hydrodynamic poles cease to exist at some value of k when the coupling is finite [48]. This implies that for |k| larger than some critical value of |k| c (dependent on the coupling), the hydrodynamic (ω, k) correlator in kinetic theory (from Ref. [46]). At low k/T , there are two hydrodynamic poles and a logarithmic branch cut from −k − i τ R to k − i τ R . Increasing k/T , the hydrodynamic poles acquire larger and larger negative imaginary part until at kc/T ≃ 4.53 they merge with the branch cut and disappear from the physical Riemann sheet. For k > kc, only the branch cut remains. Grey lines are projections of trajectories to the k = 0 plane. component vanishes from the spectrum and is replaced by purely non-hydrodynamic behavior. At least for |k| > |k| c , hydrodynamics has broken down.
One may criticize that N = 4 SYM is a very special microscopic theory, and worry about drawing general conclusions based exclusively on N = 4 SYM. However, it turns out that when calculating G ab,cd R (ω, k) in kinetic theory in the relaxation time approximation [46], similar conclusions apply. In kinetic theory, G 00,00 R (ω, k) generally has two hydrodynamic poles which are located at ω h = ±c s |k|− 2iη|k| 2 3s when |k| ≪ 1. In addition to these hydrodynamic poles, G 00,00 R (ω, k) exhibits a logarithmic branch cut which may be taken to run from −|k| − i τR to |k| − i τR where τ R = 5 η sT is the relaxation time in kinetic theory. It is interesting to note that when increasing |k| beyond some critical value |k| c , the hydrodynamic poles pass through the logarithmic cut onto the next Riemann sheet, and effectively cease to exist (see Fig. 2). Only the non-hydrodynamic branch cut remains, implying that for |k| > |k| c , hydrodynamics has broken down. I will summarize the above observations in the form of a Lemma. Given the existence of a local rest frame, hydrodynamics offers a valid and quantitatively reliable description of the energy-momentum tensor even in non-equilibrium situations as long as the contribution from all nonhydrodynamic modes can be neglected.
Proof: Consider matter possessing a local rest frame everywhere in the space-time patch of interest, such that the local energy density is non-negative in any frame. Pick a convenient global frame ("laboratory frame") and consider the Fourier decomposition of the energy-momentum tensor in this frame. Now consider real time perturbations δT ab (t, k) around the Fourier zero mode T ab background in the laboratory frame. At some initial time t 0 , the difference between the local energy-momentum tensor and the background can be viewed as an initial perturbation S ab (t 0 , k). In the limit of small perturbation amplitude |S ab | → 0, linear response theory applies, cf. Eq. (7). Furthermore, the retarded two-point function G R is known to be given by Navier-Stokes hydrodynamics [45] in the small wavenumber limit k → 0. The two-point correlator in Navier-Stokes hydrodynamics possesses hydrodynamic poles (shear and sound poles) in the complex frequency plane. Contour integration as in Eq. (7) will pick up these poles and lead to a hydrodynamic contribution to δT ab . As k is increased, the location of the hydrodynamic poles may shift, and they may even disappear from the spectrum completely at some critical wave-number. In addition to the hydrodynamic poles, new, non-hydrodynamic singularities may appear in the complex frequency plane. These nonhydrodynamic singularities, upon contour integration in Eq. (7) will lead to a non-hydrodynamic contribution that has to be added to the hydrodynamic part of δT ab as in the example given in Eq. (8). As the amplitude S ab is increased, other, non-linear structures will contribute to δT ab which can be expressed as a sum over integrals of n-point functions with the appropriate powers of the source S ab . In the limit of small wave-number, these non-linear corrections to the hydrodynamic part will, upon resummation, shift the hydrodynamic poles locally, and in addition contribute new structures familiar from the hydrodynamic gradient expansion, cf. Eq. (2). Away from k = 0 the non-linear hydrodynamic part of δT ab is analytically connected to this familiar structure, giving rise to the generalized hydrodynamic attractor form, until some critical wave-number k = k c is reached. In addition to the hydrodynamic part, there will be non-linear corrections to the non-hydrodynamic contribution of δT ab . Thus, the global energymomentum tensor can be written as T ab = T ab background + δT ab = T ab hydro + T ab non−hydro . If the non-hydrodynamic contribution can be neglected, the lemma follows trivially.
(Mathematicians probably would want to see a more formal proof than this, so the above lemma should probably be called a conjecture).
The above lemma may seem trivial: once non-hydrodynamic modes are absent, how can the energy-momentum tensor be described by anything else but hydrodynamics? However, when phrased in this fashion, hydrodynamics neither requires equilibrium, nor isotropy, nor infinitesimally small gradients. (It does require the presence of a local rest frame, though [34]). The applicability of hydrodynamics is exclusively determined by the relative importance of non-hydrodynamic modes.
As it stands, the above lemma has at least one potentially important consequence, which is phrased as a Dilemma. The phenomenological success of hydrodynamics in describing experimental data from high energy nuclear collisions does not imply near-equilibrium behavior of the matter. Experiments do not directly probe the equilibrium QCD phase diagram at finite T, µ, but explore trajectories in a space with at least one more (non-equilibrium) direction.
Incomplete equilibration in nuclear collisions is a fact well known to all heavy-ion hydro practitioners in the field and has been pointed out more than a decade ago by Bhalerao, Blaizot, Borghini and Ollitrault in Ref. [49]. The above lemma allows me to go one step further since not even near-equilibrium is required for hydrodynamics. The second part of the dilemma is a direct consequence of the first, yet it probably is not as widely appreciated. I have tried to visualize the last point of the dilemma in Fig. 3. To generate the hypothetical trajectories I have matched the pressure  [46]) and strong coupling ('gauge-gravity', [48]) frameworks (the curve labeled 'approx' is a fit to the exact numerical values from Ref. [48], see text). The regime of applicability for kinetic theory is λ ≪ 1, while for gauge-gravity duality λ ≫ 1 is required. Note that despite the difference in weak-coupling and strong-coupling frameworks, the resulting hydrodynamic breakdown scale is quantitatively similar for moderately strong couplings λ ≃ 10 − 20.
anisotropy P L /P T to the momentum anisotropy parameter ξ, first defined in Ref. [50]. I use ξ ∈ [0, ∞) to express the degree of non-equilibrium ("non-equilibriumness"), where ξ = 0 corresponds to the case of local equilibrium. There is an extensive literature in anisotropic hydrodynamics which makes the connection between P L /P T and ξ precise (see e.g. the instructive lecture notes by Strickland, Ref. [51], Eq. (3.19)). For the pressure anisotropy itself, I used Navier-Stokes hydrodynamics in Bjorken expansion (cf. [24]) and a viscosity that is given by η s = 1 4π for T > 0.17 GeV and rises linearly as the temperature is lowered to reach η s = 1 at T = 0.1 GeV (cf. [52,53]). The temperature and chemical potential dependence result from a drastically simplified version of Hung and Shuryak, Ref. [54], with just one massive degree of freedom in the hadron gas phase, and assuming constant baryon density to entropy ratio. Trajectories are started at τ = 1 fm/c with multiplicities representative of experimental measurements [55] converted to temperature values T (τ = 1fm/c) as in Ref. [56]. Clearly, none of these choices do justice to the much more accurate descriptions that are currently available. However, I expect the sketch in Fig. 3 to be qualitatively reliable.

VI. BREAKDOWN OF HYDRODYNAMICS AND TESTS
According to the central lemma in the previous section, hydrodynamics can be used to describe a system if a local rest frame exists and non-hydrodynamic modes are sub-dominant. I have nothing new to say about how to test for the presence of a local rest frame, so I will simply follow everyone else's approach and assume that a local rest frame exists. (This is very likely wrong [34] and actually should be thought about more, but I leave this task to dedicated readers).
Contrary to the existence of local rest frames, quite a bit of knowledge now exists on those non-hydrodynamic modes. For the case of kinetic theory with relaxation time τ R , we know the analytic structure of non-hydrodynamic modes in the two-point function of the energy-momentum tensor, and we know that hydrodynamic modes completely vanish from the spectrum at |k| c τ R ≃ 4.5 [46]. Using τ R = 5η sT and η s ≃ 5 g 4 for N = 3 color QCD [57], one finds |k| c ≃ 0.02λ 2 T in terms of the 't Hooft coupling λ ≡ g 2 N . From gauge gravity duals at large but finite λ, non-hydrodynamic modes imply the breakdown of hydrodynamics at |k| c ≃ πT 4 ln 6.65λ 3/2 1+105λ −3/2 (see Fig. 4 for exact numerical results from Ref. [48]). These kinetic theory and gauge-gravity results for the hydrodynamic breakdown scale k c are compared in Fig. 4. It is curious to note that -despite the difference of the kinetic theory and gauge-gravity duality approaches -  [59]). No hydrodynamic curves are shown, but it is known that hydrodynamics well describes the experimental data in the regime indicated as 'hydro' in the plot [60], possibly extending up to pT ≃ 3 GeV. By contrast, for pT > ∼ 4 GeV, the experimental data seems to deviate systematically from the low-momentum behavior, and I have labeled this region 'non-hydro'.
the results for k c turns out to be quantitatively similar in both approaches for moderate values of λ ≃ 10 − 20.
Naively applying the results for k c to QCD with α s ≡ g 2 4π ≃ 0.5 leads to kc(λ≃19) T = 4 − 7, where the higher value is actually coming from the kinetic theory result. Let's be optimistic and take |k| c ≃ 7T . Thus, I claim that no hydrodynamic description is possible for QCD systems smaller than k −1 c ≃ 0.15 fm at a typical QCD-scale temperature of T ≃ 200 MeV. The prediction that hydrodynamics must break down for |k| −1 < 0.15 fm is most likely somewhat useless, because it is very hard to falsify experimentally in nuclear collisions given that the mean proton radius is 0.86 fm. However, the limit of 0.15 fm at least constitutes an actual numerical conjecture for the lower bound of the smallest possible droplet of QCD liquid.
It should be pointed out that the results |k| c ≃ 4.5 τR and the numerical gauge-gravity results shown in Fig. 4 are quantitative upper bounds on the domain of hydrodynamic applicability in weak and strong coupling scenarios. However, it is likely that hydrodynamics breaks down before reaching these values of |k|. One reason why hydrodynamics may break down earlier would be that while hydrodynamic modes do not vanish for |k| < |k| c , they may become subdominant to certain non-hydrodynamic modes, namely those which happen to be closer to the origin of the complex ω plane. For N = 4 SYM at λ → ∞, this seems to happen at around |k| ≃ 2πT , which is at around the same value as k c obtained at finite λ. Yet another reason why hydrodynamics may break down at scales below |k| c could be non-linear effects. For a particular initial condition, numerical studies of N = 4 at λ → ∞ including full non-linear effects by Chesler [58] seem to indicate a breakdown of hydrodynamics at a scale |k| ≃ T . These results are fully consistent with the "most optimistic" result k c ≃ 7T and the resulting hard upper bound for the hydrodynamic breakdown scale, but clearly leave room for sharpening the prediction for |k| c .
A much more direct route to experimentally constrain k c could be provided by high-momentum data on flow coefficients, see Fig. 5. Experimental data for collective flow harmonics in Pb+Pb collisions suggests a change in behavior in the regime between p T = 3 GeV to p T = 4 GeV. The low momentum region is well described by hydrodynamics [60]. Assuming that measured particles originated from a constant-temperature freeze-out surface at T = 0.17 GeV, this would indicate a breakdown of hydrodynamic behavior at pT T = 18 − 23. In order to relate this scale to the hydrodynamic breakdown scale k c < ∼ 7T , quantitative calculations of the location of non-hydrodynamic modes in an expanding system are needed.
In the hadronic phase, kinetic theory would predict hydrodynamic modes to dominate for k < k c ∝ 1 τR , while non-hydrodynamic (particle) modes dominate for k > k c . As the temperature is lowered, τ R ∝ η sT increases strongly [52] until k c falls below the typical system wave-number. From this point onward, most of the system dynamics proceeds according to the non-hydrodynamic particle kinetics, providing a qualitative understanding of the transition from hydrodynamic to particle cascade dynamics ("freeze-out").
The above statements involve hard lower bounds on the smallest scales at which hydrodynamics applies, and a qualitative understanding of the freeze-out transition. However, a quantitative test of the applicability of hydrodynamics, e.g. through testing sensitivity of results with respect to non-hydrodynamic modes, is desirable in the case of nuclear collisions . Fortunately, the workhorse of relativistic viscous hydrodynamics simulations, "causal relativistic viscous hydrodynamics" (which goes by many names and acronyms but is usually associated with the work of Müller, Israel and Stewart [61,62]) does contain a non-hydrodynamic mode buried within, which may be exploited for testing  purposes. Specifically, besides the usual hydrodynamic modes, the energy-momentum tensor two point function contains a pole located at ω nh = − i τπ , where τ π is the "viscous relaxation time" that also controls the size of the second order gradient term in the one-point function of T ab . Any current numerical hydrodynamics simulation of the matter produced after a relativistic nuclear collision needs a specific value for τ π as an input. Simulators choose values of τ π as they see fit, given that the "correct" value for τ π for QCD is not known, and that primary interest is in extracting information about η s , ζ s , not some obscure second-order transport coefficient. However, varying τ π around some "fiducial" value does vary the decay-time of the non-hydrodynamic mode inherent to causal relativistic hydrodynamics, thus offering a direct handle on the sensitivity of final results on the nonhydrodynamic mode. This can be implemented in practice in relativistic viscous hydrodynamic simulations by running simulations at multiple values of τ π and expressing final results in terms of a mean value and a systematic error bar covering the variations of final results from changing τ π . Examples are shown in Fig. 6 for the case of central p+Au collisions at various values of √ s and p+p collisions at √ s = 7 TeV. While the sensitivity on non-hydrodynamic modes is not vanishingly small, the error bars do seem to signal the applicability of hydrodynamics to both p+Au and p+p collisions in general. However, in the case of p+p collisions and dN dY < 2, the error bars become large, signaling strong sensitivity of the result to non-hydrodynamic modes. This empirical result seems to indicate that hydrodynamics breaks down in p+p collisions for dN dY < 2. This multiplicity value of the hydrodynamic breakdown corresponds well to the results derived by Spalinski [63]. It would be interesting to repeat these sensitivity tests for hydrodynamics with a different non-hydrodynamic mode structure, for instance along the lines suggested in Ref. [64].

VII. CONCLUDING REMARKS
1. I would argue that there is hard experimental evidence, e.g. through the phenomenon of jet modifications, for the presence of strongly interacting QCD matter created in nuclear collisions. As argued in this note, I am doubtful about the hard evidence for this matter to be equilibrated.
2. The physics of non-hydrodynamic modes is a rich and a barely studied subject. Given that non-hydrodynamic modes play an important role in the applicability and breakdown of a hydrodynamic descriptions, I believe those non-hydro modes should receive more attention, from theorists and experimentalists alike.
3. The central lemma in section V also applies to the case of diffusion, not only momentum transport. In particular, this implies that a constitutive equation of the form J = σE could hold in the early-time, out-of-equilibrium regime following a nuclear collision, if non-hydro modes are sub-dominant. This could potentially explain a longer-than-expected life-time of the magnetic field which is critical to experimental detection of the Chiral Magnetic Effect [67] (see also Ref. [68]). 4. As outlined in the central dilemma in section V, the experimental search for the QCD critical point will necessarily explore trajectories in some non-equilibrium space (cf. Fig. 3). This implies that the standard equilibrium theory of critical fluctuations strictly speaking does not apply, and one should try to understand non-equilibrium effects (see e.g. Ref. [69]) in order to correctly interpret the experimental data.
5. In view of the 'QGP drop size lower bound' of 0.15 fm, it is maybe not surprising that the matter created in p+p collisions would behave hydrodynamically. At this scale, however, p+p collisions may not be the ultimate drop size test. QCD-QED couplings allow fluctuations of electrons to e.g. quark pairs, thus opening up the possibility of local energy deposition reminiscent of p+p collisions occurring in e + -e − collisions (cf. Refs. [70][71][72]). Data on e + -e − collisions taken at e.g. LEP should be re-analyzed with modern tools in order to find (or rule out) hydrodynamic behavior in these systems. 6. The fact that experimental data shows a qualitative change in trend from hydrodynamic behavior at low momenta to non-hydrodynamic behavior at high momenta suggests a potential experimental handle on the hydrodynamic breakdown scale k c in QCD. This potential connection should be made quantitative in further studies.
7. The entire discussion in this note ignores the presence of hydrodynamic thermal fluctuations, which arise in SU (N ) gauge theories at any finite number N . The subfield of relativistic hydrodynamics with thermal fluctuations is still in its infancy, but potentially can have important phenomenological consequences [45,[73][74][75][76][77][78][79][80]. 8. A recurring problem of non-standard cosmology (so-called "viscous cosmology", cf. [81,82]) seems to be that "interesting" deviations from standard cosmology occur when gradient corrections become order unity. In the "old-fashioned" picture of hydrodynamics, this was not acceptable since order unity corrections heralded the breakdown of applicability of the theory. In view of the central lemma in section V, it could be interesting to determine the relevant non-hydrodynamic modes in cosmology and re-evaluate the regime of applicability of viscous cosmologies.