# What is the correct way to calculate deviatoric stress tensor in lattice Boltzmann method?

Due to my previous question, where I asked about flux calculation in lattice Boltzmann (LB) method here, I have more or less same question for deviatoric stress tensor calculation due to pseudo-compressibility of LB method. In fact, strain rate tensor is calculated in LB by using non-equilibrium part of distribution functions ($$f_{i}^{neq}$$) as:

$$\hat{\varepsilon}_{\alpha\beta} = -\frac{1}{2\hat{\tau} \hat{\rho} \hat{c}_{s}^{2}} \sum_{i} f_{i}^{neq} c_{i\alpha}c_{i\beta}$$

Where hat quantities are dimensionless quantities, $$\hat{\tau}$$ is dimensionless relaxation time ($$\hat{\tau} > 0.5$$), $$\hat{\rho}$$ is dimensionless instantaneous density of the fluid, $$\hat{c}_{s}^{2} = \frac{1}{3}$$, and $$c_{i\alpha}$$ is the $$i$$th discrete velocity in $$\alpha$$ direction. Deviatoric stress is defined as:

$$\hat{\sigma}_{\alpha\beta} = 2 \hat{\mu} \hat{\varepsilon}_{\alpha\beta}$$

But, we have for dimensionless viscosity: $$\hat{\mu} = \hat{c}_{s}^{2} (1-\frac{1}{2\hat{\tau}}) \hat{\tau}$$

So, finally:

$$\hat{\sigma}_{\alpha\beta} = -(1-\frac{1}{2\hat{\tau}}) \frac{1}{\hat{\rho}} \sum_{i} f_{i}^{neq} c_{i\alpha}c_{i\beta}$$

This is in contrast with what usually LB people use, such as this one from Krüger et. al. as:

$$\hat{\sigma}_{\alpha\beta} = -(1-\frac{1}{2\hat{\tau}}) \sum_{i} f_{i}^{neq} c_{i\alpha} c_{i\beta}$$

I understand at incompressible limit ($$Mach \rightarrow 0$$), $$\hat{\rho}$$ should be close to 1, but for my simulations where $$Mach \sim 0.06$$ and $$Re \sim 600$$, $$\hat{\rho}$$ may fluctuate quite a bit. So my question is which of these formulas should be used to calculate deviatoric stress? Should I assume $$\hat{\rho} \sim 1$$ even with my high $$Mach$$ number? Any suggestion is truly appreciated.