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.