Hydrodynamic and Non-hydrodynamic Excitations in Kinetic Theory – A Numerical Analysis in Scalar Field Theory

: Viscous hydrodynamics serves as a successful mesoscopic description of the Quark-Gluon Plasma produced in relativistic heavy-ion collisions. In order to investigate, how such an effective description emerges from the underlying microscopic dynamics we calculate the hydrodynamic and non-hydrodynamic modes of linear response in the sound channel from a first-principle calculation in kinetic theory. We do this with a new approach wherein we discretize the collision kernel to directly calculate eigenvalues and eigenmodes of the evolution operator. This allows us to study the Green’s functions at any point in the complex frequency space. Our study focuses on scalar theory with quartic interaction and we find that the analytic structure of Green’s functions in the complex plane is far more complicated than just poles or cuts which is a first step towards an equivalent study in QCD kinetic theory.


Introduction
Hydrodynamics play a central role in describing the collective behavior of macroscopic systems in the real world.The premise of hydrodynamics is that it describes the late time and long wavelength limit of a system.The Quark-Gluon-Plasma (QGP) found in ultrarelativistic Heavy Ion Collisions (HIC) is very short lived and the average system size is small, yet data collected from heavy-ion collisons at RHIC and LHC suggests that the space-time dynamics can be well modelled by hydrodynamic theories [1][2][3][4].Considering that the QGP is initially far from equilibrium and is subject to large gradients, naturally the question arises on what time and distance scales hydrodynamic theories can provide a realistic description of HICs [5][6][7][8][9][10][11][12][13][14][15][16][17].Since in practice, hydrodynamic theories are based on expansions around local thermal equilibrium [18][19][20], it has long been believed that proximity to equilibrium is a necessary criterion for the applicability of a fluid dynamic description.However, in recent years various studies have indicated that -at least for certain classes of microscopic systems -viscous hydrodynamics can provide a rather accurate description even when the system is significantly out of equilibrium, featuring for example pressure anisotropies of order unity [21,22].Hence it has become customary in the field of high-energy nuclear physics, to carefully distinguish equilibration from hydrodynamization, which merely refers to the applicability of viscous relativistic hydrodynamics, and we refer the interested reader to [23][24][25][26][27][28][29][30] for recent perspectives on this issue.One field of study that emerged relatively recently to tackle the question of (in)applicability of hydrodynamics, consists of analyzing so called hydrodynamic and non-hydrodynamic modes [31][32][33][34][35]. Hydrodynamic modes, which correspond to isolated singularities in the Fourier transformed evolution equations of hydrodynamic theories, are calculated by linearly perturbing evolution equations and calculating the system response.These modes fulfill the hydrodynamic limit lim k→0 ω = 0, corresponding to conserved quantities.Non-hydrodynamic modes in turn do not fulfill this large wavelength limit and are present on any length scale, albeit varying in importance [36][37][38].Although they are called non-hydrodynamic, they are also found in hydrodynamic theories such as Müller-Israel-Stewart [39][40][41][42][43].So when these nonhydrodynamic modes are ever present, the question still stands when they dominate the dynamics of the system.The regime, where they are non-negligible, will help understand the boundaries of the applicability of hydrodynamics.A common procedure to calculate the dynamics of the QGP are multi-stage models that use kinetic theories to describe the pre-equilibrium phase which then transition smoothly into the hydrodynamic regime [44][45][46][47][48][49][50][51].Kinetic theories are also sometimes used to directly derive new hydrodynamic theories, such as the DNMR hydrodynamic theory [18,43,52,53].Because of this ability to capture both non-equilibrium and non-hydrodynamic behavior kinetic theory is a fitting model to find non-hydrodynamic modes that are so natural to the system.Solving the Boltzmann equation analytically is not feasible for non-trivial kinetic theories, thus we have to solve the problem numerically.The problem then is to find the modes we are looking for.Analytically the modes are extracted as non-analytic structures of the Green's functions of the energy-momentum tensor, for example poles or cuts, but numerically it is hard to do the same.Going into frequency space becomes problematic numerically when the integral is not restricted to the real frequency axis.Therefore we need some other method to determine the complex structure of energy-momentum Green's functions.
In this paper, we present a method to calculate eigenfunctions and eigenvalues of the Boltzmann equation with non-zero gradients k in a discretized fashion by using moments of the distribution function.Based on the method described in Sec.2, we can discuss the location of singularities and the analytic structure of Green's function in the complex frequency plane.In Sec.3 the procedure is tested within the Relaxation Time Approximation (RTA) and subsequently applied to scalar ϕ 4 kinetic theory in Sec. 4. The RTA can be compared to exact calculations and reproduce known results.The scalar calculations expand the known knowledge of zero gradient results into a finite k regime.We find an analytic structure that goes beyond poles and cuts which was predicted by some previous works [33,34].

Methodology
Kinetic theory is an effective mesoscopic theory, where the time evolution of the phasespace distribution f (p, ⃗ x, t), where p = (p 0 , ⃗ p) = (p, ⃗ p) is the four-momentum and x = (x 0 , ⃗ x) = (t, ⃗ x) are the space-time coordinates, is governed by a Boltzmann type equation.Within kinetic theory, it is comparatively straightforward to calculate the time evolution of the phase-space distribution, as well as the behavior of macroscopic quantities, such as e.g.retarded correlation functions of the energy-momentum tensor G(t), which can be obtained from moments of the phase-space distribution.In particular, there are several numerical studies of kinetic theories, including scalar ϕ 4 -theory or even QCD kinetic theory [34,38,45,48,49,[54][55][56][57][58], which explore the real-time dynamics of the system.However, as mentioned before, the problem of analyzing the structure of response functions in the complex frequency space originates from the need for a Laplace transform which becomes numerically ill behaved for complex frequencies with Im(ω) < 0, as the above integral does not converge beyond the first singularity in the lower half-plane.While for analytic solutions of the energy momentum Green's function it is still possible to calculate the Laplace transform directly, as has been done e.g. in RTA [31,33], this is evidently not possible with numerical data for G(t), and thus a different approach is required.

From Collision Integral to Matrix
In our method we will discretize the Boltzmann equation in order to calculate discrete spectra of eigenvalues and their eigenfunctions, which tell us exactly how the system can respond to perturbations.The general form of the Boltzmann equation is where C[f ] is the the collision operator, also called collision integral or collision term in the following, of the distribution function f .Throughout this work we use the mostly minus metric.We will study linear perturbations on top of an equilibrium background, which is defined by ultrarelativistic bosons with zero mass obeying the equilibrium distribution Our perturbation is then of the form 3) The collision integral, in this study only for 2 ↔ 2 scatterings, also has to be linearized around the equilibrium background, yielding where |M | 2 is the scattering matrix element squared and δF is the linearized statistical factor Note that we abbreviated the distribution functions in a way to show their dependence on a specific momentum, e.g.n 1 = n(p 1 ).For RTA the collision term takes a much simpler form, which is discussed in the respective section.For the spatial part of the Boltzmann equation we switch into the Fourier space by transforming according to where gradients are replaced by the wave vector ⃗ k Fourier transformed quantities are marked with a k in the subscript.In this work we will focus only on perturbations in the sound channel.For this we will initially perturb the system by a temperature perturbation, such that where δT 0 is the magnitude of our initial perturbation in T .Since our background is isotropic and we will consider only scalar perturbations the perturbations will only depend on the absolute value of the wave vector | ⃗ k|.Without loss of generality we orient the wave vector along the z-axis to receive the linearized Boltzmann equation where cos θ = ⃗ k•⃗ p kp is the longitudinal momentum angle and k = | ⃗ k|.The collision and gradient term will be discretized in order to solve it numerically.How we then extract and evaluate the eigenvalues depends on the diagonalization of said discretized collision integral. 1The discretized version of a linear operator is a finite dimensional matrix, which is easily decomposed into its eigenvectors.The eigenvectors of a matrix can be divided into right and left eigenvectors, which are usually not the same.Hermitian or symmetric matrices do not have this separation of eigenvectors.In the case of k = 0 the evolution matrix would be symmetric, but the addition of gradient terms in the matrix makes it neither hermitian nor symmetric, meaning we have to respect the distinction of left and right eigenvectors.If ⃗ a i is the right and ⃗ b i the left eigenvector to the eigenvalues λ i of the matrix C they satisfy The eigenvector can also be used to denote the diagonalization of C via the modal matrix P where D = diag(λ i ) is the diagonal form of C. The modal matrix and its inverse are related to the eigenvectors via the following relations Since the discretized version of our distribution function will be a vector let us call it ⃗ f here to illustrate the basic principle of the method.After the discretization process the Boltzmann equation will transform into a vector matrix equation which is solved by ⃗ f (t) = e Ct ⃗ f (0).Using the diagonal form of C and the representation of the modal matrix we receive This form gives the solution to the Boltzmann equation as a sum of contributions from each eigenvalue with some weight, that is dependent on the initial condition, which makes it incredibly easy to discuss the individual influence of each eigenvalue.From there it is also straight forward to go into frequency space .16)This form shows that each eigenvalue directly corresponds to a singularity in frequency space located at ω i = −Im(λ i ) + iRe(λ i ), revealing the complex structure of the theory.

Moments of the Distribution Function
We just showed how an operator equation may be rewritten into a much simpler and easier analyzable form.We will now discuss how the distribution function and collision integral are transformed into a discrete space.We start by the introduction of moments of the distribution function where the w i (p) are some weight functions.This way one can construct a vector ⃗ N , where the components are the individual moments N i .Our approach follows [56,59] in the construction of their weight functions.Since we will only study sound mode excitations, it is sufficient to discretize the momentum p and the polar angle cos θ.We do this with the use of so called wedge functions w i (x), which are defined as where the x i are discrete grid points of the quantity that is discretized in the wedge function.
The momentum grid contains momenta from p min = 0 to p max evenly spaced and N p in number.The cos θ grid contains N cos θ points and goes from cos θ min = −1 to cos θ max = 1.This means the grid has a total size of N tot = N p N cos θ entries.The wedge functions fulfill following relations So the moments of the distribution function, also sometimes called wedge moments in the following, are then defined as where i is a combined index in the form of i = i θ + i p N cos θ .Using the properties of the wedge function one can gather information about particle number δn k = δJ 0 k , energy δe k = δT 00 and longitudinal momentum δπ k = δT 03 k by simply summing over the wedge moments in the form of where we used the particle four current δJ µ k and the energy momentum tensor δT µν , calculated in kinetic theory by taking moments of the distribution function ) We still need a way to calculate the matrix described in our method.We do this in the same moment approach and take moments of the collision integral In order to construct a matrix from this we have to invert the wedge moments back into the form of a distribution function.We do this by approximation of the original moment integral.When we take δf k (p,t) feqe p/T to be constant between nodes we can rewrite the Integral as where we introduced a new area function A i that is calculated as Eq.(2.25)only holds for the surrounding of the i−th node point, thus in order to recover the whole distribution function we have to sum over all wedges where we introduced the "Cardinal Function" K i (p) With the introduction of this Cardinal Function the Collision integral moments become linear functions of the distribution moments and thus we can rewrite the equation as a matrix vector multiplication The matrix entries are then easily calculated with the functional derivative

.30)
With these preparations we can rewrite the Boltzmann equation for k = 0 as a matrix equation which is the desired form.
In addition to the collision integral the Boltzmann equation contains parts with nonzero k.These parts also have to be translated into the moment space.Like the collision integral they are linear in the distribution function and thus we can calculate a matrix in a similar fashion.The moments of the gradient contribution are given as from which we calculate the matrix M as Note that the k dependent part is independent of p, thus one only needs to construct angular wedge moments.Then the matrix elements can be written as where M cosθ ij are now the angular wedge moments of the gradient part.The Kronecker Delta part is there to assure that no momenta are mixing.The angular wedge moments are defined as which can be easily calculated by hand.The A cos θ j are area functions for the angular wedge moments just like the recent A j but are much simpler ∆ cos θ = cos θ i+1 −cos θ i is the distance between each angular grid point, which is constant throughout because we choose an evenly spaced grid.This allows for a fast calculation of the gradient contribution of the equation, since only k has to be multiplied to known matrix elements.In contrast the calculation of the collision integral matrix is done with a Monte Carlo scheme, where we update all matrix elements simultaneously with each sampling, such that particle number and energy momentum conservation are ensured with the help of the wedge function properties in Eq. (2.19).
Both collision integral and gradient term are combined as matrices in the ordinary differential equation with the solution Calculating the collision integral matrix in high precision takes some time but can be saved for reuse in the same discretization because for different k only the gradient contribution changes.This also saves a lot of computation time where the remaining computation time is due to the numeric calculation of the eigenspace.

Observables and Complex Frequencies
We now have the means to completely discretize the Boltzmann equation including gradient contributions.In the following we discuss the construction of observables from moments and the possibility to go into the complex frequency plane.
Observables that are linear moments of the distribution function δf , like energy, can be easily retrieved (see Eq.(2.21)) from wedge moments via a sum of moments in addition to some weight.We can write a representation of the energy in our moment space as a vector because we get the energy by forming the scalar product The Green's function is the time evolution of an observable as a response to an initial perturbation of some sort.This initial condition is already encoded in the time evolution of ⃗ N (t) as ⃗ N (0)(see Eq.(2.38)).Thus we then define the Green's function as the scalar product of the corresponding observable vector with the distribution vector The full form of the Green's function can then be calculated using the solution to the ODE in Eq.(2.15) as Each eigenvalue is behaving as a complex exponential function with a certain weight, we call this the contribution or residue to the Green's function.Let us define the contribution of an eigenvalue as µ i The Laplace transform of the Green's function is then With this we have a direct way to study the Green's function as a function of time or in the complex frequency plane, where each eigenvalue is an individual singularity.

Scaling Behavior
Since we are studying the theory with an approach to hydrodynamics in mind it makes sense to use hydrodynamic scaling properties.In first order viscous hydrodynamic theory one can rescale the wavenumber k and time t by the viscosity to receive a universal description independent of viscosity [38].Thus we define the rescaled wavenumber, time and frequency as Notice that the combined kt is still equal to kt and ω is defined according to t, being the Fourier counterpart of t.Except for the scaling variables, all dimensionful variables, like p max , are expressed in units of the only dimensionful scale, the background temperature T .

Benchmark with RTA
Before we apply the method to scalar theory we should put it to the test within a well studied kinetic theory, the RTA [31,33].In RTA all excitations of different momenta decay equally with a rate τ R , the relaxation time, correlating to a branch cut in the retarded Green's function.The RTA Boltzmann Equation is of the form where the local rest frame velocity u µ is determined by Landau matching as eigenvector of the background energy momentum tensor T µν eq u ν = e eq u µ . (3.2)

Analytical Results
To study the linear response to a perturbation we first use the linearized version of Eq.(3.1) For positions we switch into Fourier space, without loss of generality orienting the wave vector ⃗ k only along the z-axis.For time we do a Laplace transform instead of Fourier transform since we have strictly positive times and want to include initial conditions.This results in where δf k (p, t = 0) is the initial condition of the distribution function.The perturbed equilibrium distribution is given by gradients of temperature and velocity emerging from a change in δf as With this we can get the solution to Eq.(3.4) Via Landau matching we can relate the perturbed temperature and velocity to δT µν as δe k e = 1 4e ) Since our ⃗ k only lies in z-direction the only relevant components of the Energy-Momentum Tensor are δT 00 k = δe k and δT 03 k .The energy and velocity perturbation are related to the Green's functions we want to find by From the above equations one receives coupled equations for δe k and δT 03 k .The details of the solution for those are found in the Appendix A. When we define the perturbed energy and longitudinal momentum are Evidently, the above Green's functions feature a logarithmic branch cut extending between the branch points ω = −k − i/τ R and ω = k − i/τ R in the complex frequency plane.By expanding the inverse of the Green's function G 03 00,k to second order in ωτ R one also finds a pair of hydrodynamic poles for small frequencies and gradients, while for large gradients they disappear behind the cut as discussed in detail in [31].

Results from Numerical Method
Now that we have obtained the analytic solution, we can use it to benchmark our numerical method.Within the numerical method, we write the RTA Boltzmann equation as which is then expanded in the moment approach as discussed in detail in Appendix B. We can now compare the results from our numeric method with the known analytic results in RTA.In Fig. 1 we present the real and imaginary part of the energy Green's function plotted as a function of real frequency.The colored lines represent the numeric results from our new approach and the black dotted lines are the prior calculated analytic form.The numeric approach fully reproduces the analytic Green's function over the whole range.As increasing the gradient does not change this fact, we conclude that the new approach is a suitable method to calculate kinetic theory response functions.Of course our main goal is to analyze the analytic structure in the complex frequency plane.Hence we still need to check if the new approach reproduces the well studied RTA structure, namely a cut and two poles, for complex frequencies.In Fig. 2 we present each eigenvalue in the complex frequency plane for three gradients k.We see clearly two isolated modes on top which correspond to the two well known hydrodynamic poles.The eigenvalues below form a line, which represents a branch cut in discretized fashion.The cut is located at ω = −0.2iwhich coincides with the expected RTA cut located at ω = −i/τ R as the viscosity in RTA is given by η/s = 1 5 τ R T [55,60].When k increases the hydrodynamic poles wanders towards the cut until they further disappear behind the cut.This behavior is also known from analytical results of RTA [31].The deviations from a straight line occurring at higher k stems from discretization effects, which are discussed in the next chapter.Thus we conclude that the new approach shows the analytic structure clearly in the complex plane, which accurately describes the analytic knowledge.

Effects of Discretization
As we saw in the previous chapter the eigenvalues of the RTA should be two hydrodynamic poles and one cut located at −i/τ R .Additionally to the cut there are some aberrations at the edge of the cut that lean upwards in imaginary axis direction.To check if the deviations from the expected behavior are under control we analyze the effects of the discretization in the following.The increase of N p converges very quickly into a final solution with N p = 64 both in the time evolved Green's function G(t) and the eigenvalue picture in the complex frequency plane.Thus we will only discuss the influence of N cos θ by varying N cos θ and holding N p constant at N p = 16.The results are presented in Fig. 3, where the eigenvalues for various discretizations are plotted in the complex frequency plane for k = 1.We can directly see that by increasing N cos θ we push the deviation of the eigenvalues at the edges of the cut down towards the cut, telling us that these deviations from the expected cut are purely discretization effects.We also have to look at the effects of N cos θ on the time dependent Green's function G(t).The results of this are found in the inset plot of Fig. 3, where the same coloring as the eigenvalue plot is used to display different discretizations.For N cos θ = 64 the time evolution already reached its limit because further increase N cos θ does not change the function or at least only minimally for very late times.We conclude that discretization effects can be erased in the time evolved Green's function entirely by using a feasible discretization.The eigenvalue picture however requires a high N cos θ to come close to the analytic cut expected in RTA, which is not numerically feasible under the consideration of also using N p = 64.The effects of the discretization are rather easy to discern from the physics here, thus we use a discretization of N p × N cos θ = 64 × 64 and keep in mind that the slight arcs of the cut are artificial.

Scalar Theory
Previously we have shown that our method is suitable to analyze kinetic theories in both time and frequency domain.Now we will apply the method to a kinetic theory where the behavior for non-vanishing gradients ( k ̸ = 0) has not been studied yet, the scalar field theory for quartic interaction, also known as ϕ 4 theory.Studies for k = 0 have shown a complex structure in the form of a cut on the imaginary axis, meaning that except for the conserved quantities all excitations of the system decay on different time scales [34].With the addition of gradients one expects the emergence of additional complex structures, including two isolated poles, which coincide with hydrodynamic poles for very small gradients k ≪ 1.These hydrodynamic poles obey a dispersion relation, which for example in second order hydrodynamic theory [41] is given by ω where c 2 s = 1/3 is the speed of sound and τπ ≈ 6.1 is a second order transport coefficient [55].The Lagrangian of scalar ϕ 4 theory is defined by the real scalar field ϕ and the coupling strength λ s.t.following [34] the collision integral in kinetic theory then takes the form As explained before we want to use scaling variables in order to describe the system universally and independent on viscosity.The viscosity in scalar ϕ 4 theory has been calculated numerically multiple times [34,54,61] and for this work we choose the result from [61], which gives the viscosity in terms of temperature T and coupling λ as To write this in a usable manner as viscosity η over entropy density s we use the thermodynamic relations of a massless gas of ultrarelativistic bosons [19] and receive In this form it is also straightforward to see that, by use of the scaling variables defined in Eq.(2.46), the coupling strength λ can be completely scaled out of the Boltzmann equation Eq.(2.9), such that the results presented in the forthcoming sections are independent of the coupling strength.We emphasize at this point, that the following results are obtained with an effective kinetic description of the scalar ϕ 4 quantum field theory.By using an effective kinetic description, which can be derived from a combined weak-coupling and gradient expansion in quantum field theory (see e.g.[62]), we impose at least some restrictions to the possible eigenvalue picture in the complex frequency plane, such that for example the real part of the complex eigenvalues is restricted to the range −k < Re(ω) < k [33].Hence it is possible that the analytic structure of the real quantum field theory will differ from our kinetic theory results, as microscopic information is lost upon construction of the kinetic description [63].

Spectrum for k=0
Before we address the analytic structure of the Green's functions at finite wave-number, the first result we should reproduce concerns the analytic structure of the Green's function in the absence of gradients ( k = 0).As already shown in [34] the Green's function for vanishing gradients exhibits a branch cut located on the imaginary axis.In order to analyze this behavior, we investigate the spectrum of eigenvalues of the evolution operator, which in the absence of gradients only consists of the collision operator.Our results show three distinct zero modes which correspond to the conservation of energy, momentum and particle number in the sound channel.Other than that the spectrum only contains eigenvalues which are located on the imaginary frequency axis.In order to analyze whether this spectrum is discrete or continuous, one has to investigate how the eigenvalue density and their respective contribution to physical observables behaves in the continuum limit.
As objects to study these properties with we choose the eigenvalue density ρ(ω) and residue density ρ µ (ω).Since all eigenvalues have no real frequency Re(ω) the densities are calculated with respect to the imaginary frequency Im(ω).The eigenvalue density is calculated by counting all eigenvalues λ i in a frequency bin of size ∆ω with the help of the step function Θ j (ω i , ω) which is defined as where ωi is the frequency of each eigenvalue, as seen in Eq. (2.16).Normalizing this by the total number of eigenvalues, ignoring the three zero modes, we receive the eigenvalue density In order to study a non-trivial residue density a suitable observable has to be chosen.All conserved charges have vanishing residue for ω ̸ = 0, hence a non conserved charge has to be used.As observable we consequently choose the momentum squared The residue µ i of this observable per eigenvalue is given by Eq. (2.43).We have to sum this, analogously to the eigenvalue density, in a frequency bin ∆ω and normalize it.Here we don't normalize by the total number of eigenvalues but the total residue of eigenvalues with ω ̸ = 0, which is subject to small variations due to the change in discretization2 .If we define this total non-conserved residue as we can write the residue density of the non-conserved observable as where we also ignore the three zero modes in the sum.Since for our choice of initial condition (c.f Eq.(2.8)), the angular structure is trivial for k = 0, the discretization analysis is carried out by varying the number of momentum points N p , which refines the discretization and the momentum cut-off p max , which can be used to explore the addition of high momentum modes.We illustrate the results in Fig. 4 by plotting the eigenvalue density ρ(ω) in the bottom panel and the residue density ρ µ (ω) in the top panel.Each discretization is presented by a different color.The first thing to notice in the bottom panel is that the increase in p max corresponds to an addition of eigenvalues in the low frequency regime close to the origin, such that the number of eigenvalues decreases in higher frequency regions compared to discretizations with the same N p .Vice versa when we fix p max and increase the number of momenta N p the density seems to approach a smooth continuum limit for each value of the momentum cut-off.Since the variation in p max induces strong changes in the eigenvalue density we have to check if this influences the physical properties of the system.For this see the top panel of Fig. 4.Here we also see that additional low frequency eigenvalues get added when we increase p max but all of them have exponentially small residue, thus making them negligible for the linear response of the system.The increase in N p has the same influence on the residue density as on the eigenvalue density.By simultaneously increasing p max and N p , the residue density seems to approach a genuine continuum limit with better and better discretization.We conclude that in the continuum limit N p , p max → ∞ the eigenvalues are continuously distributed along the negative imaginary frequency axis, representing a branch cut, which is in line with the prior works of Moore [34].While the branch cut terminates only at the real frequency axis, the contribution of points close to the real axis to physical observables is exponentially suppressed, as the dominant contribution originates from points with Im(ω) ∼ 0.1.

Effects of discretization
Now that we have established that the zero gradient spectrum is a branch cut we will discuss the effects of adding gradients to the Boltzmann equation.Before we discuss the results we should analyze how the discretization plays into it.We denote a certain discretization by N p × N cos θ .For the results later in this work we choose a discretization of 64 × 64 and p max /T = 16.In Fig. 5 we present the results for three different discretizations for k = 0.8 in comparison with the final 64 × 64 and p max /T = 16 discretization.Eigenvalues of the evolution matrix are plotted in the top row and bottom left panel, with the hydrodynamic eigenvalues highlighted as big dots because they carry the largest individual residue for k = 0.8.By increasing the number of angles we see that we get a finer distribution of eigenvalues in the real range of frequencies, as seen in the top right panel.Equivalently we see a refinement in the imaginary frequencies by increasing the number of momenta, as seen in the bottom left panel.Although this yields more eigenvalues the overall picture of the complex plane stays the same, especially the positions of the hydrodynamic poles do not move.When we increase the maximum momentum cut-off, as in the top left panel, we see that we gather more eigenvalues closer to the real axis, as we have already seen for the case of k = 0.But again this does not change the position of hydrodynamic eigenvalues.Despite the various refinements in the complex frequency plane, the real time Green's functions, plotted in the bottom right panel of Fig. 5, shows the same behavior, where all curves of various discretizations lie on top of each other.This discretization test shows that the discretization of 64 × 64 with p max /T = 16 is well suited for the study of the complex structure of scalar theory.The further increase in discretization only added eigenvalues with small residues which did not change the complex structure at all and left the eigenvalues with the largest residues untouched.The physical description of observables, represented by the Green's function, remains completely unchanged upon further refinement.

Spectrum for finite gradients
When one switches to non zero gradients the appearance of hydrodynamic modes is one of the first things happening for very small k.The zero modes observed for k = 0, responsible for conservation laws, become hydrodynamic modes in the finite k case.All the other modes describe the non-hydrodynamic behavior.Since in our calculation we only receive the location of eigenvalues λ i and their contribution µ i to the Green's function, we first need to come up with a definition for what we call hy-drodynamic mode and how we isolate it in our data.We define the hydrodynamic modes, where possible (more on that later), as the complex conjugated pair of modes with the largest residue, which is obviously true for small k.By identifying these modes for various k, we can then obtain the dispersion relation ω( k) of hydrodynamic sound modes, which is shown in Fig. 6 for both scalar theory and RTA as a comparison.While at small k, the dispersion relations for the real (left) and imaginary (right) part of the sound mode agree with second order viscous hydrodynamics, sizeable deviations start to occur for k ≳ 0.15, and increase with increasing gradient strength, which is in line with the results reported in [38].Beyond k ≈ 1.2 the residue of the hydrodynamic mode extracted from scalar theory shrinks dramatically, as the mode disappears into a continuum of non-hydrodynamic excitations, and it is no longer meaningful to distinguish this mode from other excitations of the system, which is the reason that the curves in Fig. 6 terminate at this value of k.Strikingly, a very similar behavior can also be seen in the Relaxation Time Approximation, where the hydrodynamic modes exhibit essentially the same dispersion relations up to k ∼ 0.9, where the hydrodynamic modes in RTA disappear behind the branch cut.Now that we have established the behavior of the hydrodynamic sound mode, we continue to further analyze the behavior of the Green's functions and in particular investigate the structure and impact of non-hydrodynamic excitations.In order to perform this analysis, we monitor the behavior of the Green's functions in the real-time and real and complex frequency space, while varying the wave-number k, which characterizes the magnitude of spatial gradients.Our results are compactly summarized in Figs.7 and 8, where each row shows the behavior of the Green's function for a fixed value of k, with increasing k from top to bottom.In each row, the left most column shows the behavior of the real time Green's functions G( t) (see Eq.(2.42)), along with a decomposition into the contributions from the hydrodynamic sound mode (i.e. the complex conjugated pair of modes with the largest residue) and the non-hydrodynamic modes (i.e.all other), which can be reconstructed by summing only the contributions of the respective modes in Eq.(2.42).In the second column from the left, we depict the real and imaginary parts of the Green's function G(ω) (see Eq.(2.44)) as a function of real frequencies ω/ k.The behavior in the complex frequency plane is elucidated in the third column of Figs. 7 and 8, where black circles show the location of individual eigenvalues calculated from the matrix, which are scaled in size by their contribution µ i to the Green's function, while the color code in these plots corresponds to the logarithm of the absolute square of the Green's function. 3Finally, the right most plot is a histogram showing the summed contributions ρ µ (ω) from non-hydrodynamic modes as defined in Eq.(4.10).Since there can be both positive and negative contributions, we distinguish them by a color coding, where positive contributions are plotted in red and negative contributions in blue.Contributions of the non-hydrodynamic modes are further compared to the contribution of the the hydrodynamic sound mode, which is indicated by a green bar.
When considering the behavior of the real-time Green's function G(t) (left column), we find that for small gradients k ≪ 1 the evolution is almost purely hydrodynamic, meaning that only the complex conjugated pair of modes with the largest residue play a significant role for the time evolution of the system.With increasing k, the contribution of non-hydrodynamic modes becomes visible at early times, but exhibits a much faster decay than the hydrodynamic contribution.When k ∼ 1, the contributions from hydrodynamic and non-hydrodynamic modes become of comparable size, until for k ≈ 1.2 the contribution of the hydrodynamic mode begins to disappear as non-hydrodynamic modes start to dominate the behavior of the real-time Green's function.By k = 1.5, corresponding to the largest value shown in Fig. 8, it is then no longer possible to identify a clear signature of a hydrodynamic mode in the system.Next we consider the Green's function G( w) as a function of real frequency ω (second left column), which for small gradients ( k ≪ 1) features two distinct peaks in its real part and at the same position two steep inflection points in the imaginary part.Clearly, these features are a sign of the contribution of two very distinct regions in the complex frequency plane, in this case the hydrodynamic poles, which for k ≪ 1 are located close to real frequency axis.By increasing k the influence of the hydrodynamic poles diminishes, as their residues shrink and they move away from the real-frequency axis, leading to a broadening of peaks and inflections.Nevertheless, both the real-time and the real-frequency Green's functions exhibit a remarkably smooth behavior with increasing k, and do not immediately indicate a transition from a hydrodynamic to a non-hydrodynamic regime.Now, if we consider the behavior of the Green's function in the complex frequency plane (right two columns), for small k ≪ 1 we can clearly see the dominance of the hydrodynamic modes, which correspond to simple poles in the complex frequency plane with large residues.However, in addition to the hydrodynamic poles, additional singularities indicated by black circles appear throughout the region where − k < Re(ω) < k and Im(ω) < 0, which most likely signals the presence of an entire region of non-analyticity as discussed in [33] for the momentum dependent relaxation time approximation.Despite the fact that the singularities cover a large region in the complex frequency plane, it appears that the dominant contribution to the Green's function of the energy-momentum tensor seems to be located around Im(ω) ≈ −0.2, as can be seen from the residue density in the right most panel, which reminiscent of the k = 0 spectrum in Fig. 4 is also strongly peaked around a single imaginary part of the frequency.Hence, the typical time scale for the relaxation of contributions of non-hydrodynamic modes to the energy-momentum tensor is still given by the inverse of this characteristic frequency ∼ 1/Im(ω) ≈ 5, which is in line with the value of the second order transport coefficient τπ ≈ 6.1 [55].Due to this rather strong peak in the residue density, the behavior of the Green's function in scalar ϕ 4 theory is actually not to different from the behavior in the conformal relaxation time approximation (see also [38]) where -instead of a spread out region in the complex frequency plane -the non-hydrodynamic contributions originates from a single branch-cut located at Im(ω) = −1/τ R = 0.2.
When increasing k, the hydrodynamic sound mode moves further out into the complex frequency plane and its residue decreases; at the same time the overall contribution of the non-hydrodynamic modes increases, without any dramatic changes in the spectrum of the residue density of the contributing modes.Even though for k = 0.75, the characteristic time scales for the decay of hydrodynamic and non-hydrodynamic modes is comparable, the contribution of the hydrodynamic mode still stands out.With higher k the complex conjugated pair of hydrodynamic modes slowly dives deeper and deeper into the nonanalytic continuum until around k = 1.2, where it gets absorbed into it and is no longer distinguishable.Eventually, for large k the remnant is a large non-analytic region with fairly uniform impact on the Green's function.Nevertheless, one can still see which parts contribute more via the color coding of the Green's function, showing that the largest contributors are located at the edge of the region, towards Re(ω) = ±k.

Conclusions & Outlook
We developed a new method to numerically calculate linear response functions of the energy-momentum tensor in kinetic theory.By carefully discretizing the Boltzmann equation on a momentum grid, this method allows to extract eigenvalues and eigenfunctions of the evolution operator to access the behavior of the Green's function in the complex frequency plane, without the need to perform a numerical Laplace transform.The formalism was tested within the well studied conformal Relaxation Time Approximation and could reproduce the known analytic structure of its Green's functions, including the nonhydrodynamic cut in the lower complex plane.
Subsequently, we explored for the first time the analytic structure of the Green's function of the energy momentum tensor in the sound channel, for a scalar field theory with quartic self-interaction.We find that for small gradients the response of the system is dominated by hydrodynamic modes, which are embedded as single poles within a continuum of non-hydrodynamic excitations corresponding to a non-analytic region in the complex frequency plane that extends arbitrarily close to the real frequency axis, albeit with exponentially suppressed contributions.Nevertheless, for sufficiently small wave numbers, the contribution of non-hydrodynamic excitations to physical observables is strongly peaked around a characteristic value of the imaginary part of the frequency Im(ω) ∼ −τ −1 π .With increasing gradient strength the non-hydrodynamic modes gain more and more influence in the response, until around k ≳ 1.2 the hydrodynamic modes disappear into the continuum and the response of the energy-momentum tensor is entirely determined by non-hydrodynamic excitations.
Based on our analysis, it is thus conceivable that generalized hydrodynamic theories accounting for higher gradient corrections to properly capture the dispersion relations can provide a valid effective description of kinetic theory for sufficiently small gradients k ≲ 1, where Green's functions of the energy-momentum tensor in the sound channel are dominated by the hydrodynamic modes and a subset of non-hydrodynamic excitations with Im(ω) ≈ −1/τ π which could be captured by effective non-hydrodynamic modes.However, in the presence of large gradients k ≳ 1, where a continuum of non-hydrodynamic modes contributes to the Green's function, thus reflecting the underlying microscopic dynamics, it is not clear how hydrodynamics can provide a meaningful description of the dynamics of the system.Since our analysis of the modes is of numerical nature, it is hard to determine the precise location, where the hydrodynamic modes are completely hidden behind the non-hydrodynamic continuum.Moreover, our results indicate that the transition between the two regimes is not particularly sharp, but rather a smooth transition from one to another where hydrodynamic and non-hydrodynamic modes exchange their relative weights.Nevertheless, the location of the threshold kc ∼ 1 is in rather good agreement with the value kc ∼ 0.9 previously obtained in RTA [28,31], which is also reproduced by our numerical RTA results.By converting this estimate into coordinate space, one obtains a critical length scale l c ∼ 1/k c ≈ 0.16 fm 200MeV T η/s 0.16 for typical values of the temperature and transport properties of the QGP.Below this length scale hydrodynamic modes become suppressed and non-hydrodynamic excitations govern the dynamics.The critical length scale is extremely small and should be seen as a lower bound for the applicability of hydrodynamics.
Evidently this study provides a first step to an analogous calculation in QCD kinetic theory, where on general grounds one also expect a non-analytic structure that is far more complicated than just poles and cuts [20,33,34].Beyond the extension to QCD kinetic theory, it would also be interesting to extend the present study in scalar field, from a kinetic description to a genuine QFT treatment, which for weakly coupled theories could be achieved based on n-particle irreducible effective action techniques [64].which can finally be simplified to the form The cosine integral is calculated analytically and only has different values for combinations of i and j, which are tabulated.The p integral is calculated numerically and tabulated for all possible i.

Figure 1 .
Figure 1.Energy Green's function G(ω) in RTA as a function of real frequency ω for two different gradients k.Analytic solution to the Boltzmann equation in black and the solution with our numerical method in red and blue.

Figure 2 .
Figure 2. Eigenvalues of the RTA evolution operator for finite gradients k in the complex frequency plane.Hydrodynamic modes are highlighted in red.

Figure 3 .
Figure 3. Eigenvalues of the RTA evolution operator for various discretizations of N cos θ in the complex frequency plane.The inset plot shows the real time Green's function G( t) as a function of kt for the same discretizations.

Figure 4 .
Figure 4. Top: The residue density ρ µ (ω) for non-conserved observable in scalar ϕ 4 theory for various discretizations as a function of imaginary frequency Im(ω).Bottom: The eigenvalue density ρ(ω) for the same discretizations as a function imaginary frequency Im(ω).

Figure 5 .
Figure 5. Top row and bottom left panel: Eigenvalues of the scalar ϕ 4 evolution operator in the complex frequency plane for different discretizations.Hydrodynamic modes are highlighted as larger points.Bottom right panel: Real time Green's function G( t) as a function of kt for the same discretizations.All plots have been for a gradient of k = 0.8.

Figure 6 .
Figure 6.The real (left) and imaginary (right) part of the dispersion relation ω( k) of the hydrodynamic sound mode as functions of k in scalar ϕ 4 theory (blue) and RTA (magenta) .Black curves show the dispersion relation in second order hydrodynamics in Eq.(4.1) for comparison.

Figure 7 .
Figure 7. From left to right: 1.The real time Green's function G( t) in scalar ϕ 4 theory decomposed into hydrodynamic and non-hydrodynamic contributions as a function of kt .2. The real frequency Green's function G(ω) split into real and imaginary part as a function of the real frequency Re(ω)/ k. 3. The eigenvalues of the evolution operator as black circles in the complex frequency plane.The eigenvalue circle sizes are scaled by their respective residues µ i .The coloring in the plane is the logarithm of the absolute square of the Green's function log |G(ω)| 2 at that position in frequency space.4. The residue density summed out over the real frequency range and plotted as function of Im(ω), the hydrodynamic mode residue is separate as green bar as comparison.The results in each row are obtained for one gradient k.

Figure 8 .
Figure 8. From left to right: 1.The real time Green's function G( t) in scalar ϕ 4 theory decomposed into hydrodynamic and non-hydrodynamic contributions as a function of kt .2. The real frequency Green's function G(ω) split into real and imaginary part as a function of the real frequency Re(ω)/ k. 3. The eigenvalues of the evolution operator as black circles in the complex frequency plane.The eigenvalue circle sizes are scaled by their respective residues µ i .The coloring in the plane is the logarithm of the absolute square of the Green's function log |G(ω)| 2 at that position in frequency space.4. The residue density summed out over the real frequency range and plotted as function of Im(ω), the hydrodynamic mode residue is separate as green bar as comparison.Each row shows the results for a different gradient k.