Shakhov-type extension of the relaxation time approximation in relativistic kinetic theory and second-order fluid dynamics

We present a relativistic Shakhov-type generalization of the Anderson-Witting relaxation time model for the Boltzmann collision integral to modify the ratio of momentum di ff usivity to thermal di ff usivity. This is achieved by modifying the path on which the single particle distribution function f k approaches local equilibrium f 0 k by constructing an intermediate Shakhov-type distribution f S k similar to the 14-moment approximation of Israel and Stewart. We illustrate the e ff ectiveness of this model in case of the Bjorken expansion of an ideal gas of massive particles and the damping of longitudinal waves through an ultrarelativistic ideal gas.


Introduction
The relativistic Boltzmann equation, k µ ∂ µ f k = C[ f ], determines the space-time evolution of a local invariant distribution function f k = f (k µ , x µ ), where C[ f ] is the collision term specifying the interaction among constituents.The relaxation time approximation (RTA), introduced by Anderson and Witting (AW) [1,2] as a proper relativistic generalization of the Bhatnagar-Gross-Krook (BGK) model [3], approximates the main features of C[ f ] for both analytic and numeric computations [4,5,6,7,8,9,10,11,12,13,14,15,16].In the AW model the approach to local equilibrium is given by a momentum independent relaxation time, τ R , while the conservation of a given quantum charge, and of energy and momentum is strictly valid choosing the local rest frame and matching conditions of Landau [17].Such relaxation-time models share the common caveat that all first-order transport coefficients are related to τ R .In the BGK model, the Prandtl number Pr, representing a ratio of viscosity to thermal diffusivity, is fixed at 1, while most ideal gases have Pr = 2/3.This limitation was remedied in the extension proposed by Shakhov [18,19] by introducing another parameter that allows Pr to be controlled independently of τ R .This simple modification already leads to a remarkably good agreement with the solutions of the Boltzmann equation and with experimental data, while the BGK model shows significant deviations from both [20,21,22,23].
In this paper, we introduce the relativistic Shakhov model as an extension of the AW model, offering additional energyindependent relaxation times that are related to the corresponding first-order transport coefficients, the bulk viscosity ζ, particle diffusion κ and shear viscosity η.This modified collision term, significantly extends the applicability of kinetic RTA Email address: victor.ambrus@e-uvt.ro(Victor E. Ambrus , ) models and allows to study different parameterizations of the first-order transport coefficients in various physical phenomena such as the baryon diffusion [24], and the study of bulk and shear viscosities [25,26] in high-energy heavy-ion collisions.
The improvements compared to the AW model will be shown by comparing the numerical solutions of the kinetic model with solutions of Israel-Stewart-type second-order fluid dynamics in two different benchmarks: the (0+1)-dimensional longitudinally boost-invariant flow of massive particles highlighting the ζ/η ratio, and the damping of sound waves for an ultrarelativistic fluid concerning the κ/η ratio.
We note that an energy-dependent relaxation time related to the microscopic details of the collisional processes [27,28,29] will also modify the first-order transport coefficients [30,31,32,33,34].For example, setting τ k in the AW model provides an extra parameter ℓ that can be used to tune either ζ/η or κ/η, while it is not excluded that more elaborate models could tune both of these ratios simultaneously.The Shakhov model introduced in this paper provides an alternative to such AW models with energy-dependent relaxation times.
The detailed derivation of second-order fluid dynamics and the entropy production from the Shakhov model as well as the numerical methods are relegated to the supplementary material.

The Anderson-Witting model and transport coefficients
The Boltzmann equation in the Anderson-Witting relaxationtime approximation [1,2] for the binary collision integral reads where E k = k µ u µ is the comoving energy of a particle with four-momentum k µ and rest mass is the local momentum-independent relaxation time.The fluid four-velocity Here, δ f k ≡ f k − f 0k denotes the deviation of f k from its local equilibrium form, the Jüttner distribution [35], where α = µβ and β = 1/T is the inverse temperature, µ is the chemical potential, while a = ±1 for fermions/bosons and a = 0 for classical Boltzmann particles.The fluid four-velocity defining the local rest frame is chosen as u µ = T µν u ν /(u α T αβ u β ), where T µν = dKk µ k ν f k is the energy-momentum tensor, with dK = [g/(2π) 3 ]d 3 k/k 0 denoting the Lorentz-invariant integration measure in momentum space and g is the degeneracy factor.With respect to u µ , T µν and the particle four-flow N µ = dKk µ f k are decomposed as where ∆ µν = g µν − u µ u ν is the projector orthogonal to u µ and g µν = diag(+, −, −, −) is the metric tensor.The bulk viscous pressure Π, the particle diffusion current V µ , and the shearstress tensor π µν represent dissipative corrections due to δ f k .N µ 0 and T µν 0 correspond to the local equilibrium state defined as (4) Here, n ≡ N µ u µ = N µ 0 u µ is the particle density, e ≡ T µν u µ u ν = T µν 0 u µ u ν is the energy density, while P ≡ −T µν 0 ∆ µν /3 = P(e, n) is the thermodynamic pressure defined by an equation of state.
To evaluate the first-order transport coefficients ζ, κ and η, we apply the Chapman-Enskog method [36], hence from Eq. (1), δ f k is considered to be of first order with respect to τ R , The irreducible moments of δ f k are defined as [37]: where r denotes the power of energy are the irreducible tensors forming an orthogonal basis [37,38].The symmetric and traceless projection tensors of rank 2ℓ, are orthogonal to u µ since they are constructed using the ∆ µν projector operators.

The relativistic Shakhov-type model
The main idea behind the Shakhov-type model [18,19] is to replace the collision term (1) with where f Sk = f 0k + δ f Sk replaces f 0k in the AW model.This collision term drives f k towards f Sk , on a time-scale given by τ R , and ultimately towards f 0k on a modified path compared to C AW [ f ] from Eq. (1).We will construct f Sk such that δ f Sk vanishes equivalently in local equilibrium when δ f k = 0.As in the AW model, we also consider that the relaxation time is independent of the four-momentum of particles.Note that an energy-dependent generalization of the RTA was studied in Refs.[30,31,32,33,34].
Applying the Chapman-Enskog method while considering both δ f Sk and δ f k of first-order with respect to τ R , Eq. ( 5) is modified to Alike to Eq. ( 6), we define the irreducible moments of δ f Sk , The first-order relations corresponding to Eq. ( 13) now read Note that Π S ≡ −m 2 0 ρ S ,0 /3, V µ S ≡ ρ µ S ,0 , and π µν S ≡ ρ µν S ,0 are so far unspecified.We aim to use these new irreducible moments to modify the transport coefficients to the following form: where τ Π , τ V and τ π are the new τ R -independent relaxation times of bulk viscosity, particle diffusion and shear viscosity, representing parameters of the Shakhov model.In other words, we seek to replace Eqs. ( 7) by and hence we solve Eqs. ( 15) and ( 17) by setting [18,19] The conservation equations ∂ µ N µ = 0 and ∂ µ T µν = 0 are fulfilled when where N µ S = dKk µ f Sk and T µν S = dKk µ k ν f Sk .Imposing Landau's matching conditions n = n S and e = e S defines the local chemical potential and temperature.The local rest frame and u µ is chosen according to Landau's definition and leads to ∆ µ ν u λ T νλ S = 0. Note however that in the Shakhov model other definitions for u µ are possible, as pointed out using a similar construction in Ref. [40].
The conditions from Eqs. ( 18) and ( 19) can be written in terms of the irreducible moments while all other irreducible moments of δ f Sk are unconstrained.
In order to satisfy Eqs.(20), we construct the relativistic Shakhov distribution f Sk similarly to the near-equilibrium expansion of Grad [41] and Israel-Stewart [42].Starting from we define the Shakhov term S k as: The polynomials H (ℓ) k0 in energy E k are constructed to fulfill Eqs.(20) and were introduced in Ref. [37].In the 14-moment approximation, these lead to:

Second-order fluid dynamics
The Chapman-Enskog method defines the dissipative fields relating the thermodynamic forces through first-order transport coefficients, resulting in the Navier-Stokes relations, Eqs.(7) and (17).The method of moments by Grad [41], and Israel and Stewart [42], based on an expansion around equilibrium, alike the Shakhov distribution, leads to relaxation-type equations of motion for the dissipative fields Π, V µ and π µν .
Using the well known 14-moment approximation of dissipative fluid dynamics [37,42,43,44], the second-order relaxation equations provide closure for the conservation laws.Such relaxation equations can be derived for the Shakhov model, where the additional microscopic time scales introduced through δ f Sk correspond to the relaxation times τ Π , τ V , and τ π from a linearized collision integral, e.g. in the case of binary hard-sphere interactions [37,45].Consequently, some second-order transport coefficients acquire different values than in the case of the AW model, as shown in the Supplementary Material [46].
The equations of motion for the irreducible moments ρ µ   35)-( 37) of Ref. [37]: The irreducible moments of the Shakhov collision term C S [ f ] from Eq. ( 12) are where the first term on the right hand side corresponds to the AW collision term, while the second term involves the irreducible moments of δ f Sk .
The collision matrix A (ℓ) rn is defined by [37,45] where the summation over n is unrestricted from below or from above.Using Eq. ( 25), the collision matrix evaluates to where according to Eq. (66) of Ref. [37] is such that F (ℓ) −r,n = δ rn for all 0 ≤ r, n ≤ N ℓ , with H (ℓ) kn given in Eqs. ( 23) and (N 0 , N 1 , N 2 ) = (2, 1, 0).The different relaxation times are denoted by τ (ℓ) S , with τ (0) S = τ Π , τ (1) S = τ V and τ (2) S = τ π .The inverse of the collision matrix is the relaxation-time matrix, satisfying m τ (ℓ)  rm A (ℓ) mn = δ rn , and it is given by where the elements corresponding to the matching conditions and the choice of frame, ρ 1 = ρ 2 = 0 and ρ µ 1 = 0, are excluded.Multiplying Eqs. ( 24) by the inverse collision matrix τ (ℓ)  nr and summing over r leads to ζ n θ + higher-order terms, (30a) nr ρ⟨µ⟩ r + ρ µ n = κ n ∇ µ α + higher-order terms, (30b) r τ (2)   nr ρ⟨µν⟩ r + ρ µν n = 2η n σ µν + higher-order terms, where the first-order transport coefficients are defined as Now using Eq. ( 29) in the above formulas, we get thus the first-order transport coefficients corresponding to r = 0 reduce to in agreement with Eqs. ( 8) and (16).Contrasting these results with the transport coefficients of the AW model, the Shakhov model allows a more general approach analogous to the 14-moment approximation of second-order fluid dynamics.Note that, aside from the first-order transport coefficients in Eqs. ( 16), the Shakhov term in Eq. ( 22) also affects some second-order transport coefficients.These are discussed in the Supplementary Material [46].We will illustrate the capabilities of the relativistic Shakhov model by comparing it to the Anderson-Witting model and to second-order fluid dynamics in two different settings.
The transverse P ⊥ and longitudinal P l pressure components are related to the pressure P, the bulk pressure Π, and the shearstress tensor component π = π η s η s via In the Shakhov model, the time evolution of P, Π and π is obtained by directly solving the Boltzmann equation as described in the Supplementary Material [46].For direct comparison to the kinetic model, we also compute the fluid-dynamical evolution, where the energy conservation equation u ν ∂ µ T µν = 0 reduces to and it is closed by the following relaxation equations [39,48] with the first-order transport coefficients given by The second-order transport coefficients δ ΠΠ , λ Ππ , δ ππ , τ ππ , and λ πΠ are derived from the Shakhov collision term, as shown in the Supplementary Material [46].In order to validate the kinetic Shakhov model, we performed two sets of simulations with particles of mass m 0 = 1 GeV/c 2 .At initial time τ 0 = 0.5 fm/c and initial temperature β −1 0 = 0.5 GeV, we set f k (τ 0 ) = exp(−β 0 E k ).The results obtained using the AW model with τ π = τ Π = τ R = 0.5 fm/c, leading to 4πη/s ≃ 2.6 at the initial time, are used as reference.The solutions of the corresponding second-order fluid-dynamical equations (36)- (37) are shown in Fig. 1 with dashed black lines.Panels (a1)-(c1) of Fig. 1 show the results when τ π = 0.5 fm/c is fixed, while τ Π = τ R is varied.The red, green and blue lines with symbols correspond to the kinetic equation with τ Π /τ π = 3, 1 and 1/3, respectively.Panels (a1) and (c1) show that the evolution of P l /P ⊥ and π remains independent of τ Π .Panel (b1) shows that the amplitude of Π scales roughly with the ratio τ Π /τ π .Panels (a2)-(c2) depict the results for fixed τ Π ≡ τ R = 0.5 fm/c and variable τ π .During an initial period of time, (b2) indicates a non-linear dependence of Π on τ π .However, at late times, all curves converge towards an asymptotic solution governed by the value of τ Π .Panels (a2) and (c2) display the expected dependence of P l /P ⊥ and π on τ π , with a larger dip and a slower approach to equilibrium, P l /P ⊥ = 1 and π = 0, for larger values of τ π .The small discrepancies observable in Fig. 1 between the kinetic and fluid-dynamical results are expected due to the large values of the relaxation times and their associated viscosities [48,49,50].The results at smaller/larger viscosities are in better/worse agreement.

Example II: Longitudinal waves
We now consider sound waves propagating through an ultrarelativistic, classical ideal gas at rest, with temperature β −1 0 = 0.6 GeV and α 0 = 0. Since ζ = 0 when m 0 = 0, we are left with only two transport coefficients: κ and η.Setting τ π = τ R and τ V τ R in Eq. ( 22) and using J nq = Pβ 2−n (n + 1)!/[2(2q + 1)!!] leads to where f 0k = e α−βE k .At the initial time, the system is in equilibrium and the density n = n 0 + δn and velocity u z = δv are perturbed via δn = δn 0 cos(kz) and δv = δv 0 sin(kz), while the pressure is unperturbed, P = P 0 = n 0 T 0 .Here, k = 2π/L and L = 6.4 fm, while δn 0 /n 0 = δv 0 /c = 10 −3 .We then track the time evolution of the amplitudes δπ and δV, computed as where δV = V z and δπ ≃ π zz correspond to the diffusion current and shear-stress tensor, respectively.For more details, see Sec.V of Ref. [44].We now compare the numerical solution of the full Shakhov model described by Eqs. ( 12) and (39), obtained as described in the Supplementary Material [46], with the solution of the linearized second-order fluid dynamical equations, computed using the methods of Refs.[15,44].Panel (a) shows the ratio δV/(δn 0 kτ V ) for three values of τ V /τ π , with kτ π = 0.4 fixed, such that 4πη/s ≃ 3.11.It can be seen that the amplitude of δV and the damping exponent scale almost linearly with τ V .Similarly, for the δπ/(δv 0 kτ π ) ratio shown in panel (b), we varied τ π /τ V while keeping kτ V = 0.4 fixed.This panel shows that the damping exponent of δπ is proportional to τ π .Small discrepancies between the kinetic and fluid-dynamical results when τ V /τ π = 3 in panel (a) or τ π /τ V = 3 in panel (b) are due to the departure from the domain of applicability of second-order fluid dynamics [15].

Conclusion
In this paper, we introduced the relativistic Shakhov-type model as an extension of the Anderson-Witting RTA with momentum-independent relaxation time.This model introduces new and independently-adjustable first-order transport coefficients/relaxation times for the bulk viscous pressure, diffusion current, and shear-stress tensor.Although here our primary focus was on the first-order transport coefficients, the Shakhov model can be systematically extended [51] to control all transport coefficients of second-order fluid dynamics.δV /(δn Hydro This more general kinetic equation makes it possible to study different values or parameterizations of the first-order transport coefficients in various phenomena such as the baryon diffusion or the effect of bulk and shear viscosities similarly as in the state of the art second-order fluid dynamical simulations.Therefore, the present work significantly extends the applicability of kinetic RTA models to describe systems relevant in high-energy heavy-ion collisions. Shakhov model for the Bjorken flow.-In the case of the Bjorken expansion, we considered a massive, ideal, uncharged gas, such that α (0) r is given by Eq. ( 10).The first-order transport coefficients ζ and η are listed in Eqs.(38).The second-order transport coefficients appearing in Eq. ( 37) are listed here from Ref. [39]: Since the Shakhov distribution employed in Eq. ( 34) uses τ Π = τ R , the coefficients R (0) −r,0 reduce to their corresponding values for the AW model, namely where Eq. ( 10) was employed to replace α (0) r .On the other hand, R (2)  −r,0 becomes which, in the limit of τ Π = τ π , recovers the analogous coefficient appearing in the AW model, α (2)  −r /α (2) 0 = J 3−r,2 /J 32 .Therefore, the transport coefficients λ Ππ , δ ππ , and τ ππ involving R (2)   −2,0 are modified with respect to their AW expressions, while δ ΠΠ and λ πΠ remain unchanged.

SM-2. Entropy production
We now discuss the thermodynamic consistency of the Shakhov model by considering the entropy production where S µ = − dK k µ ( f k ln f k + a fk ln fk ) is the entropy fourcurrent.As originally pointed out by Shakhov [18], asserting the sign of ∂ µ S µ for arbitrary distributions f k is difficult, but if the fluid is not far from equilibrium, quadratic terms in δ f k or δ f Sk can be neglected and the logarithm in Eq. (SM-17) can be approximated as: where ϕ k = δ f k /( f 0k f0k ).Thus, Eq. (SM-17) becomes where on the second line, we have used the relation (δ , the first term on the right-hand side of the equation vanishes due to the matching conditions in Eq. ( 19).
The second term can be estimated using Eq. ( 13), leading to and with Eq. ( 17) confirms the second law of thermodynamics,

SM-3. Numerical method for the Shakhov model
To solve the Shakhov kinetic model we employ a discrete velocity method inspired by the Relativistic Lattice Boltzmann algorithm of Refs.[16,57,58,59,60].We consider the rapidity-based moments of f k introduced in Ref. [39], which eliminates two out of the three dimensions of the momentum space for the particular case of the (1 + 1)-dimensional longitudinal waves SM-3.1, and the (0 + 1)-dimensional boost invariant expansion SM-3.2, respectively.

SM-3.1. Longitudinal wave damping problem
In the application of Sec. 6, the fluid is homogeneous with respect to the x and y directions.Parameterizing the momentum space using (m ⊥ , φ ⊥ , v z ) as in Ref. [39], (SM-22) the Boltzmann equation with the Shakhov model for the collision term reduces to where Eq. (SM-23) becomes It can be shown [51] that the macroscopic quantities N t , N z , T tt , T tz and T zz can be obtained from F 1 and F 2 via For the case of massless particles considered in Sec.6, T µ µ = 0, such that T xx = T yy = (T tt − T zz )/2.From the above, it is clear that the time evolution of both N µ and T µν is fully determined by the functions F 1 and F 2 .In order to solve Eq. (SM-25), the functions F S n must be obtained by integrating Eq. ( 39), yielding: (SM-27) The time discretization is performed using equal time steps δt = 10 −3 fm/c and the time stepping is performed using the third-order total variation diminishing (TVD) Runge-Kutta scheme [61,62].The spatial domain [−L/2, L/2] is discretized using S = 100 cells of size δs = L/S , centred on z s = (s − 1 2 )δs − L 2 , 1 ≤ s ≤ S .The spatial derivative v z ∂ z F n is approximated using finite differences: where F n;s+1/2 represents the flux at the interface between cells s and s + 1.For definiteness, we compute this flux using the upwind-biased fifth-order weighted essentially non-oscillatory (WENO-5) scheme introduced in Ref. [63,64].Finally, the v z momentum space coordinate is discretized via the Gauss-Legendre quadrature with K = 20 points, such that P K (v z j ) = 0, with 1 ≤ j ≤ K and P K (z) being the Legendre polynomial of order K.Then, integrals with respect to v z of a function g(v z ) are approximated via 1 −1 dv z g(v z ) ≃ K j=1 g j , g j ≡ w j g(v z j ), (SM-29) with w j being the Gauss-Legendre quadrature weights [65].

SM-3.2. Algorithm for the Bjorken expansion
The algorithm for the Bjorken expansion is identical to that described in Ref. [39], hence we only recall the main method here.The spatial rapidity is η s = artanh(z/t), and the parametrization of the momentum space is as in Eq. (SM-22), with (k t , k z ) replaced by (k τ , τk η s ), where (SM-30) Retaining the definition (SM-24) of the rapiditybased moments, the non-vanishing components of T µν = diag(e, P ⊥ , P ⊥ , τ −2 P l ) are given by [39] e P l = The Boltzmann equation becomes where F S n can be obtained by integrating Eq. ( 34): z .Since π = 2 3 (P ⊥ − P l ), the system is closed in terms of F 0 and F 2 .The time integration and v z discretization proceed as in the previous subsection, however now the time step δτ n = τ n+1 − τ n is determined adaptively via δτ n = min(α τ τ n , α R τ R ), (SM-34) with α τ = 10 −3 and α R = 1/2.Furthermore, the derivative with respect to v z is performed by projecting F n onto the space of Legendre polynomials, as described in Ref. [58].Considering the discretization with K = 20 discrete velocities 1 ≤ j ≤ K, we have where the kernel K j, j ′ is given in Eq. (3.54) of Ref. [58].