Spectral Analysis and Hydrodynamic Manifolds for the Linearized Shakhov Model

We perform a complete spectral analysis of the linearized Shakhov model involving two relaxation times $\tau_{\rm fast}$ and $\tau_{\rm slow}$. Our results are based on spectral functions derived from the theory of finite-rank perturbations, which allows us to infer the existence of a critical wave number $k_{\rm crit}$ limiting the number of discrete eigenvalues above the essential spectrum together with the existence of a finite-dimensional slow manifold defining non-local hydrodynamics. We discuss the merging of hydrodynamic modes as well as the existence of second sound and the appearance of ghost modes beneath the essential spectrum in dependence of the Prandtl number.


Introduction
Hydrodynamic closures derived from kinetic theory are a fruitful research direction in statistical physics [20,43] and of basic interest for the investigation of kinetic models [12].The fundamental question: What is the connection between kinetic equations (in the small relaxation regime) and the equations for the motion of continua?Or, to phrase it differently: Can the governing equations of fluid dynamics be rigorously derived from kinetic theory?This problem has a long history.Famously, in his speech at the International Congress of Mathematics's in Paris in 1900, Hilbert proposed a program to derive the passage from the atomistic view of fluids and gases to the motion of continua [29].A modern interpretation of this challenge, known as "Hilbert's sixth problem" in this context, aims to proof the convergence of kinetic models, such as the Boltzmann equation, to hydrodynamic models, such as the Navier-Stokes equation [42,43].
There are several ways to tackle this problem [47].Assuming that the collision term scales as ε −1 , a widely used approach is to expand the density function as a formal power series in ε (the Knudsen number), called Chapman-Enskog series [12].Indeed, the zeroth order PDE obtained from this (singular) Taylor expansion gives the Euler equation, while the first order PDE reproduces the Navier-Stokes equation.We stress, however, that this holds on a formal level only.
While the the solution to the underlying kinetic equation decays to the global equilibrium due to increase of entropy, higher-order terms in the Chapman-Enskog expansion might exhibit instabilities.In [5], it was first shown that an expansion in terms of Knudsen number can lead to nonphysical properties of the hydrodynamic models: At order two (Burnett equation [12]), the dispersion relation shows a change of sign, thus leading to modes which grow in energy (Bobylev instability).Indeed, as pointed out by, e.g., Slemrod [45], convergence of the singular expansion to the leading-order equation is by no means obvious: the formation of shocks might be an obstacle to global uniform convergence in the sense of solutions [46].Furthermore, the expansion of a In our previous papers, similar considerations were already carried out for the three-component Grad system [33], the one-dimensional BGK equation with mass conservation only [34], and recently for the three-dimensional BGK equation with the mass, momentum and energy conservation [36].Similar considerations have been carried out in [9,10] for the one-dimensional linear BGK with one conservation law that of mass in the context of grossly determined solutions (in the sense of [47]).In [2], (optimal) decay rates for various simplified BGK models where derived in the context of hypocoercivity [48].
In comparison with the previously mentioned kinetic models, Shakhov's equation admits more realistic behavior, i.e., mimics properties of the full Boltzmann equation more closely.Indeed, due to the presence of two time-scales (τ fast , τ slow ) in Shakhov's model does not only admit hydrodynamic modes, but non-hydrodynamic (fast) modes as well.In difference with the BGK equation, whose modal branches stay coherent until they mix with the essential spectrum [36], modal branches of the Shakhov equation can mix and produce new branches.For instance, two diffusion modes can collide and produce a second pair of acoustic modes, see Section 4.5.Moreover, Shakhov's model and its versions are widely used in gas-dynamics applications, in particular, in the lattice Boltzmann computations of compressible flows, see, e.g., [18].We emphasize that the techniques outlined in this paper can easily be applied to a wide class of similar kinetic models, such as the ES-BGK [30].
In the following, we will give a complete and (up to a solution of a transcendental equation) explicit description of the spectrum of the Shakhov model linearized around a global Maxwellian.We will show the existence of finitely many discrete eigenvalues above the essential spectrum as well as the existence of a critical wave number for each family of modes.More precisely, we prove the following: {λ N (k|)}, where Modes denotes the set of modes (branches), which might change with wave number and Prandtl number.This is due to a collision of roots and a subsequent bifurcation (for details, we refer to Section 4.5).The essential spectrum is given by the line λ = − 1 τ fast , while the discrete spectrum consists of a finite number of discrete, isolated eigenvalues.Along with each family of modes (up to merging), there exists a critical wave number k crit,N , limiting the range of wave numbers for which λ N exists.For small wave numbers and Pr close to one, the set of modes is given by while the higher wave-number and Pr closer to zero, set of modes is given by where λ shear,1 , λ shear,2 denote the (real) primary and secondary shear modes (double degenerated), {λ ac,1 , λ * ac,1 }, {λ ac,2 , λ * ac,2 } denote pairs of complex conjugated roots, the primary and secondary acoustic modes and the real roots λ diff,1 , λ diff,2 denote the (real) primary and secondary diffusion modes, respectively.
For discussion of different branches of eigenvalues in kinetic models, we refer to [14,17].Our proof is based on the theory of finite-rank perturbations (see, e.g., [51]), together with some properties of the plasma dispersion function, collected in the Appendix for the sake of completeness.Furthermore, we give a hydrodynamic interpretation of the results by considering the dynamics on the slow hydrodynamic manifold (linear combination of eigenspaces).
The paper is structured as follows: In Section 2, we introduce some notation and give some basic definitions.In Section 3, we formulate the fundamental equations (the detailed linearization around a global Maxwellian as well as the non-dimensionalization are included to the Appendix).Section 4 is devoted to the spectral analysis of the linear part, including the derivation of a spectral function describing the discrete spectrum completely.We also give a proof of the finiteness of the hydrodynamic spectrum together with a description of the modes (primary shear, primary diffusion, primary acoustic, secondary shear and secondary diffusion) in frequency space.We also comment on the merging of branches and the formation of a secondary acoustic branch, called second sound for a certain range of Prandtl numbers, see, e.g., [40] or [27] for an analogous phenomenon in solids (phonons).Finally, in Section 5, we write down the hydrodynamic manifold as a linear combination of eigenvectors and derive a closed system for the linear hydrodynamic variables.

Notation and Basic Definitions
For a wave vector k ∈ Z 3 , k = (k 1 , k 2 , k 3 ), we denote its wave number as . Let H denote a Hilbert space and let T : H → H be a linear operator with domain of definition D(H).We denote the spectrum of T as σ(T) and its resolvent set as ρ(T).The spectral analysis of the main operator L of the paper (to be defined later) will be carried out on the Hilbert space (2.2) together with the inner product where the star denotes complex conjugation.Because of the unitary properties of the Fourier expansion, we can slice the space H for each wave number k and analyze the operator L k (restriction of L to the wave number k) on the Hilbert space (2.4) together with the inner product For n = (n 1 , n 2 , n 3 ), the three-dimensional (multidimensional) Hermite polynomials {H n } n∈N 3 are defined via the generating function for the coefficients Let us denote the j-th standard basis vector of R 3 as e j .The three-dimensional Hermite polynomials obey the recurrence relation for j = 1, 2, 3, as well as the orthogonality relation (2.9) 1 (2π) where δ m,n is the Kronecker delta.In particular, the sequence of one-dimensional Hermite polynomials H n (v) obey the recurrence relation which, by (4.15), implies the differential recurrence relation We introduce the plasma dispersion function as the integral (2.12)We collect further useful properties of Z in the Appendix.
Remark 2.1.The plasma dispersion function (2.12) -as its name suggests -appears in plasma physics in the context of Landau damping [16].An detailed description of this widely used nonelementary function is presented in [19].

Preliminaries, Linearization and Non-Dimensionalization
Consider the kinetic model for an unknown distribution function F and the collision operator Here, R denotes the gas constant, while denotes the Prandtl number for the fast time scale τ fast and the slow time scale τ slow .The physical units are given as [R] = m 2 kgs −2 K −1 and [RT ] = m 2 kgs −2 respectively, while the Prandtl number is dimensionless.The number density n, the velocity u, the pressure p and the heat flux q are defined by while the temperature T is defined through the relation .
We stress that the macroscopic quantities (n, u, T, q) depend on F through the functional relationship (3.4), but are independent of the velocity v. Equation (3.1) is called Shakhov's S-model and was first derived in [44] as a generalization to the BGK equation, allowing for two time scales, which, in particular, allows to define a family of models with varying Prandtl number.To ease notation in the following calculations, we set Remark 3.1.The Prandtl number (and the parameter r) allow us to define an interpolation between the three-dimensional BGK equation (Pr = 1 r r = 0), see e.g.[49], and a model with maximal separation between fast and slow time-scale (Pr = 0 or r = 1).The different properties of the spectra, see Section 4, vary for different values of r.
For n ≥ 0, we define the linear moments for n ≥ 2 is the n th tensor power.For compactness in the presentation, we also define the special vector The hydrodynamic variables (n, u, T, q) are related to the moments (3.8) via which can be inverted to (3.12) Equation (3.1) admits a global Maxwellian, (3.13) as an equilibrium solution.Here, the equilibrium number density n 0 and the equilibrium temperature T 0 are constants.In the following, we will be interested in the dynamics of (3.1) close to the stationary solution (3.13), i.e., the linearized dynamics of (3.1) around (3.13).The Shakhov model linearized around (3.13) in non-dimensional form is given by For the details, we refer to Appendix I. Equation (3.14) will serve as the basis for the spectral analysis performed in the following section.

Spectral Analysis of the Linearized Two-Timescale Operator
In this section, we will carry out a complete spectral analysis of the right-hand side of (3.14), following the approach in [36].This will allow us to draw conclusions on the decay properties of hydrodynamic variables, the existence of a critical wave number and the hydrodynamic closure.After reformulating the problem in frequency space, we will use the resolvent calculus to formulate a condition for the discrete spectrum (Subsection 4.1).Then, we will use properties of the plasma dispersion function (see Appendix) to define a spectral function Σ |k|,τ , whose zeros coincide with the discrete, isolated eigenvalues (Subsection 4.2).These families of eigenvalues are described in more detail in Subsection 4.3.Then, in Subsection 4.4, we prove the existence of a critical wave number k crit such that Σ |k|,τ has no zeros (i.e., there exists no eigenvalues) for |k|> k crit .In Subsection 4.5, we take a closer look at the branches of eigenvalues (modes) and the merging of diffusive branches (second sound).Finally, in Subsection 4.6, we discuss the existence of eigenvalues below the essential spectrum (ghost modes) for r < 0 (Pr > 1).4.1.Description of the discrete spectrum.In the following, we rescale the density f with a global, non-dimensional Maxwellian, (4.1) f → (2π) −3/2 e − |v| 2 2 f, which allows us to divide by the Gaussian in (3.14) and interpret the moments (A.4) as projections relative to the inner product (2.5).We define the following set of basis functions which satisfy the orthonormality condition (4.6) e n , e m v = δ nm , for 0 ≤ n, m ≤ 7.
Defining (4.7) we can infer the following relations between the moments and the coefficients (4.7):We bundle the basis functions (4.2) into a vector (4.9) e = {e j } 0≤j≤7 , and define the matrix where H n is defined in (2.6) and f n is an n-tensor.Since the eight basis vectors (4.2) appear in the expansion (4.15) via an orthogonal splitting, we have that (4.16) for any f ∈ H v .Hermite expansions were famously used by Grad in his seminal paper [25] to establish finite-moment closures. From where we have assumed that f is sufficiently regular to justify the application of the divergence theorem in x in order to remove the gradient term as well as (4.16).For 0 < r < 1 (or 0 < Pr < 1). it follows that the operator L is dissipative and that (4.18) σ(L) ≤ 0.
On the other hand, for r > 0, since P 5 and P 8 are orthogonal projections, it follows that This shows that any solution to (4.14) has to converge to zero, i.e., the global Maxwellian is a stable equilibrium up to the conserved quantities from the projected modes.On the other hand, we infer that the overall convergence rate to equilibrium can be at most − 1 τ for r > 0. For r < 0, we can estimate where we have used that P 8 f 2 ≤ f 2 as well.Consequently, for r < 0, we can only infer the weaker decay rate r−1 τ as opposed to 1 τ for r > 0. Remark 4.2.The weaker decay rate in (4.20) for negative r, i.e., Prandtl number larger than one, is already indicative for the existence of parts of the spectrum located below the essential spectrum ({ λ = − 1 τ }).Indeed, we will show in Section 4.6 that for r < 0, there exist eigenvalues below the essential spectrum for a certain range of wave numbers.
Let us proceed with the spectral analysis by switching to frequency space.Since x ∈ T 3 , we can decompose f in a Fourier series as for the Fourier coefficients In frequency space, (4.13) becomes and the spectrum of L can be calculated from the corresponding operator at each wave vector k: For k = 0, we can read off the spectrum of (4.23 for 5 ≤ l ≤ 7. On the other hand, we see from (4.25) that L0 just acts like − 1 τ on the orthogonal complement of span{e j } 0≤j≤7 .This shows that and the dimensions of the corresponding eigenspaces are given by Now, for the following, let us assume that k = 0.For more compact calculation, we define the operator (4.29) In the calculation of the discrete spectrum of the operator Lk , based on the second resolvent identity and finite-rank perturbations, we also follow the presentation in [36] closely.The spectrum of Lk is then given by where we have used that σ(S k ) = R for k = 0. We define the Green's function matrices as By the second resolvent identity, for any operators A, B and z ∈ ρ(A) ∩ ρ(B), we have for Applying equation (4.34) to e m for 0 ≤ m ≤ 4 and rearranging gives for z ∈ C \ iR.Thus, the resolvent of (iτ S k − B 8,r − z) includes the resolvent of iτ S k as well as information from the matrix {G T (z, n, m)} 0≤n,m≤7 as coefficients.
Taking an inner product of (4.35) with e m gives for 0 ≤ n, m ≤ 7 and z ∈ C \ iR.System (4.36) defines sixty-four equations for the coefficients G T (z, n, m), which can be re-written more compactly as (4.37) or, equivalently, Consequently, we can solve for the entries of G T unless det(Id − D r G S ) = 0, which implies that − Id = 0 .
An eigenvalue λ of the operator Lk is then related to the zero z in (4.39) via To ease notation, we define the spectral function 4.2.Derivation of a Spectral Function.To evaluate the integral expression in (4.41), we we decompose the wave vector k into its magnitude along a coordinate direction and a rotation: where while the vector of basis functions e transforms according to . By the orthogonality of Q k , the volume element transforms as dv = dw and we can calculate det where we have used that D r only acts on the last three columns by multiplication with a constant and hence commutes with the block-rotation-matrix, together with the orthogonality of Q k in factoring the determinant.
From equation (4.46), we see that the spectral function Σ k,τ depends on the wave vector k only through τ k and we define (4.47) κ := τ k.
A lengthy but elementary calculations allows us to integrate out the variables w 2 and w 3 in (4.46), which simplifies the spectral function according to det where , and the coefficient matrices are given by M n w n , for some matrix coefficients M n ∈ R 8×8 .
To evaluate the integral expression on the right-hand side of (4.48), we will use the plasma dispersion function (2.12).To this end, we have to calculate expressions for the derivative of Z.By repeated application of (2.13), we conclude that (4.52) for some polynomials p n , q n with integer coefficients.The first few derivatives of Z can be expressed as dZ dζ = −1 − ζZ, (4.53) We claim that q n (ζ) = (−1) n H n (ζ) for the n th Hermite polynomial.Indeed, from the recurrence relation of Z in (2.13), we have that showing that q n+1 = q n −ζq n , which is -up to a sign flip -the recurrence relation of the n th Hermite polynomial (2.11).With the relation we can now conclude that for any polynomial expanded in Hermite basis P (w) =  With the help of the plasma dispersion function Z and the insight (4.56), we readily calculate , and the coefficient matrices are given by Using the change of coordinates for the polynomials Let us conclude this section with some remarks.The main result of the preceding calculationsthe derivation of the spectral function (4.62) with coefficient polynomials and factorization (4.64) -allows us to conclude specific features of the spectrum by solving for the zeros of a holomorphic function.This is a tremendous simplification compared to the study of (4.13) directly, where the transport and the collision term interact in a delicate manner.This is as close an analog of a determinant in finite-dimensional systems as one could hope for.In the following section, we will take a closer look at the structure of the zero set of (4.62).In particular, we will identify different families of zeros (branches) that relate to hydrodynamics.4.3.Hydrodynamic Modes.In this section, we will identify the branches of the zero set of the spectral function Σ k,τ in dependence on the modified wave number κ and the parameter r.First, let us show that in the limit κ → 0 (for fixed τ ), we recover the spectral structure of (4.27).To this end, we use the asymptotic expansion under the assumption that ζ > 0. For a proof of (4.65), we refer to the Appendix.The limit κ → 0 with λ < 0 satisfies the assumptions of (4.65).A lengthy calculation shows that, indeed, , for κ → 0 and λ > 0, (4.66) which is consistent with (4.27).Since, for small κ, the spectral function Σ k,τ is continuous in κ, it follows that there exist five branches of eigenvalues (indexed by κ) emerging out of the five-fold degenerate zero λ = 0, which we call primary hydrodynamic modes.Furthermore, there exist three branches of eigenvalues emerging from the zero λ = − 1 τ slow which we call secondary hydrodynamic modes.From the factorization in (4.62) and (4.64), we see that any zero of Σ 1 will occur with algebraic multiplicity of two.For small wave number, the primary hydrodynamic modes consist of a pair of complex conjugate eigenvalues, the primary acoustic modes, a real eigenvalue of algebraic multiplicity two (but geometric multiplicity one), called primary shear mode and a algebraically and geometrically simple, real eigenvalue, called primary diffusion mode.Similarly, for small wave number, the secondary hydrodynamic modes consist of a real eigenvalue of algebraic multiplicity two (but geometric multiplicity one), called secondary shear mode and a algebraically and geometrically simple, real eigenvalue, called secondary diffusion mode, see Figure 4.1.
For larger wave numbers, depending on P r, the primary and secondary shear modes may collide and produce another pair of complex conjugated eigenvalues, called secondary acoustic modes or second sound, see Figure 4.1.This occurs through a saddle-node bifurcation, see Section 4.5.
We denote the families of modes index by wave number k as where the low wave-number set of modes is given by (4.68) while the higher wave-number set of modes is given by (4.69) Modes 2 = {shear 1 , ac 1 , ac 1 * , shear 2 , ac 2 , ac 2 * }.
As mentioned before, in a certain range of Prandtl numbers, the set of modes can change from (4.68) to (4.69) if the wave number is increased.

4.4.
Existence of a Critical Wave Number and Finiteness of the Hydrodynamic Spectrum.Along the same lines as in [36], we will show that for each family of modes, there exists a critical wave number such that k crit such that (4.70) In fact, each mode has its own critical wave number, which depends on the specific properties of the branch.s Proof.The claim follows from a combination of Rouché's theorem applied to the spectral function Σ k,τ with the asymptotic expansion (4.65).Indeed, for fixed κ, we find that

19).
At small wave numbers (k = 0.4 and k = 0.5), the primary diffusion mode decreases along the real axis, while the secondary diffusion mode increases along the real axis.Around k ≈ 0.6 they collide, a bifurcation takes place and a pair of complex conjugated modes, the secondary acoustics, is created (k = 0.7).

|arg(ζ)|≤ π
2 − δ, ζ → ∞, for any real number 0 < δ ≤ π 2 .Analogously, we find that  L mf p = τ v thermal , and transforming back to physical units, we see that the critical wave number is numerically proportional to the inverse of (4.74).Indeed, we obtain that (4.75) Remark 4.5.In fact, each family of modes has its own critical wave number k crit,N .Since branches can merge, 4.5.Global Bifurcation of Eigenvalues, Merging of Branches and Second Sound.In this section, we discuss the phenomenon of branch merging already indicated by (4.68) and (4.69) in more details.Through this section, we assume 0 ≤ r ≤ 1 (for a discussion of r < 0 and the existence of ghost modes, we refer to the following section).
In general, any zero of the spectral function can be expanded in a Puiseux-Newton series with appropriately chosen exponent as already noted in, e.g., [37].Namely, the shear modes defined by Σ1 can be expanded in a Puiseux-Newton series with exponent 1 4 (four-fold degeneracy for k = 0 and r < 1).The other modes defined by Σ2 can be expanded in a Puiseux-Newton series with exponent 1 4 as well (also a four-fold degeneracy for k = 0 and r < 1).A lengthy (and cumbersome) expansion calculation shows, however, that each branch is actually analytic in k and we can therefore expand where λ 0 and λ 1 determine the branch of eigenvalues.This is consistent with the asymptotic expansions for the hydrodynamic branches derived in [14] for general kinetic equations and small wave number.Indeed, the instantaneous directional motion of an eigenvalue k → λ(k) is given by (4.77) To obtain the coefficients λ n (r), we plug (4.76) into the spectral function (4.41) and compare powers of k (using the asymptotic expansion (B.16)): Then, we can solve the equations G n for λ n successively.A lengthy calculation shows that (4.79) which is consistent with (4.27).To ease notation, we define the two families of branches as Expanding further reveals that (4.81) Furthermore, we find that which implies that the secondary diffusion mode moves to the right for 56 81 < r < 1, while it moves to the left for 0 < r < 56 81 initially.This shows that, for r sufficiently close to one, the primary and secondary diffusion branch will inevitably collide and produce a pair of complex-conjugate zeros (secondary acoustics) through a saddle-node bifurcation.On the other hand, we see that for r close to 0, the secondary diffusion mode will travel to the left from the very beginning and will leaf the domain − 1 τ < λ < 0 before it gets the chance to collide with the primary diffusion branch, see Figure 4.3.Consequently, there does not exist the phenomenon of branch merging and second sound.Indeed, in the limit r → 0, the Shakhov S-model reduces to the BGK model, for which no second sound exists.We summarize that -if the Prandtl number is close to one, i.e., close to the dynamics of the three-dimensional linear BGK equation, The instantaneous second-order directional motion at k = 0 of the secondary diffusion mode in dependence of the parameter r.For 56 81 < r < 1, the secondary diffusion mode moves towards zero until it collides with the primary diffusion mode (second sound).For 0 < r < 56 81 , the diffusion mode moves towards the essential spectrum at k = 0 -even with a singularity at r = 0 (BGK equation).For r < 0, the secondary diffusion mode is a ghost mode (below the essential spectrum) and always moves towards it instantaneously.there exists no second sound (clearly, there is no second sound for the BGK equation), while for Prandtl number close to zero, there will always be branch merging.
Remark 4.6.While the coefficient (4.83) governs the motion of the secondary diffusion mode for small wave numbers, the behavior for larger k might be different.Indeed, for a certain value of r, the secondary diffusion mode might start out by moving to the left, but then turn to the right and collide with primary diffusion mode anyways.The above considerations thus only imply that there is a certain range for which second sound exists and that there is a certain range for which it does not exist.4.6.Spectral Properties for Prandtl number greater than one: existence of ghost modes.In this section, we derive the behavior of the spectrum of (3.14) for r < 0 (or P r > 1).As already indicated by the estimate (4.20), potential eigenvalues are no longer guaranteed to exist above the essentially spectrum exclusively for r < 0. Figure (4.4) shows some typical argument plots of the spectral function for P r = 1.5 (r = −0.5)with eigenvalues below the essential spectrum.Obviously, from the considerations around the spectrum of L0 in (4.27), we see that for P r > 1 (which is equivalent to τ fast > τ slow ), the eigenvalue − 1 τ slow is indeed located below the essential spectrum  define eigenvalues of the linearized Shakhov mode (points where a small, counterclockwise loop runs through the whole rainbow at least once).All eigenvalues have negative real part and are located above the essential spectrum { λ = − 1 τ } (solid black line), which is consistent with the decay estimates (4.19).Already at small wave numbers (k = 0.4 and k = 0.5), the secondary diffusion and secondary shear mode are close to the essential spectrum ({ λ = −2}).For this Prandtl number, the secondary diffusion mode is smaller than the secondary shear mode.At larger wave numbers (k = 0.5 and k = 0.6), these modes decrease.At k = 0.7, the secondary diffusion mode has already disappeared and no merging of branches can occur at this Prandtl number.
From the instantaneous directional motion of the secondary diffusion mode (4.83), we see that which shows that the ghost (diffusion) mode increases with wave number until it is absorbed by the essential spectrum, see Figure 4.4.A similar argument holds true for the secondary shear modes.Existence of non-hydrodynamic modes below the essential spectrum suggests that Shakhov's model becomes unrealistic Pr > 1, and other kinetic models, such as the ES-BGK should be used instead for Pr > 1.

Linear Hydrodynamic Manifolds
We define a hydrodynamic manifold as the (unique) slow spectral submanifold associated to a set of eigenvectors.In the case of linear dynamics, the manifold itself is given by the invariant linear subspace spanned by the eigenvalues associated to a set of slow eigenvalues.In particular, the hydrodynamic manifold has the following properties: (1) It contains an appropriately scaled, spatially independent stationary distribution (e.g.global Maxwellian) as a base solution (2) The projection onto the hydrodynamic moments along the manifold provide a closure of the hydrodynamic moments (mass-density, velocity and temperature) (3) It attracts all trajectories in the space of probability-density functions (which are close enough to the base solution) exponentially fast, thus acting as a slow manifold (4) It is unique.We write the hydrodynamic manifold in frequency space as a superposition of eigen-functions associated to each mode   All eigenvalues have negative real part, but some of them are located below the essential spectrum { λ = − 1 τ } (solid black line), which is consistent with the decay estimates (4.20).At small wave numbers, three ghost modes, namely a shear mode and a diffusion mode, appear below the essential spectrum (k = 0.4).As the wave number is increased, the ghost modes move closer towards the essential spectrum until the diffusion mode is absorbed (k = 0.5, k = 0.8).Finally, also the ghost shear mode is absorbed (k = 1.2) and the spectrum consist of primary hydrodynamic modes above the essential spectrum only (up to the critical wave number).
The vector β describes the evolution of the hydrodynamic variables in terms of the basis of eigenvalues (spectral basis).The time-evolution of the hydrodynamics can then be written as or, solving for β in (5.8): (5.12) ĥt = AΛA −1 ĥ.
Equation (5.12) defines the exact hydrodynamics derived from the slow motion along the hydrodynamic manifold.
Remark 5.1.The evaluation of the right-hand side of (5.12) is by no means trivial and involves properties of the spectral projection.The properties of the exact, spectrally-closed hydrodynamics together with their physical properties will be discussed in a forthcoming paper [35].

Conclusion and Further Perspectives
We performed a complete spectral analysis for the Shakhov model linearized around a global Maxwellian.The discrete eigenvalues above the essential spectrum { λ = − 1 τ fast } are described as zeros of a spectral function at each wave number.In this way, we identified families of modes (branches), depending on wave number.For small wave numbers, the family of modes is given by Modes = {ac 1 , ac 1 * , diff 1 , diff 2 , shear 1 , shear 2 }, the pair of primary acoustic modes, the primary and secondary diffusion modes as well as the primary and secondary shear modes.Within a certain range of Prandtl numbers, a merging of branches can occur at a specific wave number and the modes diff 1 and diff 2 may collide, producing another pair of acoustic modes ac 2 and ac 2 * via a saddle-node bifurcation.This phenomenon is known as second sound.The approach presented in [34,36] as well as in the current paper is general enough to infer spectral properties for any finitely-truncated collision operator, such as quasi-equilibrium approximations [24] or even Maxwell molecules [47].The explicit knowledge and quantitative properties of the spectra identified for several kinetic model equations [34,36] also allow us to move to the existence theory of non-linear hydrodynamic equations for various (finitely-truncated) kinetic models.Indeed, the fact that the discrete spectrum is well separated from the essential spectrum allows us to define a spectral projection for the whole set of eigenvalues, thus giving the first-order approximation (in terms of nonlinear deformations) to the hydrodynamic manifolds.In particular, we expect that the theory of thermodynamic projectors [22] may be helpful in proving the nonlinear extension.The quantitative insights in the structure of the spectrum could also be used to derive simplified, but still non-local, approximate hydrodynamics.This could also improve present numerical methods [32].

Appendix A. Linearization and Non-Dimensionalization of the Shakhov Equation Around a Global Maxwellian
In this section we perform -for the sake of completeness -the linearization of the Shakhov model around a global Maxwellian.To obtain the linearization of (3.1) around F eq 0 , we write F = F eq 0 +εf and calculate Using the relations n[F eq 0 ] = n 0 , u[F eq 0 ] = 0, q[F eq 0 ] = 0, p[F eq 0 ] = mRn 0 T 0 , T [F eq 0 ] = T 0 , (A.2) which in particular imply that F eq [F eq 0 ] = F eq 0 , we can reduce (A.1) to the moments of the perturbation density f , the moments (3.8) transform according to which in turn implies that Consequently, the ε-derivatives of the hydrodynamic moments become (A.9) Combining (A.8) with (A.9), equation (A.3) reads Function (B.1) is called Faddeeva function and is frequently encountered in problems related to kinetic equations [16].We then have that and, by relation (B.2), we have for ζ < 0: Consequently, we obtain as can be seen from (B.6).Function (B.6) satisfies an ordinary differential equation (in the sense of complex analytic functions) on the upper and on the lower half-plane.Indeed, integrating (2.12) by parts gives for any 0 < δ ≤ π 2 , see also [31].The proof will be based on a generalized version of Watson's Lemma [50].To this end, let us define the Laplace transform (B.17 for any real number 0 < δ ≤ π 2 , where Γ is the standard Gamma function.For a proof of the above Lemma, we refer e.g. to [15].Classically, Lemma (B.1) is applied to prove that the imaginary error function admits an asymptotic expansion for x ∈ R of the form (B. 20) erfi(x) ∼ e x 2 √ πx ∞ k=0 (2k − 1)! !(2x 2 ) k , for x > 0, x → ∞, see also [39], based on the classical version of Watson's Lemma, whose assumptions are, however, unnecessarily restrictive [52].
For completeness, we recall the derivation of (B.

,f = 4 n=0f, e n e n , P 8 f = 8 n=5f
together with the following (4.11)B 8,r = (P 5 + rP 8 )f = f, e v • D r e, as the scaled sum of two finite-rank projections (4.12) P 5 , e n e n .The linear operator appearing as the right-hand side of equation (3.14) (together with the rescaling (4.1)) then takes the simple form (4.13)

(4. 45 )
Remark 4.3.The first column of the rotation matrix Q k can be chosen as1  k k, which, by the orthonormality of the columns of Q k , implies(4.44).The change of coordinates v = Q k w can then interpreted as writing the velocity vector v as the sum of a a component parallel to the wave vector k and a component orthogonal to it:

50 )
Since the entries of M are polynomials of degree at most six, we can expand (4.51) M(w) = 6 n=0

3 D
57), we can evaluate the integral expression in(4.46)  to an explicit matrix depending on polynomials in ζ and Z(ζ) in a linear way.Indeed, with the help of (4.57), we have that det R D r N(ζ) − (iκ)Id) , 62) shows that the spectral function factors into two parts, which we denote as

Figure 4 . 1 .
Figure 4.1.Argument plot of the spectral function (4.41) for relaxation time τ = 0.5, Prandtl number P r = 0.4 and different wave numbers k.The zeros of (4.41) define eigenvalues of the linearized Shakhov mode (points where a small, counterclockwise loop runs through the whole rainbow at least once).All eigenvalues have negative real part and are located above the essential spectrum { λ = − 1 τ } (solid black line), which is consistent with the decay estimates (4.19).At small wave numbers (k = 0.4 and k = 0.5), the primary diffusion mode decreases along the real axis, while the secondary diffusion mode increases along the real axis.Around k ≈ 0.6 they collide, a bifurcation takes place and a pair of complex conjugated modes, the secondary acoustics, is created (k = 0.7).

Figure 4 .
Figure 4.2.The instantaneous second-order directional motion at k = 0 of the secondary diffusion mode in dependence of the parameter r.For 56 81 < r < 1, the secondary diffusion mode moves towards zero until it collides with the primary diffusion mode (second sound).For 0 < r < 56 81 , the diffusion mode moves towards the essential spectrum at k = 0 -even with a singularity at r = 0 (BGK equation).For r < 0, the secondary diffusion mode is a ghost mode (below the essential spectrum) and always moves towards it instantaneously.

Figure 4 . 3 .
Figure 4.3.Argument plot of the spectral function (4.41) for relaxation time τ = 0.5, Prandtl number P r = 0.6 and different wave numbers k.The zeros of (4.41) define eigenvalues of the linearized Shakhov mode (points where a small, counterclockwise loop runs through the whole rainbow at least once).All eigenvalues have negative real part and are located above the essential spectrum { λ = − 1 τ } (solid black line), which is consistent with the decay estimates (4.19).Already at small wave numbers (k = 0.4 and k = 0.5), the secondary diffusion and secondary shear mode are close to the essential spectrum ({ λ = −2}).For this Prandtl number, the secondary diffusion mode is smaller than the secondary shear mode.At larger wave numbers (k = 0.5 and k = 0.6), these modes decrease.At k = 0.7, the secondary diffusion mode has already disappeared and no merging of branches can occur at this Prandtl number.

Figure 4 . 4 .
Figure 4.4.Argument plot of the spectral function (4.41) for relaxation time τ = 0.5, Prandtl number P r = 1.5 (r = −0.5)and different wave numbers k.The zeros of (4.41) define eigenvalues of the linearized Shakhov mode (points where a small, counter-clockwise loop runs through the whole rainbow at least once).All eigenvalues have negative real part, but some of them are located below the essential spectrum { λ = − 1 τ } (solid black line), which is consistent with the decay estimates (4.20).At small wave numbers, three ghost modes, namely a shear mode and a diffusion mode, appear below the essential spectrum (k = 0.4).As the wave number is increased, the ghost modes move closer towards the essential spectrum until the diffusion mode is absorbed (k = 0.5, k = 0.8).Finally, also the ghost shear mode is absorbed (k = 1.2) and the spectrum consist of primary hydrodynamic modes above the essential spectrum only (up to the critical wave number).