Azimuthal Anisotropies at High Momentum from Purely Non-Hydrodynamic Transport

In the limit of short mean free path, relativistic kinetic theory gives rise to hydrodynamics through a systematically improvable gradient expansion. In the present work, a systematically improvable expansion in the opposite limit of large mean free path is considered, describing the dynamics of particles which are almost, but not quite, non-interacting. This non-hydrodynamic"eremitic"expansion does not break down for large gradients, and may be useful in situations where a hydrodynamic treatment is not applicable. As applications, azimuthal anisotropies at high transverse momentum in Pb+Pb and p+Pb collisions at $\sqrt{s}=5.02$ TeV are calculated from the first order eremitic expansion of kinetic theory in the relaxation time approximation.


I. INTRODUCTION
Is there a simple description for transport when hydrodynamics fails?
At low momenta, relativistic hydrodynamics has been tremendously successful in offering quantitative descriptions and predictions of experimental data from high energy nuclear collisions (see Ref. [1][2][3][4] for recent reviews.) Hydrodynamics breaks down when nonhydrodynamic modes start to dominate over hydrodynamic modes, and it has been suggested that the experimentally observed peak in azimuthal anisotropies at transverse momenta p T 4 GeV indicates the transition from hydrodynamic to non-hydrodynamic transport [5]. For conformal field theories at weak (strong) coupling, this transition happens at a well-defined momentum scale k c [6,7]. For momenta above k c , the lifetimes of hydrodynamic modes are shorter than those from non-hydrodynamic modes, hence at late times bulk transport will be dominated by purely non-hydrodynamic degrees of freedom. While many studies exist that include both hydrodynamic and non-hydrodynamic modes, little is known about the phenomenological implications and observational consequences of pure non-hydrodynamic transport where all hydrodynamic mode contributions have been turned off. The present work is meant as a step in this direction by studying non-hydrodynamic transport for the case of relativistic kinetic theory.
For vanishing mean free path, kinetic theory corresponds ideal (non-viscous) fluid dynamics that is described by the Euler equation [22]. Small, but non-vanishing mean-free path corrections give rise to viscous fluid dynamics, described by the equations of Navier and Stokes [23,24]. Higher order corrections to the small mean free path regime can be systematically calculated [25]. The opposite limit of large mean free path is known as rarefied gas dynamics or high Knudsen number regime [26][27][28], and in the extreme case of infinite mean free path gives rise to non-interacting (or free-streaming) particle dynamics. For infinite mean free path, the classic kinetic equations can be solved analytically using the method of characteristics, leading to ballistic evolution. In a sense, ballistic evolution and ideal fluid dynamics are analogues of each other, corresponding to opposite extreme limits of infinite and zero mean free path, respectively. However, while the systematic small mean free path expansion has been recognized to lead to viscous fluid dynamics, the equivalent systematic expansion at large but finite mean free path seems to have received less attention in the high energy physics literature, except for two works [29,30]. The present work is meant to consider the phenomenological consequences arising from such a systematic expansion, extending in particular the pioneering work by Borghini and Gombeaud [30] to the case of large momenta. For large mean free path, particles rarely interact, similar to hermit crabs in their natural environment, hence this systematic expansion will be referred to as "eremitic" expansion in the following.

II. SETUP
I will consider a system of massless on-shell classical particles with a continuum distribution of locations x and momenta p at any given time t. Because the particles are massless, their dynamics will be governed by relativistic kinetic theory, although it should be straightforward to modify the discussion for massive particles with non-relativistic dynamics. The relativistic Boltzmann equation is given by where f (t, x, p) is the on-shell phase-space particle distribution function, p µ = (p 0 , p) is the particle's four momentum, and C[f ] is the collision kernel which has the property that it vanishes both in equilibrium as well as for non-interacting particles. The collision kernel depends on the details of the particle interactions, and is usually a (complicated) functional of the particle distribution function f . In order to give a more hands-on treatment, it will be useful to consider a concrete and simple example for the collision kernel, such as the relaxation time (or BGK [31]) approximation where In this equation, τ R is the relaxation time (proportional to the mean free path), u µ = (u 0 , u) is a collective four-velocity vector and f eq [f ] is the pseudo-equilibrium distribution function for a configuration given by f . For this work, the mostly plus metric convention g µν = diag(−, +, +, +) will be used such that p µ u µ = −p 0 u 0 + p · u.
If the system was in equilibrium with a temperature T in a local rest frame given by the four vector u µ in some global coordinate system, then the equilibrium distribution function for classical particles can be taken as Out of local equilibrium, there strictly speaking is no temperature, but for classical particles, a local rest frame and an energy density can always [32] be obtained from the local energymomentum tensor [33,34] T µν (t, where θ(x) denotes a step-function. From the energy-momentum tensor, the local energy density (t, x) and local rest-frame four vector u µ (t, x) can be obtained as the time-like eigenvector and associated eigenvalue, together with the normalization constraint u µ u µ = −1. For massless particles in local thermodynamic equilibrium where f = f eq , the energy density obtained as the time-like eigenvector of T µν is related to the temperature T as Out of equilibrium, I define a pseudo-temperature (also denoted by T ) from the energy density as T = 1/4 (cf. the discussion in Ref. [2]). This pseudo-temperature, together with the time-like eigenvector u µ of T µν , are used to define the pseudo-equilibrium distribution , where the functional dependence on the non-equilibrium particle distribution f has been denoted explicitly.
A. Review of small mean free path expansion For small mean free path, the system is close to equilibrium and the collision kernel is almost vanishing, C[f eq ] 0. In the relaxation time approximation, the mean free path in kinetic theory is proportional to τ R , cf. Ref. [2]. To set the stage, consider first the well-known hydrodynamic gradient expansion of Eq. (1) for small mean free path τ R T 1.
To simplify the analytic treatment, I will consider the case of a conformal system where τ R T = const, but one can expect results to generalize to non-conformal cases. To leading which gives rise to an energy momentum tensor of the form where P = 3 is the local equilibrium pressure for a gas of massless particles. Conservation of this energy-momentum tensor ∂ µ T µν = 0 can be recognized as the relativistic equation of continuity and the Euler equation, respectively. Thus the zeroth order expansion in small mean free path corresponds to zeroth order, or ideal, fluid dynamics.
To first order in the small mean free path expansion of Eq. (1) one has f f fluid,(0) + f fluid, (1) , with [2] f fluid, in the relaxation time approximation for C where f eq denotes the derivative of the equilibrium distribution function with respect to p µ u µ /T and σ µν = ∂ µ u ν + ∂ ν u µ − 2 3 (g µν + u µ u ν )∂ λ u λ is the shear-stress tensor. Calculating the energy-momentum tensor for f = f fluid,(0) + f fluid, (1) one finds Here η is the shear viscosity coefficient that for a conformal system is usually expressed as a ratio with respect to the pseudo-equilibrium entropy density s = +P T . For the kinetic theory at hand, one finds η s = τ R T 5 [2]. With the energy-momentum tensor given by Eq. (11), the conservation equations ∂ µ T µν = 0 can be recognized as the relativistic Navier-Stokes equations. Therefore, the first order expansion in small mean free path corresponds to first order, or viscous, fluid dynamics. Higher expansion orders may be systematically generated using this procedure.
Note that in the small mean free path expansion, the relevant expansion parameter is τ R times a typical gradient strength, cf. Eq. (10). This implies that the expansion fails for large gradients (see the discussion in Ref. [6]).
This equation can be solved analytically using the method of characteristics: where f init (x, p) is the particle distribution function at initial time t = 0. The energymomentum tensor for the zeroth order eremitic expansion is given by (4), which requires specification of f init . Some example cases will be considered below.

C. First order eremitic expansion
The first order eremitic expansion for the distribution function is given by with in the relaxation time approximation where the free-streaming result (12) has been used.
The defining equation for f hermit,(1) is similar to (12), but with a non-vanishing constant on the rhs. Using again the method of characteristics, one finds for the first order eremitic Higher order eremitic corrections can be obtained systematically by repeating this procedure. Note that in the large mean free path expansion, the relevant expansion parameter is τ −1 R times an integral over the characteristic, cf. Eq. (16). This implies that the eremitic expansion fails for late times/large distances, which corresponds to the case of small gradients. This is because the expansion is around non-interacting particles, and corrections from interactions build up coherently along the characteristic for small gradients.

D. Collective Modes
Because hydrodynamic and eremitic expansions are opposite limits of kinetic theory (1), their collective mode structure can also be expected to be different. The collective modes of Eq. (1) in the relaxation time approximation have been analyzed in Ref. [6] for constant τ R T (see Ref. [18] for the case of momentum dependent relaxation time), and results are summarized here in order to keep this work self-contained.
The collective modes can be calculated as the singularities of the retarded two-point function G R (ω, k) of the energy-momentum tensor, with ω, k the conjugate Fourier momenta to t, x (see e.g. Ref. [2]). For simplicity, I will only discuss the singularities of G R in sound channel. For fluid dynamics, the singularities of G R (ω, k) in the sound channel are simple poles located at These are the familiar sound modes.
By contrast, the collective modes for the eremitic expansion can directly be obtained by a Fourier transform of p µ ∂ µ f hermit,(1) using the setup in Ref. [6]. One finds that the collective modes in the eremitic expansion are logarithmic branch cuts emanating from the branch Not surprisingly, the analysis of the collective modes contained in (1) for general τ R contains both of these types of singularities, logarithmic cuts emanating from branch points and hydrodynamic poles located at ω ± hydro (k) 1 . The hydrodynamic poles approach the fluid dynamic results (17) in the limit of small τ R |k|, as they should. The branch points approach the eremitic results (18) in the limit of large τ R |k|, as they should. The situation is summarized in Fig. 1, which depicts the singularity structure of the sound channel two-point function in the complex frequency plane.
It is common to refer to the collective modes corresponding to the sound poles ω hydro as hydrodynamic modes, and label the modes corresponding to the branch cuts ω cut as nonhydrodynamic modes. Thus, the fluid dynamic expansion of kinetic theory contains only hydrodynamic modes, the eremitic expansion contains only non-hydrodynamic modes, while kinetic theory without any expansion contains both. 1 For a common choice of the location of the logarithmic branch cut, the hydrodynamic poles move through the cut onto the next Riemann sheet for kτ R > ∼ 4.5313912 · · · ≡ k c τ R and thus are no longer present on the fundamental Riemann sheet [6]. However, this behavior is not generic because other choices of the logarithmic cut location may be employed [18]. What is generic, however, is that for k > k c the hydrodynamic poles are farther away from the real axis then the non-hydrodynamic branch cut, implying late-time transport to be dominated by non-hydrodynamic degrees of freedom.

III. ANALYTIC EXAMPLES
While powerful methods exist to solve the Boltzmann equation (1)  functions that can be written as with F and Λ arbitrary functions. Within this class of examples, the zeroth order eremitic expansion leads to and the associated energy momentum tensor is Examples of Λ, F where T µν hermit,(0) can be calculated analytically will be discussed in the following.

A. Single Gaussian Hot-Spot
The first example is that of a single Gaussian hot-spot, with an initial pseudo-temperature distribution given by with σ, x 0 , T init controlling the width, location, and height of the hot-spot, respectively. In addition, let me consider the case where the initial particle distribution function is in local equilibrium, For this case, the zeroth order eremitic expansion leads to with energy momentum tensor components given by The pseudo-temperature T and flow vector u µ corresponding to this energy momentum tensor can be obtained via eigenvalue decomposition, cf. Eq. (5). One which define f eq [f hermit,(0) ] in Eq. (16). For further use I will introduce It is useful to consider the early time regime t σ of this case, which can be treated analytically. For small t, the eigenvalue decomposition of T µν hermit,(0) leads to and hence one finds However, one notes that to this order in t σ, the space average d 3 xf hermit,(1) vanishes.
The space average of the first order eremitic expansion is given by where the integration variable has been shifted. With the results for T, u µ from (27) one where u 0 may be obtained from Eq. (27) with u µ u µ = −1. For early times t σ one finds The mean square particle momentum, defined as remains constant as a function of time in the absence of interactions, as can readily be observed from the free-streaming result Interactions can be expected to lead to an increase of p 2 as a function of time, such that a potential experimental observation of p 2 carries indirect information about the interaction strength. The time-dependence of p 2 in first-order eremitic expansion can be found as which increases as a function of time as expected. Fig. 2 shows a comparison of this analytic early time approximation to the first order result p 2 [f hermit,(0) + f hermit, (1) ] obtained via numerical integration 2 of Eq. 32, confirming the accuracy of the approximation for t σ.
As comparison, I also show the corresponding evolution of in Fig. 2, which is calculated from (3) using (27) for the fluid temperature and velocity 3

B. Two Gaussian Hot-Spots
It is possible to generalize the above example to the case of two Gaussian hot-spots.
Specifically, consider an initial pseudo-temperature distribution given by with σ, T init again controlling the width, and height of the hot-spots, and x 0 , x 1 specifying the hot-spot locations, respectively. Taking the distribution at the initial time to be of the form (24), the zeroth order eremitic energy-momentum tensor T µν hermit,(0) is given by a simple superposition of (26) at positions x 0 , x 1 , respectively. For early times, the eigenvalue decomposition of T µν hermit,(0) leads to have been used. This leads to Without loss of generality, one may take the two Gaussian hot-spots be located at x 0 = (x 0 , 0, 0) and x 1 = (−x 0 , 0, 0), respectively. An interesting quantity to consider is the evolution of the elliptic momentum anisotropy e 2 of the system defined as which in the case of x 0 = σ for early times evaluates to (By changing the orders of integration over x, p, it is easy to verify that e 2 [f hermit,(0) ] = 0 for all times, cf. Ref. [38].) A comparison of this result to the numerical evaluation of the first order eremitic expansion is shown in Fig. 2, along with the result for the zeroth order hydrodynamic expansion e 2 [f eq ], demonstrating that the analytic approximation works well for t σ.
The present example is very similar to the case of a deformed Gaussian in two dimensions considered in Ref. [30], where many analytic results were obtained. The main difference to Ref. [30] is that by superposition, the above calculation can be generalized to the case of multiple Gaussian hot-spots without additional complications.

IV. NUMERICAL EXAMPLES WITH BOOST INVARIANCE
Consider now eremitic expansions for particles in a system having boost invariance, as is approximately the case for the high density region of relativistic nucleus-nucleus collisions.
To this end, it is useful to consider a coordinate transformation to Milne coordinates proper time τ = √ t 2 − z 2 and space-time rapidity ξ = arctanh z t . Using x T = (x, y), the Boltzmann equation in coordinates x a = (τ, x T , ξ) becomes where p τ (τ ) = p 2 T + (τ 2 p ξ ) 2 /τ 2 . Assuming that the system is invariant under boosts in the longitudinal direction leads to f = f τ, x T , p T , p ξ , i.e. independent of rapidity. Solution of the characteristic equations for the eremitic expansion to zeroth order gives [37] and the first order result is given by where p τ 0 (τ ) = p 2 T + (τ 2 p ξ ) 2 /τ 2 0 and integration was started at some finite proper time τ 0 . Let us again consider a class of boost-invariant initial particle distribution functions at proper time τ = τ 0 parametrized by f init (x T , p T , p ξ ) = F (p τ 0 (τ )/Λ (x T )), such that the associated energy-momentum tensor in zeroth order eremitic expansion is given by where I changed variables from p ξ = τ −1 p T sinh(Y − ξ) to momentum rapidity Y and used

A. Single Gaussian Hot-Spot with Boost Invariance
Taking the initial particle distribution function to be of equilibrium form F (x) = π 2 Z 3 e −x with Z parametrizing the number of degrees of freedom, and setting Λ( to be a two-dimensional Gaussian leads to simple integral expressions for the non-vanishing components of the energy momentum tensor: where S ≡ τ cosh Y − τ 2 0 + τ 2 sinh 2 Y , I n (x) denote modified Bessel functions and l, m = (x, y). Given the energy-momentum tensor, one may use (5) to calculate the local energydensity and flow vector u a = (u τ , u T , 0) in zeroth order eremitic expansion. These in turn determine C[f hermit,(0) ] in the relaxation time approximation, from which the first-order eremitic expansion (46) may be obtained.

B. Two Gaussian Hot-Spots with Boost Invariance
Similar to the case of Gaussian hot-spots considered in III B, the superposition of two hot-spots located at position x = ±σ amounts to linear superposition of the contributions of individual hot-spots (48) to obtain the energy-momentum tensor, from which the zerothorder eremitic results for , u µ , and eventually f hermit, (1) can be obtained numerically. Of particular interest will be the so-called differential flow coefficients v n [39,40], which for a boost-invariant system at mid-rapidity may be approximated as where p T = p T (cos φ, sin φ). For the case of two hot-spots with boost-invariance, the so-called elliptic flow coefficient v 2 (τ, p T ) may be evaluated numerically. Note that In order to model the flow of QCD partons, I take the number of degrees of freedom to be parametrized by  Given that the hydrodynamic expansion breaks down at large momentum, and that the eremitic expansion breaks down at low momentum, one can expect that the "true" result said to undergo a collision if their respective distance in the transverse x T plane is less than |x T | < σ N N /π, where σ N N 60 mb is the (collision-energy dependent) nucleon-nucleon cross-section at √ s = 5.02 TeV. Each location of a collision is taken to correspond to the location of one Gaussian hot-spot. The sum over these Gaussian hot-spots defines the initial energy-density distribution in the transverse plane, which is successfully used in modern hydrodynamic modeling of lead-lead collisions [2]. The number of nucleons participating in a collision is related to the total entropy of the system, which in turn translates to the number of particles observed in experiment ("multiplicity"). In the following, I will consider mid-central lead-lead collisions corresponding to the 30-40 percent highest multiplicity class.
For this case, the parameter T init was adjusted such that the hydrodynamic evolution with a 5 The hydrodynamic results have been calculated numerically using the same initial conditions and same equation of state using the numerical solver VH2+1 [41] for the hydrodynamic equations.
QCD equation of state [46] gives multiplicities that are consistent with those found in experiment [47]. Unlike the two hot-spot case treated above, for multiple hot-spots found in the Glauber modeling of Pb+Pb collisions the momentum flow coefficients v n for n = 3, 4, . . .
are also non-vanishing in general.
Last but not least, the relaxation time coefficient for QCD is expected to scale as [48,49] up to logarithmic corrections, where α s (Q 2 ) is the QCD coupling that in one-loop running is given by whereΛ 0.376 GeV to match the experimentally determined value α s (Q 2 = M 2 z ) = 0.1184 at the mass of the Z-boson M z 91.18 GeV [50]. I therefore use τ R T = 1 7 α −2 s (p 2 T ) when attempting to make comparisons to QCD (numerically, this implies η/s = τ R T 5 0.13 when evaluating Q = p T 1.5 GeV). It should be emphasized that the choice τ R T α 2 s = 1 7 is arbitrary, and has been made for illustration purposes. Nevertheless, in view of the sizable uncertainty in the value of e.g. η/s from perturbative QCD calculations at scales p T 1.5 GeV, such a choice does not seem entirely unreasonable [51].  high momenta p T ≥ 15 GeV by a type of Padé fit, suggesting a peak in v n (τ f , p T ) for specific values of p T for n = 2, 3, 4. Note that the available information at low and high momenta, respectively, is not sufficient to unambiguously determine the location or height of the peaks in v n (τ f , p T ) .
Since the results shown for v n (τ f , p T ) are for massless partons obtained when the whole system has cooled down below a pre-defined temperature, the results are not directly comparable to experimental data. However, it is tempting to inspect the relevant experimental data on differential flow coefficients for 30-40% Pb+Pb collisions for unidentified hadrons, shown in the rhs panel of Fig. 4. Interestingly, the experimental data seems to exhibit the qualitative features of the above theoretical calculations at low momenta (rise with p T as predicted by hydrodynamic expansions) and high momenta (decrease with p T as predicted by eremitic expansions). Curiously, also the magnitude of experimentally measured v n coefficients at p T < ∼ 2 GeV and p T > ∼ 15 GeV seem to be consistent with theoretical calculations shown in the lhs panel of Fig. 4. Furthermore, note that the ratio v 3 (τ f ,p ⊥ ) ( v 2 (τ f ,p ⊥ ) ) 3/2 1 exhibits near-constant behavior close to unity as a funtion of p T at large momenta in eremitic expansion, similar to what has been observed experimentally [55].

V. HIGH ENERGY p+Pb COLLISIONS
One of the unresolved questions in the context of high energy nuclear collision is the mechanism for the measured sizable v 2 coefficient at high transverse momenta p T > ∼ 10 GeV, cf. Fig. 4. It has been suggested that the measured v 2 coefficient arises from jet quenching, with highly energetic particles (jets) losing more energy when traveling through a longer path length in a medium [57,58]. However, jet quenching seems to be absent in proton-lead collisions, yet the experimentally measured v 2 coefficient exhibits the same behavior as in lead-lead collisions [3], cf. Fig. 5. Eremitic expansions offer a potential alternative explanation for the observed v 2 coefficient, namely through non-hydrodynamic transport of the initial geometry. While the momentum anisotropies in eremitic expansions arise from the dynamics of high energy particles, these particles are nevertheless part of, and flowing with, the medium, as opposed to the modeling of jets, which are by definition treated separately from the medium.
For this reason, I have simulated central p+Pb collisions through Monte-Carlo sampling positions of nucleon collisions from a Glauber model, and using these positions as the initial location of Gaussian hot-spots as explained in the preceding sections. The dynamics encountered in p+Pb is not boost-invariant, but hydrodynamic simulations seem to indi-cate that nevertheless boost-invariance is not a bad quantitative approximation in practice [59][60][61][62]. The results for the momentum anisotropies v n (τ f , p T ) averaged over 10 events for zeroth order hydrodynamic and first order eremitic expansion are shown in Fig. 5. One finds that the same qualitative features as in Pb+Pb emerge: rising v n coefficients at low p T as predicted by hydrodynamics, and falling v n coefficients at large p T as predicted by eremitic expansions. Unlike the case for Pb+Pb collisions, the magnitude for v 2 at large p T for massless partons from eremitic expansions of p+Pb collisions, while non-vanishing, is systematically below the experimentally measured values for unidentified hadrons (rhs panel of Fig. 5). Future studies involving more realistic equations of state and a confinement prescription will be needed in order to decide if eremitic expansion qualify as explanation for the observed v 2 coefficient at high momenta.

VI. SUMMARY AND CONCLUSIONS
In this work, a systematic expansion procedure for relativistic kinetic theory in the large mean-free path regime was considered. This eremitic expansion procedure is complementary to the perhaps more familiar hydrodynamic expansion scheme in that it allows controlled calculations at very high particle momenta, while breaking down at low particle momenta.
Eremitic expansions allow to probe purely non-hydrodynamic transport phenomena since hydrodynamic modes are absent in this approach. Using kinetic theory in the relaxation time approximation as an example, first order eremitic expansions for Gaussian hot-spots with and without boost invariance were calculated. Applications for these calculations to evaluating the momentum anisotropy coefficients v n (τ f , p T ) in Pb+Pb and p+Pb collisions at √ s = 5.02 TeV were presented, and it was found that eremitic expansions qualitatively describe the experimentally measured behavior of flow coefficients at high momenta. Thus, eremitic expansions offer a potential alternative to jet quenching as the source for the measured elliptic anisotropy at high momenta. Nevertheless, eremitic expansions seem to have the potential to become an interesting new tool in the study of relativistic collision systems and the phenomenology of nonhydrodynamic transport.

VII. ACKNOWLEDGMENTS
This work was supported in part by the Department of Energy, DOE award No de-sc0017905. I would like to thank Jamie Nagle, Nicolas Borghini, Ulrike Romatschke and Jürgen Schukraft for helpful comments, and Gowri Sundaresan for collaboration in the early stages of this work.