I’ve read the article in choosing the units in lattice Boltzmann simulations by Jonas Latt but there’s something still bugging me: According to the article, dt should be chosen to fullfill : dt ~ dx² . The lattice speed would then be c = dx/dt = 1/dx .

Now in other works e.g. Master thesis Goetz, the lattice speed is often chosen to be c=1 so it wouldn’t appear in the equilibrium distribution function
feq = warho(1 + 3eau/c² + 9/2*(ea*u)²/c^4 - 3/2 * u²/c²).

My question then would be, if we chose dt ~ dx² so that c is not 1, do we still have to incorporate the c in feq ?

In lieu of someone better qualified to answer your question I’d say yes, you have to. First, I didn’t find the c=1 trick in the masters thesis at first glance. However, such things are often used in physics just because it is possible (You can also define units such that all the fundamental constants - hbar and so on - become unity). However, a result such obtained would still have to be retransformed to be useful. Second, even if you want to use this trick, it doesn’t mean you cannot follow \delta x^2~\delta t. You have to proportionality factors to play with - \delta t/t_p and \delta x/x_p. Hence, if you choose \delta x to get c to unity, this just means you have to scale \delta t accordingly to keep the Mach error at bay. But, for all practical reasons, - screw c=1. Just calculate it and choose your units according to the physics of your system…

As I said, all my humble opinion only, those better informed please feel free to correct me.