Relativistic bulk viscous fluids of Burgers type and their presence in neutron stars

We consider a relativistic fluid mixture [24] with a single conserved particle current n^µ (e.g. the baryon current) and an arbitrary number of non-conserved currents $n_a^\mu$, where a is a chemical index. We assume that all the fluid tensors are isotropic in the rest frame, defined by the four-velocity $u^\mu \propto n^\mu$, so that the constitutive relations take the form [24–27]

$$\begin{equation} \begin{split} T^{\mu \nu} = {}& (\rho+P)u^\mu u^\nu + P g^{\mu \nu} \, , \\ s^\mu = {}& \mathfrak{s} n u^\mu \, , \\ n^\mu = {}& n u^\mu \, , \\ n_a^\mu = {}& Y_a n u^\mu \, , \\ \end{split} \end{equation} \tag{ 5 }$$

where $T^{\mu \nu}$ is the stress-energy tensor and s^µ is the entropy four-current. The scalar fields ρ, P, $\mathfrak{s}$ and Y_a are interpreted respectively as the energy density, the pressure, the specific entropy and the a-th non-conserved fraction (all measured in the rest frame). The thermodynamics of the fluid is most conveniently described in the 'per-particle representation'. In particular, defined the specific volume $v = n^{-1}$ and the specific energy $u = \rho/n$, we write an equation of state of the form $u = u(\mathfrak{s},v,Y_a)$, whose differential reads [28] (we adopt Einstein's convention for repeated chemical indices)

$$\begin{equation} du = Td \mathfrak{s}-Pdv -\mathbb{A}^a dY_a \, , \end{equation} \tag{ 6 }$$

where T is the temperature and $-\mathbb{A}^a$ is the chemical potential of the species a. Since the fractions Y_a are not conserved, in chemical equilibrium we must have $\mathbb{A}^a = 0$ (due to the maximum entropy principle [28–30]), so that we can interpret $\mathbb{A}^a$ as reaction affinities [31, 32]. The equations of motion of the system are the conservation laws $\nabla_\mu T^{\mu \nu} = 0$ and $\nabla_\mu n^\mu = 0$, and the particle production equations: $\nabla_\mu n^\mu_a = \mathcal{R}_a$. If the fluid is not too far from local equilibrium, we can expand the reaction rates $\mathcal{R}_a$ to first order in the affinities, $\mathcal{R}_a = \Xi_{ab}\mathbb{A}^b$ [33, 34], so that we have the field equations

$$\begin{equation} u^\mu \nabla_\mu Y_a = v \Xi_{ab}\mathbb{A}^b \, . \end{equation} \tag{ 7 }$$

Given that particle fractions are invariant under time reversal, the Onsager-Casimir principle [35, 36] requires that the reaction matrix $\Xi_{ab}$ be symmetric. Furthermore, it must also be non-negative definite, since the second law of thermodynamics ($\nabla_\mu s^\mu \geqslant 0$ [37]), combined with all other field equations, implies

$$\begin{equation} T u^\mu \nabla_\mu \mathfrak{s} = v\Xi_{ab}\mathbb{A}^a\mathbb{A}^b \geqslant 0 \, . \end{equation} \tag{ 8 }$$

Since the fractions Y_a are not conserved, we can actually assume that $\Xi_{ab}$ is positive definite, and therefore invertible. We call its matrix inverse $\Xi^{ab}$, so that $\Xi^{ab}\Xi_{bc} = \delta{^a _{{\,\,}c}}$.

Now that we have stated the equations of motion of relativistic fluid mixtures, we can proceed to prove that such substances are indeed bulk viscous fluids. In particular, we can show that, if the dynamics is sufficiently slow, solutions of the fluid equations asymptotically relax [5] to the constitutive relations of the relativistic Navier–Stokes theory [38] with only bulk viscosity.

For the analysis that follows, it is particularly convenient to treat the collection of fields $\varphi_i = \{u^\mu, \mathfrak{s},v,\mathbb{A}^a \}$ as our independent degrees of freedom. This means that we will regard all the physical tensors in (5) as functions of these fields. In particular, all thermodynamic quantities f from now on are understood as functions $f(\mathfrak{s},v,\mathbb{A}^a)$, and partial derivatives are performed accordingly. For example, if we write $\partial f/\partial v$, it is understood that the variables which are held constant are $\mathfrak{s}$ and $\mathbb{A}^a$. If we write $\partial f/\partial \mathbb{A}^a$, it is understood that we are holding constant v, $\mathfrak{s}$ and all other $\mathbb{A}^b$ (with b ≠ a). Then, if the fluid is close to local thermodynamic equilibrium (i.e. if $\mathbb{A}^a$ are small), we can make the following first-order expansions:

$$\begin{equation} \begin{split} \rho(\mathfrak{s},v,\mathbb{A}^a) \approx{}& \rho(\mathfrak{s},v,0) +\dfrac{\partial \rho}{\partial \mathbb{A}^a} \mathbb{A}^a \, , \\ P(\mathfrak{s},v,\mathbb{A}^a) \approx{}& P(\mathfrak{s},v,0) +\dfrac{\partial P}{\partial \mathbb{A}^a} \mathbb{A}^a \, . \\ \end{split} \end{equation} \tag{ 9 }$$

Invoking equation (6), we immediately see that the partial derivative $\partial \rho/\partial \mathbb{A}^a$, being evaluated at $\mathbb{A}^b = 0$, vanishes (recall that $\rho = u/v$). This is a manifestation of the minimum energy principle [28]. The partial derivative $\partial P/\partial \mathbb{A}^a$ can be rewritten in a more illuminating form. In fact, defined the thermodynamic potential $\mathcal{G} = u+\mathbb{A}^a Y_a$, we have the differential $d\mathcal{G} = Td \mathfrak{s}-Pdv+Y_a d\mathbb{A}^a$, which can be used to derive the following Maxwell relation:

$$\begin{equation} \dfrac{\partial P}{\partial \mathbb{A}^a} = -\dfrac{\partial Y_a}{\partial v} \, . \end{equation} \tag{ 10 }$$

Physically, this equation is telling us that the susceptibility of the pressure to chemical imbalances equals the susceptibility of the chemical fractions to a volume expansion. Combining these results together, we find that the stress energy tensor in (5), expanded to first order in $\mathbb{A}^a$, takes the (Eckart-frame [38–40]) bulk viscous form

$$\begin{equation} T^{\mu \nu} = (\rho_{\text{eq}}+P_{\text{eq}}+\Pi)u^\mu u^\nu +(P_{\text{eq}}+\Pi)g^{\mu \nu} \, , \end{equation} \tag{ 11 }$$

where we are adopting the notation $f_{\text{eq}}(\mathfrak{s},v) = f(\mathfrak{s},v,0)$ and we have introduced the bulk viscous stress

$$\begin{equation} \Pi = -\dfrac{\partial Y_a}{\partial v} \mathbb{A}^a \, . \end{equation} \tag{ 12 }$$

Let us now verify explicitly that the mixture indeed admits a Navier–Stokes regime where Π is given by (2). In order to do this, first we use the chain rule to rewrite equation (7) as follows:

$$\begin{equation} \dfrac{\partial Y_a}{\partial v} u^\mu \nabla_\mu v + \dfrac{\partial Y_a}{\partial \mathfrak{s}} u^\mu \nabla_\mu \mathfrak{s} + \dfrac{\partial Y_a}{\partial \mathbb{A}^b} u^\mu \nabla_\mu \mathbb{A}^b = v\Xi_{ab}\mathbb{A}^b \, . \end{equation} \tag{ 13 }$$

If we retain only the first order terms in $\mathbb{A}^a$ we have that the contribution proportional to $u^\mu \nabla_\mu \mathfrak{s}$ can be neglected, see equation (8). Furthermore, we can use the equation $\nabla_\mu n^\mu = 0$ to prove that $u^\mu \nabla_\mu v = v\nabla_\mu u^\mu$ [25], so that, contracting both sides of (13) with $n\Xi^{ab}$ (which is the matrix inverse of $v\Xi_{ab}$), we obtain

$$\begin{equation} \begin{split} & \tau{^a_{{\,\,}b}} u^\mu \nabla_\mu \mathbb{A}^b + \mathbb{A}^a = \kappa^a \nabla_\mu u^\mu \, ,\\ & \text{with} \quad \tau{^a_{{\,\,}b}} = -n\Xi^{ac}\dfrac{\partial Y_c}{\partial \mathbb{A}^b} \, , \quad \kappa^a = \Xi^{ab}\dfrac{\partial Y_b}{\partial v} \, . \\ \end{split} \end{equation} \tag{ 14 }$$

By the relaxation effect [3], we know that in the limit of a very slow process the quantities $u^\mu \nabla_\mu \mathbb{A}^b$ are negligible compared to $\mathbb{A}^a$, so that $\mathbb{A}^a \approx k^a\nabla_\mu u^\mu$. It follows that equation (12) can be approximated as $\Pi \approx -\zeta \nabla_\mu u^\mu$, with

$$\begin{equation} \zeta = \Xi^{ab} \dfrac{\partial Y_a}{\partial v}\dfrac{\partial Y_b}{\partial v} \, . \end{equation} \tag{ 15 }$$

Clearly, ζ is non-negative, because $\Xi_{ab}$ (and thus also $\Xi^{ab}$) is positive definite. This completes our proof that a 'slow' fluid mixture obeys the Navier–Stokes constitutive relations for bulk viscosity.

In what follows, we will restrict our attention to flows that are 'almost incompressible'. This means that, fixed a reference incompressible flow (i.e. fixed an arbitrary solution of the fluid equations with $\nabla_\mu u^\mu = 0$), we will study neighbouring compressible solutions to first order in perturbation theory around such reference incompressible flow. This is needed because the Burgers equation is a linear rheological model, and it holds only for small $\nabla_\mu u^\mu$ and Π.

Let us set up the perturbative expansion rigorously. We consider a smooth one-parameter family of solutions $\{\varphi_i (\epsilon),g_{\mu \nu}(\epsilon)\}$ of the fluid equations coupled with gravity, where $\{\varphi_i(0),g_{\mu \nu}(0)\}$ is an incompressible flow, i.e. $\nabla_\mu u^\mu(0) = 0$ across all spacetime. Such incompressible flow $\{\varphi_i(0),g_{\mu \nu}(0)\}$ may be both fast rotating and accelerating, and it may admit strong shear flows (and large gradients in general), but it does not expand. Note that we are not assuming that the fluid itself is incompressible: We are just considering a particular incompressible solution of the fluid equations (e.g. a star in hydrostatic equilibrium). Furthermore, we also assume that $\mathbb{A}^a(0)$ vanishes on some initial Cauchy surface. Then, equation (14) implies that $\mathbb{A}^a(0)$ vanishes everywhere, meaning that the solution ε = 0 is reversible: $[u^\mu\nabla_\mu \mathfrak{s}](0) = 0$, and also $\Pi(0) = 0$, see equation (12). Thus, since by chain rule

$$\begin{equation} u^\mu \nabla_\mu f = \dfrac{\partial f}{\partial v} v \nabla_\mu u^\mu + \dfrac{\partial f}{\partial \mathfrak{s}} u^\mu \nabla_\mu \mathfrak{s} + \dfrac{\partial f}{\partial \mathbb{A}^b} u^\mu \nabla_\mu \mathbb{A}^b \, , \end{equation} \tag{ 16 }$$

we see that all thermodynamic quantities $f(\mathfrak{s},v,\mathbb{A}^b)$ are conserved along the flow worldlines: $[u^\mu \nabla_\mu f](0) = 0$.

Now, we linearise equation (14) to first order in ε, i.e. we differentiate (14) in ε and evaluate the result at ε = 0 [6]. This corresponds to studying (14) in a regime of small compression ² . Introducing the compact notation $Q: = Q(0)$ and $\delta Q: = dQ(0)/d\epsilon$, for any field Q, we have the following linear dynamics:

$$\begin{equation} \tau{^a_{{\,\,}b}} u^\mu \nabla_\mu \delta \mathbb{A}^b + \delta \mathbb{A}^a = \kappa^a \, \delta(\nabla_\mu u^\mu) \, . \end{equation} \tag{ 17 }$$

Here, $\tau{^a _{{\,\,}b}}$ and κ^a play the role of background quantities, as they are evaluated on the reference incompressible flow (so that $u^\mu \nabla_\mu \tau{^a _{{\,\,}b}} = u^\mu \nabla_\mu \kappa^a = 0$), while $\delta \mathbb{A}^a$ and $\delta(\nabla_\mu u^\mu)$ are first-order perturbation fields. In the case of only two non-equilibrium chemical fractions, equation (17) can be expanded into the system

$$\begin{equation} \begin{split} \tau{^1 _{{\,\,}1}} u^\mu \nabla_\mu \delta \mathbb{A}^1+\tau{^1 _{{\,\,}2}} u^\mu \nabla_\mu \delta \mathbb{A}^2+ \delta \mathbb{A}^1 = {}& \kappa^1 \, \delta(\nabla_\mu u^\mu) \, , \\ \tau{^2 _{{\,\,}1}} u^\mu \nabla_\mu \delta \mathbb{A}^1+\tau{^2 _{{\,\,}2}} u^\mu \nabla_\mu \delta \mathbb{A}^2+ \delta \mathbb{A}^2 = {}& \kappa^2 \, \delta(\nabla_\mu u^\mu) \, .\\ \end{split} \end{equation} \tag{ 18 }$$

With some simple algebra, we can rewrite the equations above as follows:

$$\begin{align} \begin{split} (\det \tau) \, u^\nu \nabla_\nu (u^\mu \nabla_\mu \delta\mathbb{A}^1) + (\text{Tr} \, \tau) \, u^\mu \nabla_\mu \delta \mathbb{A}^1 + \delta \mathbb{A}^1 = {}& \kappa^1 \delta(\nabla_\mu u^\mu) +(\tau{^2 _{{\,\,}2}}\kappa^1-\tau{^1 _{{\,\,}2}}\kappa^2)\, u^\nu \nabla_\nu \, \delta(\nabla_\mu u^\mu)\, , \\ (\det \tau) \, u^\nu \nabla_\nu (u^\mu \nabla_\mu \delta\mathbb{A}^2) + (\text{Tr} \, \tau) \, u^\mu \nabla_\mu \delta \mathbb{A}^2 + \delta \mathbb{A}^2 = {}& \kappa^2 \delta(\nabla_\mu u^\mu) +(\tau{^1 _{{\,\,}1} }\kappa^2-\tau{^2 _{{\,\,}1}}\kappa^1)\, u^\nu \nabla_\nu \, \delta(\nabla_\mu u^\mu)\, , \\ \end{split} \end{align} \tag{ 19 }$$

where $\det \tau$ and $\text{Tr} \, \tau$ are respectively the determinant and the trace of the matrix $\tau = [\tau{^a _{{\,\,}b}}]$. Recalling equation (12), we can finally combine the two field equations for $\delta \mathbb{A}_1$ and $\delta \mathbb{A}_2$ to have an equation of motion for $\delta \Pi$ of the form

$$\begin{align} \lambda_2 \, u^\nu \nabla_\nu (u^\mu \nabla_\mu \delta\Pi) + \lambda_1 \, u^\mu \nabla_\mu \delta \Pi+ \delta \Pi = -\zeta \, \delta(\nabla_\mu u^\mu) -\chi\, u^\nu \nabla_\nu \, \delta(\nabla_\mu u^\mu)\, , \end{align} \tag{ 20 }$$

where ζ is given by equation (15), and

$$\begin{equation} \begin{split} \lambda_2 = {}& \det \tau \, ,\\ \lambda_1 = {}& \text{Tr} \, \tau \, ,\\ \chi = {}& (\det \tau) \dfrac{\partial Y_a}{\partial v} (\tau^{-1}){^a _{{\,\,}b}} \kappa^b \, . \\ \end{split} \end{equation} \tag{ 21 }$$

Equation (20) is the central formula of this manuscript. It tells us that, for a mixture with two non-conserved fractions, the bulk viscous stress Π obeys an equation of motion that is of second order in time. This reflects the fact that there are two algebraic non-equilibrium degrees of freedom ($\mathbb{A}^1$ and $\mathbb{A}^2$). Indeed, if one wants to solve equation (20), they need to prescribe not only the initial value of Π, but also its initial time derivative, namely $u^\mu \nabla_\mu \Pi$.

We would like to stress that equation (20) is not an expansion in powers of '$\, u^\mu \nabla_\mu \,$'. In fact, in the derivation, no assumption was made about how fast the process is. The only approximation that we made was the linearization in $\delta \mathbb{A}^a$ and $\delta (\nabla_\mu u^\mu)$, which is an assumption on the amplitude of the perturbation, and not on its frequency. Indeed, equation (20) well approximates (7) at all frequencies. Despite this, it is clear that the Burgers model remains applicable also in dynamical regimes with small frequency and large amplitude, as it reduces to Navier–Stokes by the relaxation effect [3], with the correct bulk viscosity coefficient ζ.

Now we only need to show that equation (20) is indeed the Burgers equation for viscoelastic matter. In order to do so, we must prove that its linear-response Green function is (4). The proof goes as follows. The matrices $-n\partial Y /\partial \mathbb{A} = [-n\partial Y_a/\partial \mathbb{A}^b]$ and $\Xi = [\Xi_{ab}]$ are both symmetric and positive definite. Therefore, there exist an invertible real matrix $\mathcal{N} = [\mathcal{N}{_{(a)} ^{{\,\,}b}}]$ and a diagonal positive definite matrix $\Lambda = \text{diag}\{\tau_{(1)},\tau_{(2)}\}$ such that [42]

$$\begin{equation} {-}n \dfrac{\partial Y}{\partial \mathbb{A}} = \mathcal{N}^{-1}\mathcal{N}^{-T} \, , \quad \quad \quad \Xi = \mathcal{N}^{-1} \Lambda^{-1} \,\mathcal{N}^{-T} \, . \end{equation} \tag{ 22 }$$

Then, the matrix $\Xi^{-1} = [\Xi^{ab}]$ decomposes into $\Xi^{-1} = \mathcal{N}^{T} \Lambda\mathcal{N}$, and we find

$$\begin{equation} \tau = \Xi^{-1} \bigg({-}n \dfrac{\partial Y}{\partial \mathbb{A}} \bigg) = \mathcal{N}^T \Lambda \mathcal{N}^{-T} \, . \end{equation} \tag{ 23 }$$

Introducing the notation

$$\begin{equation} \dfrac{\zeta_{(1)}}{\tau_{(1)}} = \bigg( \mathcal{N}{_{(1)}^{{\,\,}b}} \dfrac{\partial Y_b}{\partial v} \bigg)^2 \, , \quad \quad \quad \dfrac{\zeta_{(2)}}{\tau_{(2)}} = \bigg( \mathcal{N}{_{(2)}^{{\,\,}b}} \dfrac{\partial Y_b}{\partial v} \bigg)^2 \, , \end{equation} \tag{ 24 }$$

Equations (15) and (21) can be rewritten as follows:

$$\begin{equation} \begin{split} \lambda_2 = {}& \tau_{(1)}\tau_{(2)} \, , \\ \lambda_1 = {}& \tau_{(1)}+\tau_{(2)} \, , \\ \zeta = {}& \zeta_{(1)}+\zeta_{(2)} \, , \\ \chi = {}& \zeta_{(1)}\tau_{(2)}+\zeta_{(2)}\tau_{(1)} \, . \\ \end{split} \end{equation} \tag{ 25 }$$

Thus, if we work in a global coordinate system such that $u^\mu = \delta^\mu_t$ on all events, equation (20) reduces to

$$\begin{align} \tau_{(1)}\tau_{(2)}\partial^2_t \delta \Pi +(\tau_{(1)}+\tau_{(2)})\partial_t \delta \Pi +\delta \Pi = & -(\zeta_{(1)}+\zeta_{(2)}) \, \delta(\nabla_\mu u^\mu)\nonumber\\ & -(\zeta_{(1)}\tau_{(2)}+\zeta_{(2)}\tau_{(1)})\partial_t \, \delta(\nabla_\mu u^\mu) \, . \end{align} \tag{ 26 }$$

Let us now focus on the Green function (4). If we plug this choice of $G(t^{^{\prime}})$ into equation (1), we find that $\delta\Pi$ can be expressed as the sum of two contributions, $\delta\Pi_{(1)}$ and $\delta\Pi_{(2)}$, each of which obeys an independent Israel-Stewart-type equation:

$$\begin{equation} \begin{split} (\tau_{(1)}\partial_t +1)\delta\Pi_{(1)} = {}& -\zeta_{(1)} \, \delta(\nabla_\mu u^\mu) \, ,\\ (\tau_{(2)}\partial_t +1)\delta\Pi_{(2)} = {}& -\zeta_{(2)} \, \delta(\nabla_\mu u^\mu) \, .\\ \end{split} \end{equation} \tag{ 27 }$$

Applying the operator $\tau_{(2)}\partial_t+1$ to the first equation, and the operator $\tau_{(1)}\partial_t+1$ to the second equation, and adding together the resulting formulas, we indeed recover (26). This shows that the dynamics described by (20) arises from the Burgers Green function (4), which is what we wanted to prove.

If we set the right-hand side of (26) to zero (incompressible evolution), we find that the bulk stress obeys the equation $\lambda_2 \partial^2_t \delta \Pi + \lambda_1 \partial_t \delta \Pi+\delta \Pi = 0$, which describes the dynamics of a damped harmonic oscillator [2]. Such oscillator is necessarily overdamped, i.e. $\lambda_1^2-4\lambda_2 \geqslant 0$, which follows from equation (25), and from the fact that $\tau_{(1)}$ and $\tau_{(2)}$ are always positive. The implication is that the bulk stress cannot really 'oscillate'. Instead, the evolution of $\delta \Pi$ is the superposition of two exponential relaxations. Indeed, it is straightforward to show that the chemical mixture has two non-hydrodynamic modes with purely imaginary frequency gap: $\omega_{(a)}(k = 0) = -i/\tau_{(a)}$. This is not a surprise, since it is well known that chemical oscillations are forbidden in the linear regime [43]. This is a consequence of the Onsager symmetry of $\Xi_{ab}$, which forces all non-hydrodynamic gaps to lay on the imaginary axis [7].

Our goal is to prove that equation (20) exactly reproduces the results of Alford, Harutyunyan, and Sedrakian [46, 48] about bulk viscosity in hot and dense $npe\mu$ matter. To this end, we need first to express the Burgers transport coefficients λ₂, λ₁, ζ, and χ in terms of the quantities A_n , A_p , A₁, A₂, C₁, C₂, λ_e , and λ_µ introduced in [46, 48]. Comparing our formalism with that of [48], going through some simple (albeit tedious) algebra, one can show that, in the ordered chemical basis $\{e,\mu\}$, the following 'dictionary relations' hold:

$$\begin{equation} \begin{split} -n\dfrac{\partial Y}{\partial \mathbb{A}} = {}& \begin{bmatrix} A_1 & A_n+A_p \\ A_n+A_p & A_2 \end{bmatrix}^{-1}, \\ -\dfrac{\partial Y}{\partial v} = {}& \begin{bmatrix} A_1 & A_n+A_p \\ A_n+A_p & A_2 \end{bmatrix}^{-1} \begin{bmatrix} C_1 \\ C_2 \\ \end{bmatrix} ,\\ \Xi = {}& \begin{bmatrix} \lambda_e & 0 \\ 0 & \lambda_\mu \\ \end{bmatrix} . \\ \end{split} \end{equation} \tag{ 28 }$$

Plugging these formulas into (15) and (21), and introducing the compact notation $\mathfrak{D} = A_1 A_2-(A_n+A_p)^2$, we obtain

$$\begin{equation} \begin{split} \lambda_2 = {}& (\mathfrak{D} \lambda_e \lambda_\mu )^{-1} \, , \\ \lambda_1 = {}& \dfrac{1}{\mathfrak{D}} \bigg( \dfrac{A_2}{\lambda_e} + \dfrac{A_1}{\lambda_\mu} \bigg) \, , \\ \zeta = {}& \dfrac{[C_1 A_2-C_2(A_n+A_p)]^2}{\mathfrak{D}^2\lambda_e} + \dfrac{[C_2 A_1-C_1(A_n+A_p)]^2}{\mathfrak{D}^2\lambda_\mu} \, , \\ \chi = {}& \dfrac{A_2 C_1^2+A_1 C_2^2-2(A_n+A_p)C_1C_2}{\mathfrak{D}^2 \lambda_e \lambda_\mu} \, . \\ \end{split} \end{equation} \tag{ 29 }$$

It is immediate to verify that our formula for ζ coincides with equation (43) of [48]. This confirms that the Burgers model has the correct 'infrared behaviour' in the Navier–Stokes limit.

In order to compare the dynamics of the Burgers model with that of the multicomponent model of [46, 48], we need to work in the same physical setting, and compare analogous quantities. The analysis of [46, 48] focuses on small periodic oscillations, so that we need to set $\delta Q \propto e^{-i\omega t}$ (note the different sign convention for ω in [46, 48]) for all quantities $\delta Q$, where the frequency $\omega \in \mathbb{R}$ is not necessarily small, in the sense that we may also have $|\omega\tau{^a_{{\,\,}b}}| \gg 1$. Then, equation (20) becomes

$$\begin{equation} \delta \Pi = - \dfrac{\zeta - i\chi \omega}{1-i\lambda_1\omega -\lambda_2 \omega^2 } \delta(\nabla_\mu u^\mu) \, . \end{equation} \tag{ 30 }$$

If we plug (29) into (30), we obtain a complicated formula for $\delta \Pi$, which is different from the corresponding equation (37) of [48]. The reason is that, in section 2.2, we have defined Π as the deviation of the total pressure from the state of chemical equilibrium, while Alford, Harutyunyan, and Sedrakian quantify Π as its deviation from the pressure at frozen fractions (see also [49]). Hence, we are just comparing different physical quantities. However, if our analysis is correct, the effective viscosity coefficient $\zeta_{\text{eff}}(\omega)$, defined by the relation $\zeta_{\text{eff}} = -\mathfrak{Re} [\delta \Pi/ \delta (\nabla_\mu u^\mu)]$, must be the same, since its value can be inferred from the free-energy dissipation rate [13, 17, 23]. Indeed, from equation (30), we get

$$\begin{align} \zeta_{\text{eff}}(\omega) = \dfrac{\lambda_e \lambda_\mu \big\{\lambda_e [(A_n+A_p)C_1-A_1 C_2]^2+\lambda_\mu[(A_n+A_p)C_2-A_2 C_1]^2 \big\}+\omega^2 (\lambda_e C_1^2+\lambda_\mu C_2^2) }{\big\{\lambda_e \lambda_\mu [A_1A_2-(A_n +A_p)^2]-\omega^2 \big\}^2 +\omega^2(\lambda_e A_1+\lambda_\mu A_2)^2} \, , \end{align} \tag{ 31 }$$

which perfectly agrees with equation (40) of [48]. This shows that $npe\mu$ matter is accurately described, at all hydrodynamic frequencies, by the Burgers model. Needless to say that, instead, the Israel-Stewart theory, whose effective viscosity coefficient is a Lorentzian function, $\zeta_{\text{eff}}(\omega) = \zeta/(1+\omega^2 \tau^2)$ [17], cannot reproduce equation (31), unless λ_e or λ_µ vanishes or diverges.

The astrophysical example of $npe\mu$ matter discussed above shows that bulk viscous fluids of Burgers type are not a mere mathematical conjecture. They exist in nature. More importantly, deviations from both Navier–Stokes and Israel-Stewart enter the formula for $\zeta_{\text{eff}}(\omega)$, thereby modifying the dissipation rate of sound waves. Whether these effects have a measurable impact on the damping rate of neutron star oscillations is left as a future direction of investigation.