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 momenta in Pb + Pb and p + Pb collisions at s=5.02\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=5.02$$\end{document} TeV are calculated from the first order eremitic expansion of kinetic theory in the relaxation time approximation.


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 non-hydrodynamic 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 nonhydrodynamic modes, hence at late times bulk transport will be dominated by purely non-hydrodynamic degrees of freedom. While many studies exist that include both hydrodya e-mail: paul.romatschke@colorado.edu namic 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 to 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 freestreaming) 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 aris-ing 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.

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 fourvelocity 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 energy-momentum tensor [33,34] T μν (t, where θ(x) denotes a step-function. From the energymomentum tensor, the local energy density (t, x) and local rest-frame four vector u μ (t, x) can be obtained as the timelike 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 function f eq [ f ] via Eq. (3), where the functional dependence on the nonequilibrium particle distribution f has been denoted explicitly.

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 (zeroth) order in τ R T 1, Eq. (1) leads to 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 energymomentum 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]).

Zeroth order eremitic expansion: ballistic regime
Let me now consider an "eremitic" expansion of the kinetic theory in the regime of large mean free path where C[ f ] 0 or τ R T 1 in the relaxation time approximation. To leading (zeroth) order in eremitic expansion, Eq. (1) leads to This equation can be solved analytically using the method of characteristics: leading to where f init (x, p) is the particle distribution function at initial time t = 0. The energy-momentum tensor for the zeroth order eremitic expansion is given by (4), which requires specification of f init . Some example cases will be considered below.

First order eremitic expansion
The first order eremitic expansion for the distribution function is given by with in the relaxation time approximation where the freestreaming 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 correction 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 noninteracting particles, and corrections from interactions build up coherently along the characteristic for small gradients.

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 energymomentum 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, 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 points 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 twopoint 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 non-hydrodynamic 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.

Analytic examples
While powerful methods exist to solve the Boltzmann equation (1) numerically, the strength of the eremitic expansion lies in the possibility of obtaining analytic (or at least semianalytic) results. To this end, let me point out some examples where analytic treatments are possible. All of these examples will be within a class of initial particle distribution 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.

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 where a ≡ t|x−x 0 | σ 2 . The pseudo-temperature T and flow vector u μ corresponding to this energy momentum tensor can be obtained via eigenvalue decomposition, cf. Eq. (5). One finds 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 x f 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 finds 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.
in Fig. 2, which is calculated from (3) using (27) for the fluid temperature and velocity 3

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 hotspots, 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 2 Numerical algorithms employed for this work are publicly available for download from [35]. 3 Note that Eq. (27) are not the solution to ideal fluid equations of motion ∂ μ T μν fluid,(0) = 0 because σ μν = 0, but they are very close, cf. the discussion in Refs. [36,37].
where the shorthand notations e 0 = e 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.

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 shorthand notation

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 ( T /(8σ 2 ) 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 energy-density 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.

Two Gaussian hot-spots with boost invariance
Similar to the case of Gaussian hot-spots considered in Sect. 3.2, 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 zeroth-order 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 4 where p T = p T (cos φ, sin φ). For the case of two hotspots with boost-invariance, the so-called elliptic flow coefficient v 2 (τ, p T ) may be evaluated numerically. Note that v 2 (τ, p T ) [ f hermit,(0) ] is time independent and vanishes identically if v 2 [ f init ] = 0, consistent with the expectation that no elliptic flow is generated for non-interacting particles, while v 2 (τ, p T ) is in general time-dependent in first order eremitic expansion. In order to model the flow of QCD partons, I take the number of degrees of freedom to be parametrized by where here and in the following the number of colors and flavors are chosen as N c = 3 and N f = 2, respectively. The initial temperature is chosen as T init 0.535 GeV at τ 0 = 0.25 fm and the evolution is allowed to continue until a proper time τ f at which the local temperature has dropped below T c = 0.170 GeV everywhere in the system (sometimes known as constant proper-time decoupling). Numerical results for v 2 (τ f , p T ) are shown in Fig. 3 for the first order eremitic expansion for τ R T = 2.5 as a function of transverse momentum. As can be seen from this figure, the eremitic expansion suggests a near-constant behavior of the elliptic flow coefficient at large momenta. This matches the results found by Borghini and Gombeaud for the case of a two-dimensional deformed Gaussian hot-spot without boost invariance [30].
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 Boltzmann equation (1) would rise according to the hydrodynamic result at low momenta, and saturate at a constant value according to the eremitic result at high momenta for constant τ R T . It should be straightforward to test this expectation using numerical solutions to the Boltzmann equation [42,43] for a two hot-spot case.

High energy Pb+Pb collisions
The same techniques as for the two hot-spot case may be used to model high energy nuclear collisions, such as Pb+Pb collisions at √ s = 5.02 TeV. This is because the so-called Glauber model [44,45] provides initial conditions for the matter distribution deposited after the collision as a sum over Gaussian hot-spots corresponding to the locations of collisions of the individual nucleons (see Ref. [2] for a recent review of relativistic nuclear collision modeling). For the purpose of this work, hot-spot locations are generated by first Monte-Carlo sampling nucleon positions for two lead nuclei from a suitably normalized Woods-Saxon probability distribution function ρ(x) ∝ (1 + e (|x|−r 0 )/a 0 ) −1 with r 0 = 6.62 fm and a 0 = 0.546 fm. Random sampling of impact parameters for the collision of two lead nuclei, nucleons are 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) nucleonnucleon 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  [52][53][54] for momentum anisotropy coefficients for unidentified hadrons in Pb+Pb collisions at √ s = 5.02 TeV in the 30-40% multiplicity class 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 % highest multiplicity class. For this case, the parameter T init was adjusted such that the hydrodynamic evolution with a 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 hotspots 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].
Results for the anisotropic flow coefficients v n (τ f , p T ) , averaged over 10 events of initial hot-spot locations in the 30-40% multiplicity class of Pb + Pb collisions are shown in Fig. 4. Again, results from zeroth order eremitic expansion for v n (τ f , p T ) vanish identically, while first order results differ significantly from zero for v 2 [ f hermit, (1) ]. For p T 15 GeV, the eremitic expansion seems to break down for this choice of τ R since the correction f hermit, (1) approaches the leading-order result f hermit,(0) . At high p T , v n (τ f , p T ) appears to approach a constant times the inverse of τ R T , just as what was found in the two hot-spot case above. For τ R T of the form (51), this implies v n (τ f , p T ) falling as the QCD coupling constant squared for large p T . For comparison, results for v n (τ f , p T ) [ f fluid,(0) ] are also shown in Fig. 4. The hydrodynamic gradient expansion breaks down at high p T , but v n (τ f , p T ) [ f fluid,(0) ] exhibits the rising trend familiar from full hydrodynamic simulation studies [12]. The hydrodynamic gradient results 6 at low momenta p T ≤ 2 GeV can be connected to the eremitic curves at 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   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].

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 indicate 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.

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 nonhydrodynamic 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.
Many generalizations and validations of the present work are possible. For instance, the second order correction to eremitic expansions should be straightforward to calculate for many of the examples given in this work. The quantitative reliability of eremitic expansions should be checked by direct comparison to full numerical solutions of the Boltzmann equation. The application of eremitic expansions to relativistic collision systems should be made more realistic by including a QCD equation of state and a hadronization procedure.
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.