Study of transport coefficients in ultrarelativistic kinetic theory

TL;DR

This paper investigates transport coefficients in relativistic kinetic theory using analytical solutions, numerical Boltzmann equation solutions, and polynomial approximations, confirming consistency with Chapman-Enskog predictions at small relaxation times.

Contribution

It provides a comprehensive analysis of transport coefficients in relativistic hydrodynamics, combining analytical, numerical, and approximation methods to validate theoretical predictions.

Findings

Transport coefficients agree with Chapman-Enskog predictions at small relaxation times.

Numerical solutions of the Boltzmann equation confirm analytical results.

Orthogonal polynomial approximations recover Chapman-Enskog results exactly.

Abstract

A spatially-periodic longitudinal wave is considered in relativistic dissipative hydrodynamics. At sufficiently small wave amplitudes, an analytic solution is obtained in the linearised limit of the macroscopic conservation equations within the first- and second-order relativistic hydrodynamics formulations. A kinetic solver is used to obtain the numerical solution of the relativistic Boltzmann equation for massless particles in the Anderson-Witting approximation for the collision term. It is found that, at small values of the Anderson-Witting relaxation time $\tau$, the transport coefficients emerging from the relativistic Boltzmann equation agree with those predicted through the Chapman-Enskog procedure, while the relaxation times of the heat flux and shear pressure are equal to $\tau$. These claims are further strengthened by considering a moment-type approximation based on orthogonal polynomials under which the Chapman-Enskog results for the transport coefficients are exactly recovered.

Tables

Table 1. Table 1: Reference quantities employed for the non-dimensionalisation procedure employed in this paper.

Reference quantity	Conventional value
${\hat{L}}_{\mathrm{ref}}$	$\hat{L}=2\pi /\hat{k}$
${\hat{v}}_{\mathrm{ref}}$	$\hat{c}$ (speed of light in vacuum)
${\hat{n}}_{\mathrm{ref}}$	${\hat{n}}_0$
${\hat{T}}_{\mathrm{ref}}$	${\hat{P}}_0/{\hat{K}}_B{\hat{n}}_0$
${\hat{t}}_{\mathrm{ref}}$	${\hat{L}}_{\mathrm{ref}}/{\hat{c}}_{\mathrm{ref}}$
${\hat{p}}_{\mathrm{ref}}$	$\widehat{K_B}{\hat{T}}_{\mathrm{ref}}/\hat{c}$
${\hat{f}}_{\mathrm{ref}}$	${\hat{n}}_{\mathrm{ref}}/{\hat{p}}_{\mathrm{ref}}^3$

Equations

$$p^{\mu }{\partial }_{\mu }f=\frac{p\cdot u_L}{\tau }(f-f_L^{(\mathrm{eq})}),$$

$$f_L^{(\mathrm{eq})}=\frac{n_L}{8\pi T_L^3}\exp \left(\frac{p\cdot u_L}{T_L}\right).$$

$$N^{\mu }=\int \frac{d^{3}p}{p^0}f{p}^{\mu },T^{\mu \nu }=\int \frac{d^{3}p}{p^0}f{p}^{\mu }{p}^{\nu }.$$

$$N_{(\mathrm{eq})}^{\mu }=n_L{u}_L^{\mu },T_{(\mathrm{eq})}^{\mu \nu }=(E_L+P_L)u_L^{\mu }{u}_L^{\nu }+P_L{\eta }^{\mu \nu }.$$

$$T^{\mu }{}_{\nu }{u}_L^{\nu }=-E_L{u}_L^{\mu }.$$

$$n_L=-N_{\mu }{u}_L^{\mu }.$$

$${\partial }_{\mu }{N}^{\mu }=0,{\partial }_{\nu }{T}^{\mu \nu }=0.$$

$$u^{\mu }=N^{\mu }/\sqrt{-N^2}.$$

$$N^{\mu }=n{u}^{\mu },$$

$$T^{\mu \nu }=E{u}^{\mu }{u}^{\nu }+(P+\overline{\omega }){\Delta }^{\mu \nu }+u^{\mu }{q}^{\nu }+q^{\nu }{u}^{\mu }+{\Pi }^{\mu \nu },$$

$$n=-u_{\mu }{N}^{\mu },E=u_{\mu }{u}_{\nu }{T}^{\mu \nu },P=\frac{1}{3}{\Delta }_{\mu \nu }{T}^{\mu \nu },$$

$${\Pi }^{\mu \nu }=\left({\Delta }^{\mu }{}_{\lambda }{\Delta }^{\nu }{}_{\kappa }-\frac{1}{3}{\Delta }^{\mu \nu }{\Delta }_{\lambda \kappa }\right)T^{\lambda \kappa },$$

$$q^{\mu }=-{\Delta }^{\mu }{}_{\nu }{u}_{\lambda }{T}^{\nu \lambda },$$

$${\partial }_{t}f+\xi {\partial }_{z}f=-\frac{{\gamma }_L(1-{\beta }_L\xi )}{\tau }(f-f_L^{(\mathrm{eq})}),$$

$$f_L^{(\mathrm{eq})}=\frac{n_L}{8\pi T_L^3}\exp \left[-\frac{p{\gamma }_L}{T_L}(1-{\beta }_L\xi )\right],$$

$$u^{\mu }{\partial }_{\mu }=\gamma ({\partial }_t+\beta {\partial }_z),q^{\mu }{\partial }_{\mu }=q(\beta {\partial }_t+{\partial }_z),$$

$${\Pi }^{\mu \nu }=\Pi \left(\begin{array}{cccc}{\beta }^2{\gamma }^2& 0& 0& \beta {\gamma }^2\\ 0& -\frac{1}{2}& 0& 0\\ 0& 0& -\frac{1}{2}& 0\\ \beta {\gamma }^2& 0& 0& {\gamma }^2\end{array}\right),$$

$$E_L=$$

$${\beta }_L=$$

$$n=n_0+\delta n,P=P_0+\delta P,$$

$$N^{\mu }\simeq (n_0+\delta n,0,0,n_0\beta {)}^T,T^{\mu \nu }\simeq \left(\begin{array}{cccc}3(P_0+\delta P)& 0& 0& 4\beta P_0+q\\ 0& P_0+\delta P-\frac{\Pi }{2}& 0& 0\\ 0& 0& P_0+\delta P-\frac{\Pi }{2}& 0\\ 4\beta P_0+q& 0& 0& P_0+\delta P+\Pi \end{array}\right),$$

$$n_L\simeq$$

$$P_L\simeq$$

$${\beta }_L\simeq$$

$${\partial }_t\delta n+n_0{\partial }_z\beta =0,$$

$$3{\partial }_t\delta P+4P_0{\partial }_z\beta +{\partial }_{z}q=0,$$

$$4P_0{\partial }_t\beta +{\partial }_{t}q+{\partial }_z\delta P+{\partial }_z\Pi =0.$$

$${\partial }_{t}f+\xi {\partial }_{z}f\simeq -\frac{1}{\tau }(f-f_L^{(\mathrm{eq})}),$$

$$\begin{gathered} f_L^{(\mathrm{eq})}\simeq \frac{n_0}{8\pi T_0^3}{e}^{-p/T_0}[1+\frac{p\xi }{T_0}\left(\beta +\frac{q}{4P_0}\right) \\ +\frac{4\delta n}{n_0}-\frac{3\delta P}{P_0}+\frac{p}{T_0}\left(\frac{\delta P}{P_0}-\frac{\delta n}{n_0}\right)]. \end{gathered}$$

$$\left(\begin{array}{c}\beta \\ q\end{array}\right)=\left(\begin{array}{c}\tilde{\beta }\\ \tilde{q}\end{array}\right)\sin kz,\left(\begin{array}{c}\delta n\\ \delta P\\ \Pi \end{array}\right)=\left(\begin{array}{c}\widetilde{\delta n}\\ \widetilde{\delta P_{\alpha }}\\ \widetilde{{\Pi }_{\alpha }}\end{array}\right)\cos kz,$$

$$\tilde{M}=\sum _{\alpha }{M}_{\alpha }{e}^{-\alpha t},$$

$$\alpha \delta n_{\alpha }-k{n}_0{\beta }_{\alpha }=0,$$

$$3\alpha \delta P_{\alpha }-4k{P}_0{\beta }_{\alpha }-k{q}_{\alpha }=0,$$

$$4\alpha P_0{\beta }_{\alpha }+\alpha q_{\alpha }+k\delta P_{\alpha }+k{\Pi }_{\alpha }=0.$$

$$\tilde{\beta }(t=0)=$$

$$\widetilde{\delta n}(t=0)=$$

$$\widetilde{\delta P}(t=0)=$$

Full text

I Introduction

The relativistic Boltzmann equation is known to reduce to the equations of relativistic hydrodynamics in the limit when the mean free path of the particle constituents is negligible compared to the typical length scales of the system cercignani02 . The transition from kinetic theory to relativistic hydrodynamics is traditionally performed following two approaches: the Chapman-Enskog procedure and Grad’s 14 moments approximation cercignani02 . These two approaches yield different expressions for the transport coefficients appearing in the constitutive equations of the underlying hydrodynamic equations. While the non-relativistic limit of these expressions coincides between the two formulations, their ultrarelativistic limits differ.

In order to check which of the two approaches (Chapman-Enskog expansion or the Grad method) correctly predicts the transport coefficients of the hydrodynamic equations, a solution of the relativistic Boltzmann equation is required. Solving the Boltzmann equation requires an explicit expression for the collision term, which in general is an integral operator taking into account local binary collisions. A considerable simplification arises by employing a single relaxation time (SRT) approximation. The most common SRT approximations are the Marle marle69 and the Anderson-Witting anderson74a ; anderson74b models, which generalise the widely-used Bhatnagar-Gross-Krook (BGK) model introduced in Ref. bhatnagar54 for the non-relativistic case. Since it is known that the Marle model is not appropriate for the study of the flow of massless particles cercignani02 ; anderson74a , only the Anderson-Witting model will be considered in this paper.

There has been recent evidence in the literature indicating that the transport coefficients predicted by the Chapman-Enskog method are closer to those recovered from solutions of the Boltzmann equation than those obtained through Grad’s 14 moment approximation.

Florkowski et al. obtained a solution of the Anderson-Witting-Boltzmann (AWB) equation in the case of Bjorken flow bjorken83 at non-vanishing relaxation time, written in integral form for the case when the distribution function depends only on proper time. This solution was restricted to the massless case in Refs. florkowski13a ; florkowski13b and extended to the massive case in Ref. florkowski14 . In Refs. florkowski13a ; florkowski13b , it was shown that the numerical solution of the Israel-Stewart equations israel79 leads to better agreement with the AWB solution (computed also numerically) when the Chapman-Enskog value for the shear viscosity $\eta$ is used compared to when the 14 moment approximation is used. The same conclusion is reached in Ref. florkowski15 for the case of Bose-Einstein and Fermi-Dirac statistics.

The solution of the AWB equation describing the Bjorken flow of massive particles obtained in Ref. florkowski14 is used in Ref. ryblewski15 to highlight that the second-order Chapman-Enskog expansion is in closest agreement to the solution of the AWB equation as compared to the Israel-Stewart and the 14 moment approaches introduced in Refs. israel79 and denicol14 , respectively.

In Ref. bhalerao14 , Bhalerao et al demonstrated that a Chapman-Enskog-like approximation of the non-equilibrium distribution function yielded a closer agreement with the inviscid limit of the Bjorken flow than the Grad approximation.

The relativistic lattice Boltzmann (LB) method has also been used as a tool to solve the AWB equation mendoza10prl ; romatschke11 ; hupp11 ; mohseni13 ; mendoza13prd ; blaga17prc ; blaga17aip ; gabbana17a ; gabbana17b ; coelho17 . The propagation of planar shock waves of massless mendoza10prl ; hupp11 ; mohseni13 ; mendoza13prd ; blaga17prc and massive gabbana17a ; gabbana17b particles was investigated using the LB method and the results were validated by comparison with the data obtained using the Boltzmann Approach to Multi Parton Scattering (BAMPS) reported in Refs. bouras09prl ; bouras09nucl ; bouras10 . In order to compare the LB results obtained in the frame of the AWB equation and the BAMPS results, the shear viscosity to entropy ratio $\eta/s$ must be kept constant. Matching within the LB method the values of $\eta/s$ employed in the BAMPS simulations requires as an input the exact expression for the shear viscosity $\eta$, which can be obtained via Grad’s 14 moment approach mendoza10prl ; hupp11 ; mohseni13 ; mendoza13prd ; gabbana17a or the Chapman-Enskog procedure blaga17prc ; gabbana17b . The authors of Refs. blaga17prc ; gabbana17b note that employing the Chapman-Enskog value for $\eta$ leads to better agreement with the BAMPS data than when the Grad value is employed.

Recently, the relativistic lattice Boltzmann model developed in Ref. gabbana17a was used in Ref. gabbana17b to study the dissipative attenuation of the relativistic equivalent of the Taylor-Green vortices. This study allowed the authors of Ref. gabbana17b to demonstrate that the correct value of the shear viscosity is that given by the Chapman-Enskog procedure rather than the Grad method for a wide range of particle masses.

In this paper, a study of the transport coefficients and second-order relaxation times arising in the Anderson-Witting model is presented, by considering the dissipative attenuation of a harmonic longitudinal wave. Three regimes of wave propagation will be analysed in this paper, corresponding to the cases when the velocity $\beta$ along the wave propagation direction (Case 1), pressure $P=P_{0}+\delta P$ (Case 2a) or density $n=n_{0}+\delta n$ (Case 2b) are perturbed harmonically about $\beta=0$, $n=n_{0}$ and $P=P_{0}$. In each case, the other two macrosopic variables are left unperturbed in the initial state. In all cases, at initial time, the fluid is assumed to be in local thermodynamic equilibrium characterised by the Maxwell-Jüttner distribution corresponding to the local values of $n$, $P$ and $\beta$. The analytic analysis of this system is restricted to the regime of small amplitudes $\beta_{0}$, $\delta n_{0}$ and $\delta P_{0}$, where the linearised form of the macroscopic equations can be solved exactly.

In the first-order description (i.e. the five field approximation), the wave amplitude predicted by the analytic solution consists of a damped, oscillatory term (with respect to time) which allows the wave to propagate at approximately the speed of sound, for which the attenuation coefficient $\alpha_{d}$ is directly proportional to the shear viscosity $\eta$. The second term is non-oscillatory (evanescent) and its attenuation coefficient $\alpha_{\lambda}$ is directly proportional to the heat conductivity $\lambda$. Depending on the initial conditions, this system allows $\eta$ and $\lambda$ to be measured separately. It can be shown that the wave corresponding to Case 1 propagates adiabatically (with no heat flux being present), allowing $\eta$ to be measured independently. Furthermore, the heat flux $q$ is purely evanescent (has no oscillatory contribution), such that its time evolution is completely determined by $\alpha_{\lambda}$, being independent of $\alpha_{d}$ and hence of $\eta$. Cases 2a and 2b will be therefore used to measure $\lambda$ independently of $\eta$.

A known fundamental limitation of the first-order theory is that it allows the non-causal instantaneous response in the heat flux $q$ and shear pressure $\Pi$ induced by changes in the gradients of the fundamental variables $n$, $\beta$ and $P$. In particular, this theory does not allow the values of $q$ and $\Pi$ to be set at initial time $t=0$ independently of the values of $n$, $\beta$ and $P$. This is incompatible with the initial local thermodynamic equilibrium state considered in this paper, in which $q=\Pi=0$, such that the first-order theory prediction for the evolution of $q$ and $\Pi$ is not accurate for a duration of time proportional to the Anderson-Witting relaxation time $\tau$. In the second-order hydrodynamics approach, both $q$ and $\Pi$ obey independent evolution equations which, for sufficiently small values of $\tau$, allow them to relax from arbitrary initial configurations to the first-order predictions on time scales given by the relaxation times $\tau_{q}$ and $\tau_{\Pi}$, respectively. These relaxation times will be studied for Case 1 and Case 2b, where it will be demonstrated that $\tau_{q}\simeq\tau$ and $\tau_{\Pi}\simeq\tau$ at small values of $\tau$.

In the analysis of Case 2b, a more subtle limitation of the first-order theory was encountered. If the initial state is prepared as described in Case 2b, the numerical simulations indicate that the pressure perturbation $\delta P$ and the shear pressure $\Pi$ remain zero throughout the evolution of the wave, for all tested values of the relaxation time, provided the initial perturbation $\delta n_{0}$ is small. Furthermore, the decay of the amplitudes of $\delta n$, $\beta$ and $q$ is strictly exponential, with no oscillations. While the above behaviour is successfully recovered in the second-order theory, in the first-order theory, $\delta P$ and $\Pi$ are non-zero. Moreover, the amplitude of $\beta$ is oscillatory and the oscillation amplitude is of the same order of magnitude as the non-oscillatory $\beta$ predicted by the second-order theory.

Next, a moment-based approach similar to the one introduced in Refs. romatschke11 ; blaga17prc ; blaga17aip ; denicol12 is considered, where the distribution function is expanded with respect to the Laguerre and Legendre polynomials, corresponding to the magnitude $p$ of the particle momentum and the $z$ component $\xi=p^{z}/p$ of its velocity, respectively. Retaining zeroth and first-order terms with respect to the Laguerre polynomials and terms up to second-order with respect to the Legendre polynomials, a system of $6$ evolution equations is obtained for the density $n$, pressure $P$, velocity $\beta$, heat flux $q$, shear pressure $\Pi$ and an extra non-hydrodynamic variable. In the frame of this moment-based model, the Chapman-Enskog predictions for the shear viscosity $\eta$ and heat conductivity $\lambda$ are exactly recovered. As highlighted in Ref. denicol12 , there is a fundamental difference between the above proposed moment-based method and Grad’s 14 moment approach, due to the fact that the former is based on an expansion with respect to orthogonal polynomials, while the latter relies on an expansion on polynomials in $p^{\mu}$, which do not constitute an orthogonal basis. An analytic solution obtained within the above model is employed to show that at large values of the relaxation time $\tau$, the moment-based approach provides a better analytic description of the evolution of the heat flux compared to the second-order hydrodynamics result.

Finally, the ballistic regime is analysed, where the particle constituents stream freely. Since the flow is now collisionless, no dissipation occurs and the wave attenuation is no longer exponential. Instead, the dispersive regime sets in, since now the wave can be regarded as a packet which consists of non-interacting constituents which propagate in the longitudinal direction at different velocities $\xi$. An analytic solution for the linearised limit of the ballistic regime is presented, with the aid of which the capability of the numerical code employed in this paper to capture the free-streaming dynamics is demonstrated.

According to Refs. blaga17prc ; blaga17aip , high order quadratures (i.e. large velocity sets) are required to obtain accurate simulation results at large values of $\tau$, when the flow is out of equilibrium and rarefaction effects become important. This is performed in a straightforward manner following the procedure described in Ref. blaga17prc and summarised in Appendix B. The numerical experiments presented in this paper were therefore conducted using the quadrature-based R-SLB models developed in Ref. blaga17prc . While the analysis presented herein is restricted to the case of massless particles, it can be easily extended to the case of massive particles, e.g. following Refs. gabbana17a ; gabbana17b ; romatschke12 .

The paper is organised as follows. The general framework for the study of the propagation of longitudinal waves is introduced in Sec. II by linearising the relativistic hydrodynamics equations with respect to the wave amplitude. The relativistic Boltzmann equation in the Anderson-Witting approximation for the collision term (the AWB equation) is also briefly presented, alongside a description of the Landau frame. In Secs. III and IV, the longitudinal wave problem is considered from the perspective of the first- and second-order hydrodynamics theories, respectively, while in Sec. V, the same problem is considered using a moment-based approach. In all cases, numerical simulations are employed to study the validity and applicability of these theories as the relaxation time is increased. In Sec. VI, the propagation of the longitudinal wave is analysed analytically and numerically in the ballistic regime. A short description of the numerical method employed in this paper is provided in Appendix B.

Throughout this paper, the metric convention $\eta_{\mu\nu}={\rm diag}(-1,1,1,1)$ is employed. The non-dimensionalisation convention is presented in Appendix A and summarised in Tab. 1.

II Relativistic fluid dynamics

In this section, the common framework used in later sections for the analysis of the evolution of longitudinal waves is presented. Subsection II.1 introduces a brief review of the connection between the relativistic Boltzmann equation and the macroscopic hydrodynamics equations, which are written in linearised form in Subsec. II.2. The equations which serve as the basis for the analysis of longitudinal waves are presented in Subsec. II.3.

II.1 Relativistic kinetic theory

This paper is focused on the relativistic Boltzmann equation for massless particles in the Anderson-Witting approximation for the collision term, which reads anderson74a :

[TABLE]

where it is assumed for simplicity that the relaxation time $\tau$ is constant. The equilibrium distribution $f^{(\mathrm{eq})}$ is taken to be the Maxwell-Jüttner distribution function:

[TABLE]

In the above, $n_{L}$ represents the particle number density, $u_{L}^{\mu}$ is the macroscopic four-velocity, $T_{L}$ is the local temperature and $p^{\mu}$ is the on-shell particle four-momentum. The quantities bearing the subscript $L$ are expressed in the Landau (energy) frame anderson74a ; landau87 .

The transition from the Boltzmann equation (1) to relativistic hydrodynamics is done by considering the macrosopic four-flow vector $N^{\mu}$ and stress-energy tensor (SET) $T^{\mu\nu}$, which are obtained by integrating the distribution function over the momentum space:

[TABLE]

Substituting $f^{(\mathrm{eq})}_{L}$ (2) into Eq. (3) gives the equilibrium four-flow vector $N^{\mu}_{(\mathrm{eq})}$ and SET $T^{\mu\nu}_{(\mathrm{eq})}$:

[TABLE]

The Landau velocity $u_{L}^{\mu}$ is defined as the eigenvector of $T^{\mu\nu}$ corresponding to the Landau energy density $E_{L}$:

[TABLE]

For massless particles, $E_{L}=3P_{L}$ and the Landau pressure $P_{L}=n_{L}T_{L}$ is used to define the Landau temperature $T_{L}$, while the Landau particle number density $n_{L}$ is obtained by contracting $N^{\mu}$ with $u^{\mu}_{L}$:

[TABLE]

Multiplying the Boltzmann equation (1) by the collision invariants $\psi\in\{1,p^{\mu}\}$ and integrating with respect to the momentum space, the following conservation equations are obtained:

[TABLE]

Due to its simplicity and pedagogical value, the Eckart (particle) frame will be employed in this paper, where the macroscopic velocity $u^{\mu}$ is defined as the unit vector parallel to $N^{\mu}$ eckart40 ; rezzolla13 :

[TABLE]

With respect to $u^{\mu}$, $N^{\mu}$ and the SET $T^{\mu\nu}$ can be decomposed as:

[TABLE]

where $\Delta^{\mu\nu}=\eta^{\mu\nu}+u^{\mu}u^{\nu}$ is the projector on the hypersurface orthogonal to $u^{\mu}$. The particle number density $n$, energy density $E$, isotropic pressure $P+\overline{\omega}$, heat flux $q^{\mu}$ and shear stress tensor $\Pi^{\mu\nu}$ can be obtained as follows bouras10 ; rezzolla13 :

[TABLE]

while the dynamic pressure $\overline{\omega}=0$ for massless particles, when $E=3P$. In general, the Eckart quantities introduced above are different from the corresponding quantities defined in the Landau frame (more details will be given in Sec. II.2).

The system (7) consisting of $5$ equations is not closed, since $q^{\mu}$ and $\Pi^{\mu\nu}$ are a priori unconstrained. The constitutive relations which close this system corresponding to the first- and second-order relativistic hydrodynamics frameworks will be discussed in Secs. III and IV.

II.2 Linearised relativistic hydrodynamics

Let us now consider a system which is homogeneous along the $x$ and $y$ directions. In this case, the AWB equation (1) reduces to blaga17prc :

[TABLE]

where $\xi=p^{z}/p$ represents the particle velocity along the $z$ axis, taking values in $[-1,1]$. Taking into account the constraints $u_{\mu}q^{\mu}=0$ and $u_{\mu}\Pi^{\mu\nu}=0$, the variables $u^{\mu}$, $q^{\mu}$ and $\Pi^{\mu\nu}$ can be taken as follows blaga17prc ; bouras10 :

[TABLE]

where $\gamma=(1-\beta^{2})^{-1/2}$ is the Lorentz factor corresponding to the velocity $\beta$. The Landau frame can be constructed analytically by solving the eigenvalue equation (5) blaga17prc ; bouras10 :

[TABLE]

In order to arrive at the linearised form of Eqs. (7), $n$ and $P$ can be written as:

[TABLE]

where $\delta n/n_{0}$ and $\delta P/P_{0}$ are quantities of order $O(\beta)$ and the limit $\beta\ll 1$ was considered. Furthermore, $\Pi$ and $q$ are also of order $O(\beta)$, since they represent non-equilibrium quantities. Neglecting the terms of order $\beta^{2}$, $N^{\mu}$ and $T^{\mu\nu}$ (9) reduce to:

[TABLE]

while the Landau quantities $n_{L}$, $P_{L}$ and $\beta_{L}$ can be approximated through:

[TABLE]

In the linearised approximation, the conservation equations (7) reduce to:

[TABLE]

Noting that $f-f^{(\mathrm{eq})}_{L}$ is also of order $O(\beta)$, the Boltzmann equation (11) can also be expressed in linearised form:

[TABLE]

where $f^{(\mathrm{eq})}_{L}$ can be linearised as follows:

[TABLE]

II.3 Longitudinal waves

Next, solutions of the following form are sought:

[TABLE]

where $k=2\pi/L$ is the wave number and $L$ is the wavelength. The quantities with a tilde $\widetilde{M}\in\{\widetilde{\beta},\widetilde{\delta n},\widetilde{\delta P},\widetilde{q},\widetilde{\Pi}\}$ depend only on time $t$. Taking this dependence in the form:

[TABLE]

Eq. (17) can be solved for each (constant) value of $\alpha$ independently, yielding a spectrum of linearly-independent modes satisfying:

[TABLE]

The imaginary part of $\alpha$ represents the propagation angular frequency, while its real part causes the dissipative dampening of the wave. In order to solve the above set of equations, the constitutive equations for $q$ and $\Pi$ must be supplied separately.

The initial conditions for Eqs. (22) are given in the form

[TABLE]

In this paper, the following sets of values for $\beta_{0}$, $\delta n_{0}$ and $\delta P_{0}$ will be considered:

• •

Case 1: $\delta n_{0}=\delta P_{0}=0$, $\beta_{0}\neq 0$;

• •

Case 2a: $\delta n_{0}=\beta_{0}=0$, $\delta P_{0}\neq 0$;

• •

Case 2b: $\beta_{0}=\delta P_{0}=0$, $\delta n_{0}\neq 0$.

In Secs. III and IV, the consitutive relations corresponding to the first- and second-order relativistic hydrodynamics will be employed. In Sec. V, a solution of Eqs. (22) will be constructed starting from the Boltzmann equation (18) written in linearised form.

III First-order hydrodynamics

The equations of first-order relativistic hydrodynamics represent the analogue of the Navier-Stokes-Fourier equations of non-relativistic hydrodynamics. In this formulation, the fields $n$, $u^{\mu}$ and $P$ are considered as fundamental variables. Since $u^{\mu}$ is normalised according to $u^{2}=-1$, the theory contains five independent fields and is sometimes referred to as the five field theory cercignani02 . In this first-order framework, the constitutive equations for the heat flux $q^{\mu}$ and shear stress tensor $\Pi^{\mu\nu}$ represent algebraic relations linking them to the gradients of the fundamental fields via the transport coefficients $\lambda$ (heat conductivity) and $\eta$ (shear viscosity), respectively. Since $q^{\mu}$ and $\Pi^{\mu\nu}$ respond instantaneously to changes in the fundamental fields, the ensuing system of equations is not hyperbolic rezzolla13 , rendering the theory non-causal. This issue can be remedied within the second-order relativistic hydrodynamics framework, as will be discussed in Sec. IV. This section is focused on determining $\lambda$ and $\eta$ by comparing the analytical and numerical results for the attenuation process occuring in the longitudinal wave problem described in Sec. II.3.

III.1 Constitutive relations

The constitutive equations for $q^{\mu}$ and $\Pi^{\mu\nu}$ can be written in the frame of the first-order relativistic hydrodynamics as cercignani02 ; rezzolla13 :

[TABLE]

where $\lambda$ and $\eta$ represent the coefficients of heat conductivity and shear viscosity $\eta$, respectively. At the level of the first-order hydrodynamics theory, it is not specified whether the macroscopic velocity $u^{\mu}$ appearing in the right hand side of the second line of Eq. (24) is defined in the Eckart or in the Landau frame. In this section, the Landau frame velocity will be considered, since this choice seems natural when the Anderson-Witting approximation is used for the collision term florkowski15 ; ryblewski15 ; denicol14 ; bouras10 ; denicol12 ; jaiswal13a ; jaiswal13b ; chattopadhyay15 ; jaiswal15 . In Sec. V, it will be shown that this choice arises naturally when a moment-based approach is used to solve the AWB equation (18).

The connection between the Boltzmann equation (1) and the constitutive equations given in Eq. (24) is commonly achieved via two paths: (a) the Chapman-Enskog expansion; and (b) Grad’s 14 moments approximation. In the ultrarelativistic regime considered in this paper, the transport coefficients $\eta$ and $\lambda$ are given by:

[TABLE]

where the dimensionless constants $\eta_{0}$ and $\lambda_{0}$ are obtained using Grad’s approximation and the Chapman-Enskog procedure as follows cercignani02 :

[TABLE]

The validity of the constitutive equations (24) and of the above expressions for the transport coefficients is limited to the hydrodynamic regime, i.e. when the relaxation time $\tau$ is sufficiently small.

In the linearised approximation introduced in Sec. II.2, the constitutive equations (24) reduce to:

[TABLE]

where $u_{L}^{\mu}\simeq(1,0,0,\beta+q/4P_{0})$ was used in the expression for $\Pi$.

III.2 Longitudinal waves solution

The modes $q_{\alpha}$ and $\Pi_{\alpha}$ appearing in Eq. (22) can be found from Eq. (27):

[TABLE]

Noting from Eq. (22a) that

[TABLE]

Eq. (22b) reduces to:

[TABLE]

According to Eq. (28), the first parenthesis vanishes only when $q_{\alpha}=0$. Thus, the solution

[TABLE]

corresponds to the only mode which dissipates heat. In this case, Eq. (22) can be used to obtain:

[TABLE]

It is remarkable that this mode induces no viscous dissipation.

Considering now that $q=0$, Eq. (22c) reduces to:

[TABLE]

The solution $\beta_{\alpha}=0$ is trivial since in this case $\delta n_{\alpha}=\delta P_{\alpha}=0$. Setting the quantity inside the parenthesis equal to zero yields the following allowed values for $\alpha$:

[TABLE]

where the dampening ($\alpha_{d}$) and oscillatory ($\alpha_{o}$) parts of $\alpha_{\pm}$ read:

[TABLE]

It is worth noting that the phase velocity $\alpha_{o}/k=c_{s}\sqrt{1-\frac{3\alpha_{d}^{2}}{k^{2}}}$ predicted in the first-order theory is smaller than the sound speed $c_{s}=1/\sqrt{3}$. The amplitudes of the density and pressure perturbations $\delta n_{\pm}$ and $\delta P_{\pm}$ are given in terms of the velocity amplitudes $\beta_{\pm}$ as follows:

[TABLE]

Taking into account the above allowed values for $\alpha$, the general solution (21) reads:

[TABLE]

In the above, $\beta_{\lambda}$, $\beta_{c}=\beta_{+}+\beta_{-}$ and $\beta_{s}=-i(\beta_{+}-\beta_{-})$ are independent integration constants with respect to which the following definitions were made:

[TABLE]

while $q_{c}=q_{s}=0$. The other constants $\delta n_{\lambda}$, $\delta P_{\lambda}$, $\Pi_{\lambda}$ and $q_{\lambda}$ were already defined in Eq. (32).

The constants $\beta_{\lambda}$, $\beta_{c}$ and $\beta_{s}$ can be obtained by substituting the solution (37) into the initial conditions (23) yielding:

[TABLE]

The solution of Eq. (39) can be written as:

[TABLE]

The analytic solution presented in this section facilitates the study of the transport coefficients corresponding to a relativistic gas. Using the numerical method described in Appendix B, this system will be considered in the following subsections for the study of the ultrarelativistic limits of the shear viscosity $\eta$ and heat conductivity $\lambda$ arising from the AWB equation (1).

III.3 Case 1: Adiabatic flow

First, an adiabatic flow is considered (i.e. $q=0$) such that the shear viscosity $\eta$ can be isolated from the heat conductivity $\lambda$. This can be achieved when $3n_{0}\delta P_{0}=4P_{0}\delta n_{0}$. This condition is equivalent to the requirement that the fugacity $\lambda_{\rm fug}=n^{4}/P^{3}$ is constant in the initial state. Indeed, combining the first two relations in Eq. (17) gives:

[TABLE]

The above equation (valid in the linearised regime) shows that if there is no heat flux present, the fugacity remains constant in time.

Furthermore, Eq. (40) indicates that $\beta_{\lambda}=0$ when $3n_{0}\delta P_{0}=4P_{0}\delta n_{0}$, while Eq. (32) implies that $\delta n_{\lambda}$, $\delta P_{\lambda}$, $q_{\lambda}$ and $\Pi_{\lambda}$ cancel. Thus, the evolution of the fluid is completely independent of $\alpha_{\lambda}$, enabling $\eta$ to be determined independently. For simplicity, the initialisation corresponding to Case 1 in Sec. II.3 will be considered (i.e., $\delta n_{0}=\delta P_{0}=0$). In this limit, Eq. (40) reduces to:

[TABLE]

such that the exact solution (37) reads:

[TABLE]

where $\alpha_{d}$ and $\alpha_{o}$ are given in terms of $\eta$ in Eq. (35).

The analytic results in Eq. (43) are represented in Fig. 1 for the initial conditions $\beta_{0}=10^{-3}$ and $\delta n_{0}=\delta P_{0}=0$ alongside the corresponding numerical results obtained using the method described in Appendix B. The relaxation time was taken to be $\tau=0.0083$, such that both the dampening and the oscillatory characteristics of the solutions can be highlighted on the same timescale. The first entry in the legend (fine dotted lines) corresponds to the analytic expressions in Eqs. (43), where $\alpha_{d}$ and $\alpha_{o}$ are computed using the Chapman-Enkog value for $\eta$. The numerical results are indistinguishable from the analytic predictions.

Also in the plots in Fig. 1, the dampening caused by the $\exp(-\alpha_{d}t)$ factor in Eqs. (43) is represented when $\alpha_{d}$ is calculated using the Chapman-Enskog and Grad expressions for $\eta$. In the amplitude of the dampening terms in $\widetilde{\beta}/\beta_{0}$ and $\widetilde{\Pi}/\beta_{0}$, the approximation $\sqrt{1+\alpha_{d}^{2}/\alpha_{o}^{2}}\simeq 1$ was used. It can be seen that the dampening predicted by the analytic solution when the Grad expression for $\eta$ is used does not match the numerical results.

III.4 Cases 2a and 2b: Non-adiabatic flow

The coefficient $\alpha_{\lambda}$ can be investigated most easily by considering the decay of the amplitude $\widetilde{q}$ of the heat flux. The system will thus be initialised accordint to Case 2a ($\delta n_{0}=\beta_{0}=0$ and $\delta P_{0}=10^{-3}$) and Case 2b ($\delta P_{0}=\beta_{0}=0$ and $\delta n_{0}=10^{-3}$) described in Sec. II.3, while $\tau=0.0083$. According to Eqs. (32) and (40), $\widetilde{q}$ takes the following form:

[TABLE]

The above analytic result is compared in Fig. 2 to the numerical results obtained using the method described in Appendix B. For each of the two cases mentioned above, three curves are represented. The numerical results (dashed lines and points) are overlapped with the analytic prediction (44) when $\alpha_{\lambda}$ is calculated using the Chapman-Enskog expression for $\lambda$ (continuous line). The analytic prediction (44) corresponding to the case when $\alpha_{\lambda}$ is obtained using the Grad expression for $\lambda$ (dashed line) is clearly not consistent with the numerical results. The points at $t=0$ indicate that in the numerical simulations, the system was initialised using an equilibrium state, in which the heat flux vanishes. In the five-field theory, the initial value of $\widetilde{q}$ does not represent a free parameter, as can be seen from the solution in Eq. (44). However, since $\tau$ is small, the system quickly relaxes towards the five-field theory prediction (44). This relaxation process will be further considered in Sec. IV.

III.5 Limits of the linearised hydrodynamics equations

Next, an assessment of the limits within which the solution of the linearised equations (17) is applicable is performed. In order to reduce the rarefaction effects, $\tau$ is fixed at $\tau=0.0083$ throughout this subsection.

The solution (43) predicts that, for Case 1, the time evolution of $\widetilde{\beta}$ is damped according to the factor $\exp(-\alpha_{d}t)$, with $\alpha_{d}$ (35) being independent of the magnitude $\beta_{0}$ of the perturbation. Figure 3 shows that this is not the case: while at small values of $\beta_{0}$, $\widetilde{\beta}$ follows closely the analytic prediction [as confirmed in Fig. 1(c)], at larger values of $\beta$, the dampening is enhanced compared to the linearised limit (35).

In order to test the versatility of the functional form of the solution corresponding to the linearised regime, the parameters $\alpha_{d}$, $\alpha_{o}$ and $\alpha_{\lambda}$ are determined using nonlinear fits of the analytic solutions (43) and (44) to the corresponding numerical data. The coefficients $\alpha_{d}$ and $\alpha_{o}$ are obtained by treating them as free parameters while performing a two-parameter nonlinear fit of $\widetilde{\delta n}$, $\widetilde{\beta}$, $\widetilde{\delta P}$ and $\widetilde{\Pi}$ given in Eq. (43) for the initial conditions corresponding to Case 1 (i.e. $\delta n_{0}=\delta P_{0}=0$ and various values of $\beta_{0}$) to the corresponding numerical results. The coefficient $\alpha_{\lambda}$ is obtained by performing a one-parameter nonlinear fit of $\widetilde{q}$ with respect to the numerical results with the initial conditions described in Case 2a ($\beta_{0}=0$, $\delta n_{0}=0$ and various values for $\delta P_{0}$) and Case 2b ($\beta_{0}=0$, $\delta P_{0}=0$ and various values for $\delta n_{0}$).

The dependence of $\alpha_{d}$, $\alpha_{o}$ and $\alpha_{\lambda}$ on the amplitude of the perturbations is presented in the plots (a), (b) and (c) of Fig. 4, respectively. The horizontal axis in Fig. 4(c) represents the amplitude of the initial perturbation, i.e. $\delta P_{0}$ for Case 2a and $\delta n_{0}$ for Case 2b. All of the above plots show the analytic predictions (35) and (31) for $\alpha_{d}$, $\alpha_{o}$ and $\alpha_{\lambda}$, specialised to the cases when the transport coefficients $\eta$ and $\lambda$ are computed using the Chapman-Enskog (26b) and the Grad (26a) expressions. The results clearly favor the Chapman-Enskog expressions. These plots also indicate that the analytic analysis performed in Sec. III.2 in the context of the linearised hydrodynamic equations loses applicability when the perturbation amplitudes $\beta_{0}$, $\delta n_{0}/n_{0}$ or $\delta P_{0}/P_{0}$ are larger than $\sim 0.05$.

III.6 Limits of the first-order hydrodynamics regime

It is known that the constitutive equations (7) are valid only when $\tau$ is small cercignani02 . This subsection is focused on investigating the validity of the analysis presented in Sec. III.2 as the relaxation time $\tau$ is increased.

In order to test the effect of increasing $\tau$, the perturbations are kept at a small value, i.e. $\beta_{0}=10^{-3}$ (for Case 1), $\delta P_{0}=10^{-3}$ (for Case 2a), or $\delta n_{0}=10^{-3}$ (for Case 2b). The plots in Fig. 5 show the dependence of (a) $\alpha_{d}$, (b) $\alpha_{o}$, (c) $\alpha_{\lambda}$ and (d) the ratio $\lambda/\eta=2\alpha_{\lambda}n_{0}/3\alpha_{d}P_{0}$ on the value of $\tau$. As before, the analytic predictions for the dependence of these coefficients on $\tau$ is also shown for the cases when the transport coefficients $\eta$ and $\lambda$ are obtained using the Chapman-Enskog (26b) and Grad (26a) expressions. For $\tau<0.1$, the numerical results clearly favor the Chapman-Enskog expression. Plot (d) confirms that for small values of $\tau$, the ratio $\lambda/\eta$ is equal to $5/3$, as predicted in the Chapman-Enskog theory (26b). This value is in agreement with the high chemical potential limit of Fig. 2 in Ref. jaiswal15 .

While for $\tau\gtrsim 0.1$, the constitutive equations (7) no longer hold, our nonlinear fit analysis seems to indicate that the dampening coefficients $\alpha_{d}$ and $\alpha_{\lambda}$ plateau at large $\tau$. This conclusion is not necessarily meaningful, since the ansatz (20) that the time dependence of $\delta n$, $\delta P$ and $\beta$ is of the form $e^{-\alpha t}$ with constant $\alpha$ is not guaranteed to be valid in the transition regime. It is certain that the time dependence of the above quantities is more complex at large values of $\tau$, since the dissipative exponential attenuation is replaced by a polynomial dispersive attenuation in the ballistic regime, as will be shown in Sec. VI.

III.7 Summary

In this section, the attenuation of a longitudinal wave was analysed using the equations of first-order hydrodynamics. The results obtained by numerically solving the AWB equation at small values of the relaxation time $\tau$ showed an excellent agreement with the analytic results when the transport coefficinets $\lambda$ and $\eta$ were obtained using the Chapman-Enskog method. The analytic results corresponding to the transport coefficients obtained using the Grad method exhibited a clear discrepancy compared to the numerical results.

The validity of the linearised form of the AWB equation (18) was further tested by comparing the numerical and analytic results at increasing values of the wave amplitudes, while keeping $\tau=0.0083$. A visible discrepancy can be seen at wave amplitudes of $\beta_{0},\delta n_{0}/n_{0},\delta P_{0}/P_{0}\gtrsim 0.1$. The applicability of the functional form of the solution of the first-order hydrodynamics equations was further considered by numerically fitting the attenuation coefficients $\alpha_{\lambda}$ and $\alpha_{\eta}$, as well as the oscillation frequency $\alpha_{o}$, to the numerical results. The above analysis shows that the best fit parameters differ significantly from the analytic prediction when the wave amplitude is $\gtrsim 0.05$. At smaller values of the wave amplitudes, the numerical results clearly favored the Chapman-Enskog predictions for the transport coefficients, while the Grad predictions showed a visible discrepancy to the numerical results.

Finally, the validity of the first-order hydrodynamics equations was investigated at increasing values of $\tau$. The dampening coefficients $\alpha_{d}$ and $\alpha_{\lambda}$ were found to be directly proportional to $\tau$ only for $\tau\lesssim 0.1$, confirming the Chapman-Enskog prediction. At larger values of $\tau$, they increase at a significantly slower rate, signaling the breakdown of the first-order hydrodynamics formulation when $\tau\gtrsim 0.1$. The Grad prediction was in clear disagreement with the numerical results for all values of $\tau$.

IV Second-order hydrodynamics

The five-field theory provides the constitutive equations for $q^{\mu}$ and $\Pi^{\mu\nu}$ in the form given in Eq. (24). Since these constitutive equations do not represent evolution equations, $q^{\mu}$ and $\Pi^{\mu\nu}$ are fully determined by the spatial and temporal gradients of $n$, $u^{\mu}$ and $P$. In particular, their initial values at $t=0$ cannot be set to arbitrary values. In this section, the second-order extension of the five-field theory will be employed in order to study the relaxation process of the heat flux and shear stress from their initial vanishing value to the value required through the constitutive equations of the five-field theory.

IV.1 Constitutive relations

There are many variations of the form in which the equations of second-order hydrodynamics (also known as extended irreversible thermodynamics rezzolla13 ) are presented, essentially due to the route adopted in deriving them israel79 ; hiscock83 ; el10 ; denicol10 ; jaiswal13b ; chattopadhyay15 . Only the form introduced in Refs. hiscock83 ; rezzolla13 will be further considered, according to which $q^{\mu}$ and $T^{\mu\nu}$ satisfy the following equations:

[TABLE]

where the thermodynamic coefficients $\beta_{1}$, $\beta_{2}$ and $\gamma_{1}$ and the viscous-heat flux coupling coefficient $\alpha_{1}$ are a priori not known. After performing the linearisation described in Sec. II.2, Eq. (45) becomes:

[TABLE]

where Eq. (17) was used to eliminate the time derivative of $\beta$ in Eq. (46a), while the relaxation times were defined as $\tau_{q}=\lambda(T\beta_{1}-1/4n_{0})$ and $\tau_{\Pi}=2\eta\beta_{2}$. In order for the constitutive equation for $q$ and $\Pi$ to match the first-order versions (27) in the limit $\tau_{q},\tau_{\Pi}\rightarrow 0$, we set $\alpha_{1}=-1/4P_{0}$ and Eqs. (46) reduce to:

[TABLE]

By analogy to Eq. (25), the reduced relaxation times $\tau_{q,0}$ and $\tau_{\Pi,0}$ can be introduced through:

[TABLE]

In general, the values of $\tau_{q,0}$ and $\tau_{\Pi,0}$ are determined by the properties of the collision term in the Boltzmann equation. When the Anderson-Witting single relaxation time approximation is used, the following values for $\tau_{q,0}$ and $\tau_{\Pi,0}$ are commonly employed within both the Chapman-Enskog- and moment-like methods denicol12 ; jaiswal13a ; bhalerao14 ; chattopadhyay15 ; tinti17 :

[TABLE]

It is interesting to test the hyperbolicity of the resulting set of equations. Equations (17) and (47) can be written in the following form:

[TABLE]

where

[TABLE]

while $\mathbb{A}$ is given by:

[TABLE]

The five eigenvalues of $\mathbb{A}$ can be found analytically and are given by $a_{0}=0$ and

[TABLE]

Since the eigenvalues of $\mathbb{A}$ are real, the system of equations (50) is hyperbolic rezzolla13 ; toro09 .

IV.2 Longitudinal waves: modes

Employing the ansätze (20) and (21) allows $\Pi_{\alpha}$ to be expressed from Eq. (47b) as follows:

[TABLE]

Furthermore, Eq. (22b) allows $\delta P_{\alpha}$ to be written as:

[TABLE]

Inserting Eq. (54) and (55) into Eq. (22c) yields:

[TABLE]

The above equation is satisfied when either of the two factors cancel. These two cases are discussed separately below.

$q_{\alpha}=-4P_{0}\beta_{\alpha}$.

This case is also recovered for $q_{\lambda}$ in the first-order approximation (32). In this case, Eqs. (54) and (55) show that $\Pi_{\alpha}=\delta P_{\alpha}=0$. The values of $\alpha$ corresponding to this regime can be found from Eq. (47a), which reduces to:

[TABLE]

The case $\beta_{\alpha}=0$ trivially corresponds to a vanishing perturbation (i.e. $\delta n_{\alpha}=\delta P_{\alpha}=0$). Setting the term inside the square bracket to [math] yields the following values for $\alpha$:

[TABLE]

It can be seen that, in the small $\tau$ limit, $\alpha_{\lambda,-}$ corresponds to the dampening coefficient $\alpha_{\lambda}=k^{2}\lambda/4n_{0}$ identified in the first-order theory (31). The term $\alpha_{\lambda,+}$ allows $\widetilde{q}$ to relax from an arbitrary initial condition at $t=0$ to the first-order expression (27) on a timescale of order $\tau_{q}$. The modes corresponding to the above values of $\alpha_{\lambda,\pm}$ can be written in terms of the amplitudes $\beta_{\lambda,\pm}$ as follows:

[TABLE]

The square root in the definition of $\alpha_{\lambda,\pm}$ (58) becomes imaginary when $\tau>\tau_{\lambda,{\rm lim}}$, where

[TABLE]

When $0<\tau<\tau_{\lambda,{\rm lim}}$, the modes corresponding to $\alpha_{\lambda,\pm}$ are overdamped (evanescent), while when $\tau>\tau_{\lambda,{\rm lim}}$, the underdamped regime settles in. In order to treat both regimes within a unitary framework, it is convenient to cast Eq. (58) in the form $\alpha_{\lambda,\pm}=\alpha_{\lambda,d}\pm\alpha_{\lambda,o}=\alpha_{\lambda,d}\pm i\overline{\alpha}_{\lambda,o}$, where $\alpha_{\lambda,d}=1/2\tau_{q}$ and

[TABLE]

$q_{\alpha}=0$.

Next, the case when the first bracket in Eq. (56) vanishes is considered. Substituting Eqs. (22a) and (55) in Eq. (47a) yields:

[TABLE]

The term inside the square bracket is identical to the one in Eq. (57). Setting this term to [math] yields the same values for $\alpha$ as considered previously, when $q_{\alpha}=-4P_{0}\beta_{\alpha}$. To obtain a different set of values for $\alpha$, Eq. (62) is now solved by setting $q_{\alpha}=0$.

The allowed values for $\alpha$ can be found by solving the following cubic equation:

[TABLE]

Equation (63) has the roots $\alpha\in\{\alpha_{\eta,r},\alpha_{\eta,\pm}\}$, where the notation $\alpha_{\eta,\pm}=\alpha_{\eta,d}\pm i\alpha_{\eta,o}$ is introduced by analogy to the first-order case (34). The exact expressions for the coefficients $\alpha_{\eta,r}$, $\alpha_{\eta,d}$ and $\alpha_{\eta,o}$ are:

[TABLE]

where

[TABLE]

The parameter $\tau_{\eta,\rm lim}$ is defined as the value of $\tau$ at which the expression under the cubic root in Eq. (65) vanishes. It is given by:

[TABLE]

The definition (65) of $R_{\eta}$ ensures that the coefficients $\alpha_{\eta,*}$ ($*\in\{r,d,o\}$), defined in Eq. (64), are real for all positive values of $\tau$, provided that $R_{\eta,{\rm aux}}>0$. In order to investigate the latter inequality, it is convenient to consider $\eta/\tau_{\Pi}P_{0}$ as an independent parameter. In this case, the roots of $R_{\eta,{\rm aux}}$ (66) are:

[TABLE]

It can be seen that $(k^{2}\tau_{\Pi}^{2})_{\pm}$ is real only when $\eta/\tau_{\Pi}P_{0}\geq 8$. In the hydrodynamic regime, the Chapman-Enskog expansion (26b) together with Eq. (49) predict that $\eta/\tau_{\Pi}P_{0}=4/5$, which is much smaller than $8$. For the sake of simplicity, the case when $(k^{2}\tau_{\Pi}^{2})_{\pm}$ are real will not be consider in this paper. Instead, the coefficients $\alpha_{\eta,*}$ (64) will be assumed to be real for all values of $\tau$ considered in this section. It is worth mentioning that setting $\eta_{0}=4/5$, $\tau_{\Pi,0}=1$ and $k=2\pi$ in Eq. (67) yields $\tau_{\eta,{\rm lim}}\simeq 0.12$, which is within the range of values of $\tau$ considered in Subsec. IV.4.

The small $\tau$ expansion of Eq. (64) reveals the first-order coefficients $\alpha_{d}=k^{2}\eta/6P_{0}$ and $\alpha_{o}=k/\sqrt{3}+O(\tau^{2})$ reported in Eq. (35). The coefficient $\alpha_{\eta,r}$ was not present in the first-order theory and in this case it corresponds to the mode that ensures the relaxation of $\widetilde{\Pi}$ from $\widetilde{\Pi}(t=0)=0$ to the value predicted through the constitutive equation of the first-order theory.

IV.3 Longitudinal waves: solution

The full solution of the longitudinal wave problem can be written in general in the following form:

[TABLE]

where $\widetilde{M}\in\{\widetilde{\beta},\widetilde{\delta n},\widetilde{\delta P},\widetilde{q},\widetilde{\Pi}\}$. The subscripts $\lambda$ and $\eta$ refer to the parts of the solution corresponding to $\alpha_{\lambda,\pm}$ and $\alpha_{\eta,*}$ ($*\in\{r,d,o\}$), respectively.

According to Eq. (59), $\widetilde{\delta P}_{\lambda}=\widetilde{\Pi}_{\lambda}=0$, while $\widetilde{\beta}_{\lambda}$, $\widetilde{\delta n}_{\lambda}$ and $\widetilde{q}_{\lambda}$ can be written for $\tau<\tau_{\lambda,{\rm lim}}$ (60) as:

[TABLE]

When $\tau>\tau_{\lambda,{\rm lim}}$, the hyperbolic functions in the above expression become trigonometric functions:

[TABLE]

The constants $\beta_{\lambda,*}$ and $\overline{\beta}_{\lambda,*}$ (where $*\in\{c,s\}$) are fixed by the initial conditions, while:

[TABLE]

The part of the solution (69) corresponding to $\alpha_{\eta,*}$ ($*\in\{r,d,o\}$) can be written as:

[TABLE]

while $\widetilde{q}_{\eta}=0$. The integration constants $\beta_{\eta,*}$ are fixed by the initial conditions, while

[TABLE]

The initial conditions (23) for $\widetilde{\beta}$, $\widetilde{\delta n}$ and $\widetilde{\delta P}$ are supplemented by

[TABLE]

Substituting the solution (69) into the initial conditions equations (23) and (75) yields:

[TABLE]

The equation $q_{\lambda,c}=0$ implies that $\beta_{\lambda,c}=0$ (this is also true when $\tau_{q}>\tau_{\lambda,{\rm lim}}$, i.e. $\overline{q}_{\lambda,c}=\overline{\beta}_{\lambda,c}=0$). Furthermore, noting that $\delta P_{\eta,*}=4P_{0}\delta n_{\eta,*}/3n_{0}$, the second equality in Eq. (76) implies:

[TABLE]

The above analytic results will now be employed to study the relaxation process for $\widetilde{q}$ and $\widetilde{\Pi}$ in the context of Case 1 and Case 2b introduced in Sec. II.3.

IV.4 Case 1: Adiabatic flow

Setting $\delta n_{0}=\delta P_{0}=0$ in Eq. (77a) yields $\beta_{\lambda,s}=0$, such that $\widetilde{M}_{\lambda}=0$ for all $\widetilde{M}\in\{\widetilde{\beta},\widetilde{\delta n},\widetilde{\delta P},\widetilde{q},\widetilde{\Pi}\}$. The integration constants $\beta_{\eta,r}$, $\beta_{\eta,c}$ and $\beta_{\eta,s}$ can be found from Eq. (77b) as follows:

[TABLE]

Substituting the above expressions in Eq. (74) yields:

[TABLE]

Figure 6(a) shows the numerical results for the evolution of $\widetilde{\Pi}$ for $\beta_{0}=10^{-3}$ and $\tau=0.056$, compared with the first- and second-order hydrodynamics solutions given in Eqs. (43) and (79), respectively. The early time disagreement between the first-order hydrodynamics and numerical results becomes negligible when $t\gtrsim 2.5$, while the second-order solution is in excellent agreement with the numerical results at all values of $t$. At $\tau=0.1$, Fig. 6(b) shows that the first-order hydrodynamics prediction remains in visible disagreement at large times, while the initial disagreement between the second-order solution and the numerical results becomes negligible for $t\gtrsim 2$.

The early-time validity of the solution (79) is inspected in Fig. 6(c). At small values of $\tau$, $\widetilde{\Pi}$ relaxes from its initial vanishing value to the value predicted by the first-order expression (43) after a time $t\sim 5\tau$. For $\tau\gtrsim 0.05$, the first-order approximation becomes non-satisfactory at small values of $t$, since it lags behind the numerical solution. At small $\tau$, the second-order approximation is overlapped with the numerical solution for all values of $t$. An early-time discrepancy between the second-order prediction and the numerical result can be seen when $\tau\gtrsim 0.1$, which however becomes negligible at large times, as shown in Fig. 6(b). At sufficiently small values for $\tau$, Fig. 6(d) shows that this discrepancy arises when the values of $\eta_{0}$ and $\tau_{\Pi,0}$ predicted by Eqs. (26b) and (49) are used in the analytic expression. If instead, these values are treated as free parameters, the functional form (79) can be used to accurately represent $\widetilde{\Pi}$, even at $\tau\simeq 0.26$.

Next, the validity of Eqs. (26b) and (49) pertaining to the Chapman-Enskog expressions for $\eta_{0}$ and $\tau_{\Pi,0}$ is examined. For this purpose, a nonlinear fit of the functional forms in Eq. (79) will be performed, where the coefficients $\alpha_{\eta,r}$, $\alpha_{\eta,d}$ and $\alpha_{\eta,o}$ are considered free parameters. In addition, the solutions for $\widetilde{\delta n}$ and $\widetilde{\delta P}$ depend explicitly on $\tau_{\Pi}$, while $\widetilde{\Pi}$ depends explicitly on $\eta/\tau_{\Pi}$. The inversion of Eq. (64) would allow $\eta$ and $\tau_{\Pi}$ to be written in terms of $\alpha_{\eta,r}$ and $\alpha_{\eta,d}$, but this operation is mathematically intractable. Thus, $\tau_{\Pi}$ will also be treated as a free parameter for $\widetilde{\delta n}$ and $\widetilde{\delta P}$, while in the case of $\widetilde{\Pi}$, $\eta/\tau_{\Pi}$ will be considered as an independent parameter. The results of the numerical fits for $\alpha_{\eta,r}$, $\alpha_{\eta,d}$ and $\alpha_{\eta,o}$ are presented in Fig. 7. In the case of the $\alpha_{\eta,r}$ coefficient, it can be seen that the analytic expression predicts a sharper decrease at large $\tau$ compared to the numerical results. The shapes of $\alpha_{\eta,d}$ and $\alpha_{\eta,o}$ remain qualitatively similar to those obtained for $\alpha_{d}$ and $\alpha_{o}$ in the first-order theory, which are shown in Fig. 5. In the first-order theory, $\alpha_{d}$ is directly proportional to $\tau$, while the numerical results seem to indicate a saturation of $\alpha_{d}$ for $\tau\gtrsim 0.1$. This saturation can be seen also for $\alpha_{\eta,d}$, but in this case, the analytic expression predicts that $\alpha_{\eta,d}$ decreases with $\tau$. Another notable difference can be seen in the analytic prediction for $\alpha_{\eta,o}$, which in the second order case qualitatively follows the numerical results ($\alpha_{\eta,o}$ increases at large $\tau$) compared to the first-order case, when $\alpha_{o}$ is predicted to decrease at large $\tau$.

Finally, the analysis of the dependence of the relaxation time $\tau_{\Pi}$ and of the ratio $\eta/\tau_{\Pi}P_{0}$ on $\tau$ is presented below. The curves shown in Fig. 8 represent three types of results. The first type corresponds to the various analytic predictions, which are represented as follows: in Fig. 8(a), $\tau_{\Pi}$ (49) is shown using a continuous line; in Figs. 8(b) and 8(c), the Chapman-Enskog and Grad predictions for $\eta/P_{0}\tau_{\Pi}$ and $\eta$ are shown using continuous and dotted lines. The second type of results are obtained using the nonlinear fit procedure described in the previous paragraph for $\tau_{\Pi}$ (obtained from $\widetilde{\beta}$, $\widetilde{\delta n}$ and $\widetilde{\delta P}$) and for $\eta/\tau_{\Pi}P_{0}$ (obtained from $\widetilde{\Pi}$), being labelled using the suffix “-fit”. The third type of results are obtained by numerically finding the values of $\tau_{\Pi}$ and $\eta/\tau_{\Pi}P_{0}$ for which the roots of Eq. (63) correspond to the values of $\alpha_{\eta,r}$ and $\alpha_{\eta,d}$ obtained through the nonlinear fit procedure (the value of $\alpha_{\eta,o}$ is not taken into account in this algorithm). As expected, the relaxation time $\tau_{\Pi}$ stops increasing linearly with $\tau$ when $\tau\gtrsim 0.1$. Figure 8(b) shows that the ratio $\eta/\tau_{\Pi}P_{0}$ increases with $\tau$, which seems to indicate that the increase of the shear viscosity of the gas with respect to $\tau$ is of higher order than the linear prediction of the first-order theory (24). However, Fig. 8(c) shows that in fact $\eta$ (obtained by multiplying the results for $\eta/\tau_{\Pi}P_{0}$ and $\tau_{\Pi}$ obtained as explained above) increases at a sub-linear rate with respect to $\tau$ when $\tau\gtrsim 0.1$. This result is consistent with the one obtained in Fig. 5 in the first-order case. When $\tau\lesssim 0.1$, the numerical results favor the Chapman-Enskog prediction for the transport coefficients, as well as the relations (49).

IV.5 Non-adiabatic flow

The non-adiabatic case can be analysed when the system is initialized according to Case 2b, when $\beta_{0}=\delta P_{0}=0$. This case is particularly simple since, according to Eq. (77b), $\beta_{\eta,*}=0$ for all $\eta\in\{r,c,s\}$. The only non-vanishing integration constants are $\beta_{\lambda,s}$ (overdamped case) and $\overline{\beta}_{\lambda,s}$ (underdamped case), which can be found from Eq. (77a):

[TABLE]

Noting that $\widetilde{\delta P}=\widetilde{\Pi}=0$ by virtue of Eq. (73), the full solution in the overdamped case reads:

[TABLE]

In the first-order theory, Eqs. (40) and (38) predict that, when $\beta_{0}=\delta P_{0}=0$, $\beta_{c}$ and $\beta_{s}$ become proportional to $\delta n_{0}$, which also implies that $\widetilde{\Pi}$ and $\widetilde{\delta P}$ are non-zero. This prediction is in contradiction with the second-order theory results. This discrepancy holds for any value of $\tau$, hence it cannot be considered a “higher order” (rarefaction) effect. Instead, it represents a fundamental flaw of the first-order theory, which can be explained as follows.

The solution $\widetilde{\Pi}=\widetilde{\delta P}=0$ is supported in the first-order theory only by the $\alpha_{\lambda}$ mode, as indicated in Eq. (32). However, according to Eq. (32), $\delta n_{\lambda}$ is proportional to $\beta_{\lambda}$, such that requiring that $\beta_{\lambda}=0$ automatically implies $\delta n_{\lambda}=0$. Thus, the initial conditions corresponding to Case 2b cannot be imposed using only the $\alpha_{\lambda}$ mode, such that the $\alpha_{\pm}$ modes also become excited. The existence of two modes $\alpha_{\lambda,\pm}$ in the second-order theory which support $\widetilde{\Pi}=\widetilde{\delta P}=0$ (59) is sufficient to allow the initial conditions of Case 2b to be imposed without exciting the $\alpha_{\eta,*}$ modes, such that $\widetilde{\Pi}=\widetilde{P}=0$ throughout the evolution.

The case discussed above is represented in the plots on the first line of Fig. 9, where the time evolution of $\widetilde{\beta}$, $\widetilde{\delta P}$ and $\widetilde{\Pi}$ (normalised with respect to $\delta n_{0}$) is represented for the initial conditions $\delta n_{0}=10^{-3}$, $\delta P_{0}=0$ and $\beta_{0}=0$. Both the second-order theory and the numerical results indicate that $\widetilde{\delta P}$ and $\widetilde{\Pi}$ remain zero throughout the evolution, while the first-order result predicts oscillations of these quantities. Moreover, in the first-order theory, $\widetilde{\beta}$ presents strong oscillations about a decaying exponential, which are not present when the second-order theory is employed. It is worth mentioning that the evolution predicted by the second-order theory and captured by the numerical method is consistent with the evolution equations (17), which reduce for $\delta P=\Pi=0$ to:

[TABLE]

The above relations are also recovered for the $\alpha_{\lambda}$ mode (32) of the first order theory. However, the initial conditions $\beta_{0}=\delta P_{0}/P_{0}=0$ and $\delta n_{0}=10^{-3}$ cannot be imposed using only the $\alpha_{\lambda}$ mode, since according to Eq. (32), $\delta n_{\lambda}=kn_{0}\beta_{\lambda}/\alpha_{\lambda}$ implies that $\delta n_{0}$ is proportional to $\beta_{\lambda}$. Thus, setting $\beta_{0}=0$ and $\delta n_{0}\neq 0$ automatically excites the modes $\delta n_{\pm}$ and $\beta_{\pm}$, which no longer satisfy Eq. (82), such that the solution becomes contaminated with the addition of oscillatory modes.

The discussion in the preceding paragraph suggests that the first-order theory cannot fully recover the hydrodynamic regime (small $\tau$) of the decaying longitudinal wave when the initialisation is performed according to Case 2b, i.e. $\beta_{0}=\delta P_{0}=0$ and $\delta n_{0}=10^{-3}$. It is instructive to further test if this conclusion holds at non-vanishing values of $\beta_{0}$, while keeping $\delta P_{0}=0$. The second and third lines of Fig. 9 show the evolution of $\widetilde{\beta}$, $\widetilde{\delta P}$ and $\widetilde{\Pi}$ (again normalised with respect to $\delta n_{0}$) when $\beta_{0}=10^{-5}$ and $\beta_{0}=10^{-3}$, respectively. At $\beta_{0}=10^{-5}$, the first-order theory predicts oscillations which are (nearly) in antiphase to the numerical and second-order theory results. When $\beta_{0}=10^{-3}$, the discrepancy between the curves corresponding to the first-order theory and the numerical and second-order theory results is no longer visible.

The ability of the second-order hydrodynamics to capture the relaxation of $\widetilde{q}$ from [math] at initial time to the value predicted by the first-order theory is investigated in Fig. 10(a). It can be seen that at small $\tau$, $\widetilde{q}$ relaxes to the value predicted through the first-order theory with $\lambda$ given by the Chapman-Enskog expansion (26b) after a time $t\sim 5\tau$. At $\tau\gtrsim 0.1$, the first-order approximation seems to no longer agree with the numerical solution (as indicated in Sec. III.6), while the second-order approximation slowly loses its validity. In Fig. 10(b), the numerical result for the evolution of $\widetilde{q}$ at $\tau=0.1$ is compared to the first (44) and second (81a) order hydrodynamics predictions, specialised to the Chapman-Enskog case, when $\lambda_{0}=4/3$ and $\tau_{q,0}=1$. It can be seen that there is significant discrepancy between the analytic and numerical results. The curve labelled Hydro 1 (fit) represents the best fit of the functional form of the analytic solution (44) to the numerical results, with $\alpha_{\lambda}$ considered as a free parameter. The second-order fits are performed on the two functional forms (81a) and (81b), corresponding to the overdamped and underdamped cases, respectively. In the overdamped case, Eq. (81a) is fitted to the numerical data by considering $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}$ as free parameters. In the underdamped case, $\overline{\alpha}_{\lambda,d}$ and $\overline{\alpha}_{\lambda,o}$ are found by performing a nonlinear fit of Eq. (81b) with respect to the numerical data. It can be seen in Fig. 10(b) that, at $\tau=0.1$, the first- and underdamped second-order fits still present significant deviations from the numerical curve. However, the second-order overdamped fit is in very good agreement with the numerical result, validating the functional form (81a) for $\tau\lesssim 0.1$.

Next, Fig. 11 shows an analysis of $(\alpha_{\lambda,d},\alpha_{\lambda,o})$ and $(\overline{\alpha}_{\lambda,d},\overline{\alpha}_{\lambda,o})$ as obtained by performing a two-parameter nonlinear fit of Eqs. (81a) and (81b) for the overdamped (OD) and underdamped (UD) regimes, respectively, to the numerical results obtained for $\widetilde{\beta}$, $\widetilde{\delta n}$ and $\widetilde{q}$. Only the initialisation corresponding to Case 2b (i.e. $\beta_{0}=\delta P_{0}=0$ and $\delta n_{0}=10^{-3}$) is considered here. The numerical fit confirms that at small values of $\tau$, $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}\simeq 1/2\tau$. Furthermore, the results shown in Fig. 11(a) for $\alpha_{\lambda}$ indicate that the UD regime is not valid when $\tau\lesssim 0.2$, while the OD regime loses applicability when $\tau\gtrsim 0.1$. This can also be seen in Fig. 11(b), where the strong spikes indicate the points where $\alpha_{\lambda,o}=0$, i.e. where the transition from the UD to the OD regime occurs. According to Eq. (60), this happens when $\tau=\tau_{\lambda,{\rm lim}}^{\rm CE}\simeq 0.138$ and $\tau=\tau_{\lambda,{\rm lim}}^{\rm G}\simeq 0.1779$ when the Chapman-Enskog (26b) and Grad (26a) values for $\lambda_{0}$ are used, respectively. The above numerical analysis indicates that $\tau_{\lambda,{\rm lim}}\simeq 0.199$ (in the case of $\widetilde{\delta n}$) and $\tau_{\lambda,{\rm lim}}\simeq 0.183$ (in the case of $\widetilde{\beta}$ and $\widetilde{q}$), higher than both the Chapman-Enskog and the Grad predictions.

The analysis of the relaxation time $\tau_{q}$ and the heat conductivity $\lambda$ is presented in Fig. 12. Starting from the numerical fits of $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}$, where the overdamped (OD) and underdamped (UD) values are taken when $\tau<\tau_{\lambda,{\rm lim}}$ and $\tau>\tau_{\lambda,{\rm lim}}$, respectively, the values of $\lambda$ and $\tau_{q}$ are found by requiring that the roots $\alpha_{\lambda,\pm}$ (58) satisfy:

[TABLE]

where the absolute value is interpreted in the usual sense in the case when $\alpha_{\lambda,\pm}$ are complex numbers. The results for $\tau_{q}$ and $\lambda$ are shown in Figs. 12(a) and 12(b), respectively. It can be seen that the numerical results for $\tau_{q}$ agree with the theoretical prediction $\tau_{q}=\tau$ at small values of $\tau$, while for larger values of $\tau$, $\tau_{q}$ seems to reach a plateau. The curve representing the numerical results for $\lambda$ is practically overlapped with the Chapman-Enskog prediction (26b), while at larger values of $\tau$, $\lambda$ seems to reach a plateau, in good qualitative agreement with the predictions of the first-order theory presented in Fig. 5.

IV.6 Summary

In this section, a particular form of the second-order hydrodynamics equations was employed to study the attenuation of a longitudinal wave in the linearised regime. More precisely, the choice $\alpha_{1}=-1/4P_{0}$ was made for the coupling constant between the shear pressure $\Pi$ and heat flux $q$ in the theory presented in Refs. hiscock83 ; rezzolla13 , in order to ensure consistency with the first-order hydrodynamics analysis presented in Sec. III.

The main aim of this section was to confirm that the second-order hydrodynamics correctly describes the relaxation process of $\Pi$ and $q$ from essentially arbitrary initial values (both vanish at initial time in the analysis presented in this paper) to some non-zero value which agrees with the prediction of the first-order theory at small enough values of $\tau$. This is shown in Figs. 6 and 10, while the values of the relaxation times $\tau_{\Pi}$ and $\tau_{q}$ are analysed with respect to the relaxation time $\tau$ of the Anderson-Witting model in Figs. 8 and 12.

During this analysis, a fundamental limitation of the first-order theory was pointed out, which can be summarised as follows. If at initial time, the velocity and pressure perturbations vanish (i.e. $\beta_{0}=\delta P_{0}=0$), the pressure $\widetilde{\delta P}$ and shear stress $\widetilde{\Pi}$ perturbations remain zero at all later times, while the attenuation of $\widetilde{\beta}$ is purely evanescent (non-oscillatory). This result is obtained analytically in the second-order theory and is confirmed via numerical simulations in Fig. 9. On the other hand, the first-order theory predicts that $\widetilde{\delta P}$ and $\widetilde{\Pi}$ are proportional to $\delta n_{0}$ and $\widetilde{\beta}$ has an oscillatory component the amplitude of which is proportional to $\tau$. This leads to the conclusion that the first-order theory cannot correctly describe this particular dissipative process, even when the relaxation time is small.

V Moment method

The analysis presented in the previous sections provides indication that the correct expresions for the transport coefficients $\lambda$ and $\eta$ are obtained using the Chapman-Enskog procedure. Recently denicol12 , it was shown that one of the main drawbacks of the moment method as originally employed by Israel and Stewart israel79 is that the distribution function $f$ is expanded with respect to the nonorthogonal basis formed of powers of the particle four-momentum $p^{\mu}$. Truncating this series at a finite order discards an infinite number of terms of first-order with respect to the Knudsen number ${\rm Kn}$ (in the Anderson-Witting model, ${\rm Kn}\sim\tau$). The solution proposed in Ref. denicol12 was to expand $f$ in terms of irreducible tensors with respect to the particle momentum $p^{\mu}$. The proposed scheme recovers the expressions for the transport coefficients $\lambda$ and $\eta$ obtained using the Chapman-Enskog expansion.

In this section, a moment-based method similar to the one introduced in Ref. denicol12 is considered. Spherical coordinates $\{p,\theta,\varphi\}$ are employed in the momentum space and the distribution function $f$ is expanded with respect to the generalised Laguerre polynomials for the momentum magnitude $p$, the Legendre polynomials for $\xi=\cos\theta$ and the trigonometric basis $\{\cos m\varphi,\sin m\varphi\}$ for $\varphi$. Due to the symmetries of the system, $f$ can be considered independent of $\varphi$, such that only the $p$ and $\xi$ expansions will be discussed. Truncating the system at order $N_{L}=1$ and $N_{\xi}=2$ with respect to the Laguerre and Legendre polynomials, respectively, yields a system of six equations for the five hydrodynamic variables $\delta n$, $\beta$, $\delta P$, $q$ and $\Pi$, as well as a non-hydrodynamic variable. The importance of this sixth variable in establishing the symmetry between the shear stress and heat flux solutions is illustrated, such that the underdamped (UD) and overdamped (OD) regimes discussed in Sec. IV are represented unitarily in this new solution.

V.1 Constitutive relations

In this section, the longitudinal wave problem is again approached, but this time by employing a moment method. Instead of performing the standard Grad-like expansion of $f$ in terms of polynomials in $p^{\mu}$, $f$ is expanded following Refs. blaga17aip ; blaga17prc as follows:

[TABLE]

where $L^{(2)}_{\ell}(z)$ are the generalised Laguerre polynomial of type $2$ and order $\ell$, which satisfy the following orthogonality relation:

[TABLE]

Thus, the coefficients $\mathcal{F}_{\ell}$ can be obtained as follows:

[TABLE]

Multiplying the Boltzmann equation (18), valid only in the linearised regime of the longitudinal wave problem, by $p^{2}L_{\ell}^{(2)}(p/T_{0})$ and integrating over $p$ yields:

[TABLE]

where the coefficients $\mathcal{F}^{({\rm eq})}_{\ell}$ corresponding to the equilibrium distribution function $f^{(\mathrm{eq})}$ are defined by analogy to Eq. (86):

[TABLE]

In the absence of collisions [i.e. when neglecting the right hand side in Eq. (87)], each coefficient $\mathcal{F}_{\ell}$ evolves independently. Since $f^{(\mathrm{eq})}$ is constructed only in terms of $N^{\mu}$ and $T^{\mu\nu}$, which can be written entirely in terms of $\mathcal{F}_{\ell}$ with $\ell=0$ and $\ell=1$, the evolution of $N^{\mu}$ and $T^{\mu\nu}$ is fully determined by considering Eq. (87) only for $\ell=0$ and $\ell=1$ and neglecting all higher $\ell$ terms blaga17prc .

The coefficients $\mathcal{F}_{\ell}$ are further expanded with respect to $\xi$ using the complete set of Legendre polynomials $P_{s}(\xi)$:

[TABLE]

where the coefficients $\mathcal{F}_{\ell,s}$ depend only on $z$ and $t$ and are obtained using the orthogonality of the Legendre polynomials as follows:

[TABLE]

The coefficients $\mathcal{F}_{\ell,s}^{({\rm eq})}$ corresponding to $f^{(\mathrm{eq})}$ can be defined in a similar manner:

[TABLE]

The expansion coefficients $\mathcal{F}_{\ell,s}$ can be linked to $N^{\mu}$ and $T^{\mu\nu}$ as follows:

[TABLE]

The coefficients $\mathcal{F}_{0,2}$ and $\mathcal{F}_{1,2}$ on their own have no correspondent with respect to $N^{\mu}$ and $T^{\mu\nu}$. The equilibrium coefficients $\mathcal{F}_{\ell,s}^{({\rm eq})}$ ($\ell=0,1$ and $s=0,1,2$) can be found from Eq. (19):

[TABLE]

while $\mathcal{F}_{0,2}^{({\rm eq})}\simeq 0$ and $\mathcal{F}_{1,2}^{({\rm eq})}\simeq 0$.

Using the recurrence relation:

[TABLE]

Eq. (87) can be projected on the space of the Legendre polynomials as follows:

[TABLE]

The above procedure produces an infinite system of equations corresponding to various values of $(\ell,s)$, where knowledge of $\mathcal{F}_{\ell,s+1}$ is required in order to determine the evolution of $\mathcal{F}_{\ell,s}$. As also discussed in Ref. blaga17prc , the above system can be closed at an order $Q$ by imposing $\mathcal{F}_{\ell,Q}=0$. This procedure is intimately related to the numerical method employed in this paper (described in detail in Ref. blaga17prc and also summarised in Appendix B). In particular, $Q$ represents the quadrature order of the model and the resulting system of equations is guaranteed to be hyperbolic. Since only the study of $N^{\mu}$ and $T^{\mu\nu}$ is of interest in this section, only the case $Q=3$ will be considered henceforth, such that $\mathcal{F}_{0,3}=\mathcal{F}_{1,3}=0$. The resulting set of equations can be written as:

[TABLE]

The above system is closed. Substituting Eqs. (92) and (93) into Eq. (96), the conservation equations (17) can be obtained, together with the following constitutive equations:

[TABLE]

Comparing the above equations to the second-order hydrodynamics constitutive equations (47), it can be seen that the transport coefficients have the following expressions:

[TABLE]

The above relations confirm the Chapman-Enskog prediction (26b) for $\lambda$ and $\eta$ and agree with the second-order hydrodynamics values for $\tau_{q}$ and $\tau_{\Pi}$ given in Eq. (49). Furthermore, the constitutive equation (97a) contains an extra term compared to the second-order hydrodynamics version (47a). A simple power counting shows that this term is cubic in the relaxation time $\tau$, hence it cannot be present in the second-order hydroynamics theory.

V.2 Longitudinal waves: modes

Considering now the propagation of a wave with wave number $k$, the ansatz (20) can be applied to the new variables $\mathcal{F}_{0,2}$ and $\mathcal{F}_{1,2}$ as follows:

[TABLE]

The mode decomposition (21) can be applied to $\widetilde{\mathcal{F}}_{0,2}$ and $\widetilde{\mathcal{F}}_{1,2}$ as follows:

[TABLE]

Substituting the above expansions into Eqs. (96c) and (96f) gives:

[TABLE]

The moment method introduced in this section bears many similarities with the second-order hydrodynamics method discussed in Sec. IV. In particular, since the constitutive equations (97b) and (47b) for $\Pi$ are the same in the two theories, Eqs. (54), (55) and (56) remain unchanged. The latter equation again can be solved either by setting $q_{\alpha}=-4P_{0}\beta_{\alpha}$ or by setting the square bracket to [math]. In the latter case, the allowed values for $\alpha$, namely $\alpha_{\eta,r}$ and $\alpha_{\eta,\pm}$ can be written as in Sec. IV.2, being given in Eq. (64). For completeness, these expressions are reproduced below, specialised to the values of $\eta$, $\lambda$, $\tau_{q}$ and $\tau_{\Pi}$ given in Eq. (98):

[TABLE]

where $\alpha_{d}=2k^{2}\tau/15$ is the first-order coefficient given in Eq. (35) and $R_{\eta}$ (65) becomes:

[TABLE]

The function $R_{\eta,{\rm aux}}$ (66) reduces to:

[TABLE]

Since the roots $(k^{2}\tau^{2})_{\pm}=\frac{1}{9}(5\pm 10i)$ of $R_{\eta,{\rm aux}}$ have a non-vanishing imaginary part and $R_{\eta,{\rm aux}}(\tau=0)=1$, $R_{\eta,{\rm aux}}>0$ for all values of $\tau$. The threshold value $\tau_{\eta,{\rm lim}}$ appearing in Eq. (103) is

[TABLE]

which coincides with Eq. (67) when $\eta$ and $\tau_{\Pi}$ are replaced according to Eq. (98).

When $q_{\alpha}=-4P_{0}\beta_{\alpha}$, Eq. (57) is replaced by:

[TABLE]

The square bracket cancels when $\alpha\in\{\alpha_{\lambda,r},\alpha_{\lambda,\pm}\}$, where $\alpha_{\lambda,\pm}$ can be written as $\alpha_{\lambda,\pm}=\alpha_{\lambda,d}\pm i\alpha_{\lambda,o}$. The exact expressions for the coefficients $\alpha_{\lambda,r}$, $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}$ read:

[TABLE]

where $\alpha_{\lambda}=k^{2}\tau/3$ is defined in Eq. (31) and

[TABLE]

In the above, $R_{\lambda,{\rm aux}}$ is defined as:

[TABLE]

Since the roots $(k^{2}\tau^{2})_{\pm}=\frac{5}{72}(11\pm 3i\sqrt{15})$ of $R_{\lambda,{\rm aux}}$ are complex, $R_{\lambda,{\rm aux}}>0$ for all values of $\tau$. The parameter $\tau^{\rm mom}_{\lambda,{\rm lim}}$ is defined as the value of $\tau$ at which the expression under the cubic root in Eq. (108) vanishes. It is given by

[TABLE]

being identical to $\tau_{\eta,{\rm lim}}$ (105). The definition (108) of $R_{\lambda}$ ensures that the coefficients $\alpha_{\lambda,*}$ ($*\in\{r,d,o\}$), defined in Eq. (107), are real for all positive values of $\tau$.

V.3 Longitudinal waves: solution

The solution can be split as in Eq. (69), i.e. $\widetilde{M}=\widetilde{M}_{\lambda}+\widetilde{M}_{\eta}$. In this case, the $\lambda$ and $\eta$ sectors of the solution have symmetric expressions, i.e.:

[TABLE]

The coefficients $M_{\times,*}$ (where $\times\in\{\lambda,\eta\}$ and $*\in\{r,c,s\}$) for $M\in\{\delta n,\mathcal{F}_{0,2}\}$ are

[TABLE]

On the $\lambda$ sector, $\delta P_{\lambda,*}=\Pi_{\lambda,*}=0$, while

[TABLE]

On the $\eta$ sector, $q_{\eta,*}=0$ and

[TABLE]

The initial conditions (23) and (75) refer only to $\widetilde{\beta}$, $\widetilde{\delta n}$, $\widetilde{\delta P}$, $\widetilde{q}$ and $\widetilde{\Pi}$. In the moment approach considered in this section, the coefficients $\mathcal{F}_{0,2}$ and $\mathcal{F}_{1,2}$ are also free to evolve [in fact, they contribute only one degree of freedom, since $\Pi=2T_{0}(\mathcal{F}_{0,2}-\frac{1}{3}\mathcal{F}_{1,2})$ is taken as an indepedent variable]. Since at $t=0$, the system is initialised with the equilibrium distribution $f^{(\mathrm{eq})}$, the initial conditions for $\mathcal{F}_{0,2}$ and $\mathcal{F}_{1,2}$ can be read from Eq. (93):

[TABLE]

Imposing the initial conditions (23), (75) and (115) on the solution (111) yields the following solution for the integration constants $\beta_{\lambda,*}$:

[TABLE]

while $\beta_{\lambda,c}=-\beta_{\lambda,r}$. The coefficients $\beta_{\eta,*}$ are the same as in Eq. (77b), which were obtained in the frame of the second-order theory.

Since the integration constants $\beta_{\eta,*}$ obtained in the moment approach coincide with those obtained in the second-order theory, the analytic solutions for $\widetilde{\delta P}$ and $\widetilde{\Pi}$ are the same in these two approaches. Moreover, for the initial conditions corresponding to Case 1 (i.e. $\delta n_{0}=\delta P_{0}=0$), when $\beta_{\lambda,*}=0$, the full analytic solution is identical in the moment approach and in the second order theory, being given in Eq. (79). Case 1 will therefore not be analysed in this section. Instead, Case 2b will be analysed in the following subsection.

V.4 Numerical results (Case 2b)

Setting $\beta_{0}=\delta P_{0}=0$ in Eq. (116) yields:

[TABLE]

while $\beta_{\lambda,c}=-\beta_{\lambda,r}$. Noting that $\widetilde{\beta}_{\eta}=0$, the heat flux is given simply by $\widetilde{q}=-4P_{0}\widetilde{\beta}$, while

[TABLE]

The full solution can be written as:

[TABLE]

while $\widetilde{q}=-4P_{0}\widetilde{\beta}$ and $\widetilde{\delta P}=\widetilde{\Pi}=0$.

The functional form of $\widetilde{q}$ obtained using the moment method is more convenient to use compared to the one given in the second-order hydrodynamics case (81), since in the former case, there is no distinction between the overdamped and the underdamped regimes. In Fig. 13, the validity of the solution (119) for $\widetilde{q}$ is tested at $\tau=0.22$ and $\tau=1.0$. Since, according to Eq. (60), $\tau>\tau_{\lambda,{\rm lim}}\simeq 0.14$ (the Chapman-Enskog value $\lambda_{0}=4/3$ and $\tau_{q,0}=1$ were used), Eq. (81b) is used to represent the analytic solution obtained in the frame of the second-order hydrodynamics theory. At $\tau=0.22$, the solution corresponding to the moments method is much closer to the numerical result than the second-order hydrodynamics one. When $\tau=1.0$, both theories give solutions which deviate considerably from the numerical results. In this regime, the validity of the functional form of the analytic solutions discussed above can be further tested. In the second-order hydrodynamics case, a nonlinear fit of the solutions (81a) and (81b) is performed by considering the coefficients $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}$ ($\overline{\alpha}_{\lambda,d}$ and $\overline{\alpha}_{\lambda,o}$) as free parameters. Fig. 13 shows that, at $\tau=0.22$, the fit corresponding to the overdamped form (81a) is less accurate than the fit corresponding to the underdamped form (81b). At $\tau=1.0$, the fit corresponding to the UD form also starts to present visible deviations from the numerical result. In the moment method solution (119), the nonlinear fit is performed by considering $\tau$, $\alpha_{\lambda,r}$, $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}$ as free parameters. The resulting fit is in much better agreement with the numerical results.

Next, the $\tau$ dependence of the coefficients $\tau$, $\alpha_{\lambda,r}$, $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}$ as obtained by performing a nonlinear fit of Eq. (119) to the numerical data is considered. Figure 14(a) shows that $\tau$-fit [i.e., the best fit value for the parameter $\tau$ appearing in Eq. (119)] depends non-monotonically on $\tau$, i.e. it reaches a maximum value around $\tau\simeq 0.22$ and $\tau\simeq 0.32$ when considering the evolution of $\widetilde{\delta n}$ and $\widetilde{q}$, respectively, after which it decreases with $\tau$.

The analytic expression for $\alpha_{\lambda,r}$ (107), represented using a continuous line in Fig. 14(b), reduces at small values of $\tau$ to $\alpha_{\lambda}$ (31) defined within the first-order theory. While the first-order theory predicts a linear increase of $\alpha_{\lambda}$ with $\tau$, the moment method predicts a maximum of $\alpha_{\lambda,r}$ at $\tau\simeq 1.19/k\simeq 0.19$, after which it decreases according to the asymptotic behaviour $\lim_{\tau\rightarrow\infty}\alpha_{\lambda,r}=5/9\tau$. The fitted values also are non-monotonic, exhibiting a slight decreasing trend when $\tau\gtrsim 0.38$ in the case of $\widetilde{q}$ and $\widetilde{\beta}$ and $\tau\gtrsim 0.68$ when the nonlinear fit is performed on $\widetilde{\delta n}$.

The coefficient $\alpha_{\lambda,d}$ is analysed in Fig. 14(c). The analytic expression (107) predicts a monotonic decrease of $\alpha_{\lambda,d}$. The curve corresponding to the nonlinear fit also decreases monotonically over the range of $\tau$ considered in Fig. 14(c), but at a lesser rate compared to the analytic prediction.

Finally, the oscillation frequency $\alpha_{\lambda,o}$ (which has no analogue in the first-order theory) is represented in Fig. 14(d). It can be seen that both the analytic and the numerical curves indicate that $\alpha_{\lambda,o}$ starts to increase when $\tau\gtrsim 0.1$.

V.5 Summary

In this section, a moment-based method was employed to study the attenuation of a longitudinal wave. The moment equations were obtained by projecting the distribution function $f$ on the space of the generalised Laguerre and Legendre polynomials. The same expressions for the transport coefficients as those obtained through the Chapman-Enskog expansion are found. The difference between the approach taken in this section and the traditional Grad moment method introduced by Israel and Stewart israel79 is that the truncation of $f$ is performed with respect to orthogonal polynomials, while in the latter approach, a nonorthogonal polynomial basis is employed denicol12 . In this sense, the present approach is similar to that employed in Ref. denicol12 .

The minimal set of moment equations which gives access to the evolution of the macroscopic four-flow $N^{\mu}$ and stress-energy tensor $T^{\mu\nu}$ is obtained by retainig in the expansion of $f$ the zeroth and first order terms with respect to $p$ (expanded using generalised Laguerre polynomials) and zeroth, first and second order terms with respect to $\xi=p^{z}/p$ (expanded using Legendre polynomials). We note that an expansion with respect to $p$ and $\xi$ is also performed in denicol16 . The system contains 6 equations for the five hydrodynamic variables $\delta n$, $\beta$, $\delta P$, $q$ and $\Pi$, as well as for a non-hydrodynamic variable not present in the second-order hydrodynamics theory.

The solution of the set of moment equations is identical on the shear stress sector with the one obtained from the second order hydrodynamics equations discussed in Sec. IV. On the heat flux sector, the second-order hydrodynamics solution is improved in the moment method approach, where the functional form of $q$ allows for a smooth transition from the overdamped to the underdamped regimes highlighted in Sec. IV.5. Furthermore, the range of validity of the analytical solution of the moment equations is larger than the one corresponding to the second-order hydrodynamics equations. Moreover, the functional form of the former can be fitted to the numerical data with remarkable accuracy even at $\tau=1$.

VI The ballistic limit

This section ends the analysis of the longitudinal wave problem by considering the free-streaming limit. In this case, the relativistic Boltzmann equation (1) reduces to:

[TABLE]

where $\xi=p^{z}/p$. The solution of Eq. (120) is $f(z,\xi,t)=f(z-\xi t)$, subject to the following initial condition:

[TABLE]

In the case of the longitudinal wave, the initial conditions for the macroscopic fields are:

[TABLE]

In the above, $k=2\pi/L$ represents the wave number for a longitudinal wave having the wavelength equal to $L$. Assuming that $\delta n_{0}$, $\delta P_{0}$ and $\beta_{0}$ are small, Eq. (121) can be linearised as follows:

[TABLE]

The solution at $t>0$ is given by replacing the product $kz$ in Eq. (123) by $k(z-\xi t)$.

The time evolution of the macroscopic quantities $n$, $P$, $\beta$, $q$ and $\Pi$ can be obtained from $N^{\mu}$ and $T^{\mu\nu}$ (3), which reduce to:

[TABLE]

where $v^{\mu}=p^{\mu}/p=(1,\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$ and $\cos\theta=\xi$. The result is:

[TABLE]

The leading order term in all of the above expressions is damped according to a factor of $t^{-1}$.

Figure 15 illustrates the close agreement between the numerical results and the analytic solution (125) when $k=2\pi/L$ and $L=1$. Each plot in Fig. 15 contains three pairs of curves, each pair corresponding to the initial conditions described in Cases 1 ($\delta n_{0}=\delta P_{0}=0$, $\beta_{0}=10^{-3}$), 2a ($\delta n_{0}=\beta_{0}=0$, $\delta P_{0}=10^{-3}$) and 2b ($\delta P_{0}=\beta_{0}=0$, $\delta n_{0}=10^{-3}$). The plots illustrate the time evolution of $\widetilde{\delta n}$, $\widetilde{\delta P}$, $\widetilde{\beta}$, $\widetilde{q}$ and $\widetilde{\Pi}$, where the numerical results are represented with dashed lines and points, while the analytic expressions (125) are represented using solid lines. The quantities on the vertical axis are divided by the amplitude of the perturbation, namely $\beta_{0}$ for Case 1, $\delta P_{0}$ for Case 2a and $\delta n_{0}$ for Case 2b. For simplicity, the amplitude was taken equal to $10^{-3}$ in all cases. It can be seen that the agreement between the numerical results and the analytic expressions is excellent.

Each plot in Fig. 15 displays two non-trivial curves and a line corresponding to a vanishing value. This is because all of the expressions in Eq. (125) have on the right hand side only two terms, e.g. $\widetilde{\delta n}$ vanishes when $\delta n_{0}=\beta_{0}=0$ for all values of $\delta P_{0}$, etc. It is worth pointing out that $\widetilde{\delta P}$ and $\widetilde{\Pi}$ vanish when $\beta_{0}=\delta P_{0}/P_{0}=0$, as also predicted by the second-order hydrodynamics theory discussed in Sec. IV, as well as by the moment method presented in Sec. V.

A fundamental difference between the hydrodynamic and the free-streaming regimes is that the attenuation of the wave perturbation in the former case is exponential (dissipative), while in the latter case, it is of the form $t^{-1}$ (dispersive).

VII Conclusion and outlook

In this paper, the attenuation of a longitudinal wave in a medium formed of ultrarelativistic (massless) particles was studied from the following perspectives: the first- and second-order hydrodynamics equations, the moment method, the free-streaming regime and by employing the numerical method introduced in Ref. blaga17prc . These investigations were carried out by considering the linearised limit of the hydrodynamics equations, which can be solved analytically. The analytic solutions were confronted with the numerical results in order to highlight the properties of the transport coefficients (and relaxation times in the second-order hydrodynamics and moment-method cases) in this system, for three particular cases: in Case 1, the initial density and pressure perturbations vanish ($\delta n_{0}/n_{0}=\delta P_{0}/P_{0}=0$); in Case 2a, the initial density and velocity perturbations vanish ($\delta n_{0}/n_{0}=\beta_{0}=0$); finally, in Case 2b, the initial pressure and velocity perturbations vanish ($\delta P_{0}/P_{0}=\beta_{0}=0$). Since in Case 1, the flow is adiabatic (i.e. the heat flux vanishes at all times) for all tested values of the initial wave perturbation $\beta_{0}$, this case was considered for the study of the shear viscosity $\eta$ by following the attenuation of the density, velocity, pressure and shear pressure perturbations. Cases 2a and 2b were considered in order to study the heat conductivity $\lambda$ by following the attenuation of the heat flux.

Throughout this paper, two types of tests were performed: (a) comparisons of the time evolution of the amplitudes $\widetilde{\delta n}$, $\widetilde{\delta P}$, $\widetilde{\beta}$, $\widetilde{q}$ and $\widetilde{\Pi}$ of the wave perturbations obtained from the numerical simulations with the analytic predictions corresponding to the Chapman-Enskog and Grad expressions for the transport coefficients; and (b), nonlinear numerical fits of these analytic expressions to the numerical results by considering the transport coefficients as fitting parameters. The analytic solutions for the evolution of the amplitudes were considered in terms of the modes allowed by the hydrodynamic equations.

In the first-order theory, three modes were highlighted: an evanescent mode (corresponding to the dampening coefficient $\alpha_{\lambda}$) and two modes undergoing oscillatory attenuation described by the dampening coefficient $\alpha_{d}$ and the oscillation angular frequency $\alpha_{o}$. Since there is no contribution to the heat flux from the oscillatory modes, the evanescent mode can be regarded as describing the heat flux sector (also, $\alpha_{\lambda}$ is determined exclusively in terms of the heat conductivity $\lambda$). Conversely, the oscillatory modes describe the shear pressure sector, since $\alpha_{d}$ and $\alpha_{o}$ depend only on the shear viscosity $\eta$. Thus, the numerical fits of $\widetilde{\delta n}$, $\widetilde{\beta}$, $\widetilde{\delta P}$ and $\widetilde{\Pi}$ were performed by considering $\alpha_{d}$ and $\alpha_{o}$ as free parameters, while $\alpha_{\lambda}$ was considered as a free parameter during the nonlinear fit of the expression for $\widetilde{q}$ to the numerical data.

In the analysis based on the second-order hydrodynamics equations, the heat flux sector is described by two modes ($\alpha_{\lambda,+}$ and $\alpha_{\lambda,-}$) which are evanescent for $\tau$ smaller than some value $\tau_{\lambda,{\rm lim}}$, while for $\tau>\tau_{\lambda,{\rm lim}}$, their attenuation is oscillatory. The coefficients $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}$ now depend on $\lambda$ and also on the heat flux relaxation time $\tau_{q}$. The shear pressure sector is also enlarged by the addition of an evanescent mode corresponding to the dampening coefficient $\alpha_{\eta,r}$, while the other two modes describe an oscillatory attenuation with dampening coefficient $\alpha_{\eta,d}$ and angular frequency $\alpha_{\eta,o}$. These three coefficients depend on $\eta$, as well as on the shear pressure relaxation time $\tau_{\Pi}$. The nonlinear fit of the heat flux was performed by considering $\alpha_{\lambda,d}$ and $\alpha_{\lambda,o}$ as free parameters with the expression of $\widetilde{q}$ written in both the evanescent (overdamped, OD) and in the oscillatory (underdamped, UD) forms. The nonlinear fits of the other amplitudes was performed by considering $\alpha_{\eta,*}$ ($*\in\{r,d,o\}$) as free parameters, as well as a fourth parameter ($\tau_{\Pi}$ for $\widetilde{\delta n}$, $\widetilde{\delta P}$ and $\widetilde{\beta}$ and the ratio $\eta/\tau_{\Pi}$ for $\widetilde{\Pi}$), which was considered as a free parameter due to the mathematical form of the analytic solution.

In the case of the moment-based method, the shear pressure sector was found to be identical to that obtained within the second-order hydrodynamics approach. Due to the addition of a sixth (non-hydrodynamic) mode, the heat flux sector was enlarged by the addition of a purely evanescent mode damped by the coefficient $\alpha_{\lambda,r}$, while the other two modes are of oscillatory type, damped by the coefficient $\alpha_{\lambda,d}$ and having oscillation frequency $\alpha_{\lambda,o}$. In comparison to the second-order hydrodynamics result, this solution behaves as nearly evanescent at small values of $\tau$, since $\alpha_{\lambda,d}\sim\tau^{-1}$ quickly suppresses the oscillatory contributions. At larger values of $\tau$, the oscillatory modes contribute significantly to the time evolution of $\widetilde{q}$, which explains the better agreement to the numerical data observed at large values of $\tau$. The nonlinear fit of $\widetilde{q}$ in the case of the moment method was performed by considering $\alpha_{\lambda,*}$ ($*\in\{r,d,o\}$) and $\tau$ as free parameters. The nonlinear fit of $\widetilde{\Pi}$ was performed by considering $\alpha_{\eta,*}$ ($*\in\{r,d,o\}$) and $\tau$ as free parameters.

Both tests described above support the conclusion that, at small values of the Anderson-Witting relaxation time $\tau$ (typically, $\tau\lesssim 0.05$), the expressions for the first-order transport coefficients ($\lambda$ and $\eta$) corresponding to the Anderson-Witting collision term are those predicted through the Chapman-Enskog procedure, which differ from the expressions obtained using Grad’s 14 moment approach. This conclusion is supported by various evidence in the literature florkowski13a ; florkowski13b ; florkowski15 ; ryblewski15 ; denicol14 ; bhalerao14 ; blaga17prc ; gabbana17b . Also in the limit $\tau\lesssim 0.05$, the above tests confirmed that the relaxtion times for the shear pressure $\tau_{\Pi}$ and heat flux $\tau_{q}$ are $\simeq\tau$. These results were also obtained analytically when a moment-based method was employed in order to construct the solution of the AWB equation. As remarked in Ref. denicol12 , this moment-based method is capable of reproducing the Chapman-Enskog transport coefficients since the distribution function is expanded with respect to orthogonal polynomials, while in the standard Grad method, the expansion is performed with respect to the nonorthogonal basis consisting of powers of the particle momentum $p^{\mu}$.

By performing the test (a), it was highlighted that the analytic solution obtained using the first-order hydrodynamics theory loses applicability when $\tau\gtrsim 0.05$. Since the constitutive equations for the heat flux and shear pressure tensor do not allow initial conditions to be specified for these fields, the first-order approximation is always inaccurate for a time scale $t\simeq 5\tau$. On this interval, the solution of the second-order hydrodynamics equations reproduce with good accuracy the numerical results for $\tau\lesssim 0.1$. While the evolution of the amplitude of the shear pressure $\widetilde{\Pi}$ is the same within the frames of the second-order hydrodynamics and the moment-based method considered in this paper, numerical experiments show that the evolution of the heat flux amplitude $\widetilde{q}$ is better captured by the moment method, which offers a reasonable agreement with the numerical results up to $\tau\simeq 0.22$.

Furthermore, the viability of the functional form of the analytic solutions obtained using the various hydrodynamic theories described above was considered. To this end, the nonlinear fits corresponding to test (b) were performed and the results were analysed in three ways, as described below.

First, a comparison was considered between the numerical results for $\widetilde{q}$ and $\widetilde{\Pi}$ and their analytic expressions corresponding to the best-fit values of the free parameters. At $\tau=0.26$, this analysis was used to highlight that the analytic expression for $\widetilde{\Pi}$ obtained using the second-order hydrodynamics theory (also from the moment-based method) was indistinguishable from the numerical results when the best-fit values of $\alpha_{\eta,*}$ and $\tau_{\Pi}$ were used, compared to the analytic prediction for these coefficients. The improvement of the analytic expression obtained using the first-order hydrodynamics equations was not significant, since this expression does not permit the value of $\widetilde{\Pi}$ to be fixed at $t=0$. In the case of the heat flux, the moment-based method provided a much more robust analytic expression for $\widetilde{q}$, which could be fitted remarkably well to the numerical results even at $\tau=1$, while the solution obtained within the second-order hydrodynamics formulation corresponding to the best fit parameters was in visible disagreement compared to the numerical result, although the overall evolution was still in reasonable agreement with the numerical data.

Furthermore, the dependence on $\tau$ of the free parameters used in the nonlinear fitting procedure was considered. In all approaches, the dampening coefficient of the oscillatory modes on the shear pressure sector ($\alpha_{\eta}$ in the first-order theory, $\alpha_{\eta,d}$ in the second-order theory and in the moment-based approach) was predicted analytically to grow (almost) linearly with $\tau$. However, the numerical fits indicate that this coefficient increases at a much slower rate when $\tau\gtrsim 0.05$. A similar behaviour was highlighted for the coefficient governing the evanescent mode on the heat flux sector ($\alpha_{\lambda}$ in the first-order theory, $\alpha_{\lambda,-}$ in the second-order theory and $\alpha_{\lambda,r}$ in the moment-based method). Thus, the above analysis indicates that the hydrodynamic theories considered in this paper break down when $\tau\gtrsim 0.05$.

A similar analysis of the dependence of the best fit parameters with respect to the wave amplitude at fixed $\tau=0.0083$ was performed. Significant deviations from the analytic predictions were found when the amplitude ($\beta_{0}$ for Case 1, $\delta P_{0}/P_{0}$ for Case 2a and $\delta n_{0}/n_{0}$ for Case 2b) exceeded $\sim 0.05$, indicating the inapplicability of the analysis in the linearised approximation at wave amplitudes larger than this value.

In order to gain some insight on the reason for the failure of the hydrodynamic theories to describe the attenuation of the longitudinal wave at larger values of $\tau$, the ballistic (free molecular flow) limit of this problem was investigated. In this case, the analytic solution for the distribution function $f$ indicates that the attenuation of the longitudinal wave is dispersive (i.e. the dampening is polynomial in $t^{-1}$) instead of dissipative (i.e. there is no exponential dampening). This behaviour was exactly recovered numerically, confirming the applicability of the numerical method in this regime. Thus, the solution of the hydrodynamics equations cannot describe correctly the attenuation of the longitudinal wave in the transition regime, where the dispersive component becomes important, since the functional form of these solutions does not include terms which are polynomial in $t^{-1}$.

It is worth presenting the particular case when at initial time, the velocity and pressure perturbations of the wave vanish, i.e. $\beta_{0}=\delta P_{0}/P_{0}=0$, while $\delta n_{0}=10^{-3}$. In this case, the second-order hydrodynamics theory and the moment method predict that the pressure and shear pressure remain constant in time and the attenuation of the density, velocity and heat flux is purely evanescent (i.e. non-oscillatory). This prediction is confirmed by the numerical simulations. However, the first-order theory always predicts an oscillatory attenuation of all variables (including the pressure and shear pressure, but excluding the heat flux), where the amplitude of the oscillations is of the same order of magnitude as the evanescent component. The above behaviour persists at small values of $\tau$ (the tests were performed at $\tau=0.0083$), indicating a fundamental flaw of the first-order hydrodynamics equations.

It is worth noting that the generalisation of the conclusions presented in this paper to higher-order extensions of the Chapman-Enskog procedure or of the moment methods is not straightforward, since Ref. cercignani02 warns that the higher orders in the Chapman-Enskog expansion can introduce spurious steady-state solutions, and in certain circumstances, the Chapman-Enskog series may exhibit a divergent behaviour denicol16 .

The present work can be naturally extended to the analysis of dissipative phenomena in fluids composed of massive particles gabbana17a ; gabbana17b or which obey quantum statistics coelho17 . Another extension can be made towards the analysis of dissipation in flows on curved spaces, such as the Bjorken flow in the Milne universe, where the solution of the Boltzmann equation is known semi-analytically florkowski13a ; florkowski13b , as well as the homogeneous and isotropically expanding flow on a background Friedmann-Lemaître-Robertson-Walker (FLRW) space, which was studied analytically and numerically in Refs. bazow16 and tindall17 , respectively.

Acknowledgements

The author is grateful to Prof. Amaresh Jaiswal for useful discussions, as well as to an anonymous Phys. Rev. C referee for carefully reading the manuscript and for many useful comments. This work was supported by a grant of the Romanian National Authority for Scientific Research and Innovation, CNCS-UEFISCDI, project number PN-II-RU-TE-2014-4-2910.

Appendix A Non-dimensionalisation convention

The dimensional form of the Boltzmann equation (1) is:

[TABLE]

where the convention that dimensional quantities are written using a hat was employed. The above equation can be non-dimensionalised using fundamental reference quantities, which are chosen as:

[TABLE]

where $\hat{L}$ is the wavelength of the longitudinal wave, $\hat{c}$ is the speed of light in vacuum and $\hat{n}_{0}$ and $\hat{P}_{0}$ represent the average density and pressure of the medium. From the above fundamental reference quantities, the reference time $\hat{t}_{\rm ref}$, reference momentum $\hat{p}_{\rm ref}$ and reference particle distribution function $\hat{f}_{\rm ref}$ can be derived, as summarised in Tab. 1. The non-dimensional form (1) of the Boltzmann equation can be obtained by multiplying (126) by $\hat{t}_{\rm ref}/\hat{f}_{\rm ref}$, while the non-dimensional relaxation time $\tau$ is given by

[TABLE]

Appendix B Numerical method

In order to solve the AWB equation (1), we employ the relativistic spherical lattice Boltzmann (R-SLB) models introduced in Ref. blaga17prc as an extension of the non-relativistic spherical lattice Boltzmann (SLB) models introduced in Ref. ambrus12 .

A number of $N=100$ nodes are chosen along the $z$ axis, where periodic boundary conditions apply, while the flow is assumed to be homogeneous along the $x$ and $y$ directions. The advection and time evolution are performed using the fifth-order weighted essentially non-oscillatory (WENO-5) rezzolla13 and third-order TVD Runge-Kutta (RK-3) trangenstein07 schemes, as presented in Ref. blaga17prc . The lattice spacing is $\delta z=10^{-2}$. The time step was set to $\delta t=10^{-3}$ for the analysis performed within the frame of the first-order hydrodynamics theory in Sec. III. In the case of the second-order hydrodynamics theory and moment-based method considered in Secs. IV and V, a time step $\delta t=10^{-4}$ was employed to allow an increased temporal resolution (i.e. more data points) for the study of the early time evolution of the longitudinal wave.

The momentum space is factorised using spherical coordinates $p$, $\theta$ and $\varphi$, which are discretised using $Q_{L}$, $Q_{\xi}$ and $Q_{\varphi}$ quadrature points, respectively. The quadrature order along the $p$ direction is set to $Q_{L}=2$, while the azimuthal quadrature order is $Q_{\varphi}=1$. The model thus employs $Q_{L}\times Q_{\xi}\times Q_{\varphi}=2Q_{\xi}$ velocities. The value of $Q_{\xi}$ is chosen depending on the value of the relaxation time $\tau$. As discussed in Ref. blaga17prc , $Q_{\xi}=6$ is sufficient to obtain accurate results at $\tau<0.01$. For values of $\tau$ between $0.01$ and $0.1$, $Q_{\xi}=20$ was employed, while for $\tau\geq 0.1$, $Q_{\xi}$ was set to $Q_{\xi}=200$.

The system is initialised with an equilibrium distribution $f^{(\mathrm{eq})}$ (2) at each point $z_{\ell}=-0.5+(\ell-\frac{1}{2})\delta z$, truncated to $N_{L}=1$ and $N_{\Omega}=5$ with respect to $p$ and $\xi=\cos\theta$, as explained in Ref. blaga17prc .

At a later time $t_{s}=s\,\delta t$ ($s=1,2,\dots T$), the quantities with tilde defined in Eqs. (20) are obtained as:

[TABLE]

In the case of the analysis using the first-order theory performed in Sec. III, $T=20,000$ and the resulting values $\widetilde{\beta}_{s}$, etc. are stored at intervals of $10\delta t=0.01$ resulting in a number of $2,000$ values which are further processed using MathematicaTM to obtain a nonlinear fit, based on the analytic solution for the flow, as described in Secs. III.3 and III.4. Since in the first order theory, $q$ and $\Pi$ cannot be imposed at initial time, the analytic expressions for the evolutions of $\widetilde{q}$ and $\widetilde{\Pi}$ are not accurate at small values of $t$. Thus, the first $50$ points (i.e. up to $t=0.5$) in the data sets were always ignored when performing the nonlinear fits for these quantities.

In the case of the analysis using the second-order theory and the moment method, performed in Secs. IV and V, respectively, the values of the field amplitudes $\widetilde{\beta}_{s}$, etc. were stored at intervals of $\delta t=10^{-4}$. The value of $T$ was obtained using the following algorithm. For $0.001\leq\tau\leq 0.1$, a number of time steps equal to $T=100(\tau/\delta t)$ was considered (i.e. up to $t=100\tau$). This ensured a balanced coverage of the initial stage corresponding to the relaxation of the nonequilibrium parameters $\widetilde{q}$ and $\widetilde{\Pi}$ from their initial vanishing values towards the values predicted by the first order theory (up to $t\sim 5\tau$), as well as of the later stage of the wave evolution, where the attenuation effects dominate. For $0.1<\tau<0.32$, a number of $T=1/(\tau\delta t)$ points was chosen (i.e. corresponding to $t\simeq\tau^{-1}$), while for all $\tau\geq 0.32$, $T=32,000$ time steps were performed (i.e. up to $t=3.2$).