Thermal LBM convection instability

Hi, I am trying to understand how thermal LBM works, for that purpose I wrote simple code in Python from formulas I found in nice review paper by Sharma et al. My code use double distribution for momentum and temperature calculation.

It looks like that it is working fine for momentum only or temperature only, but not for both, e. g. I quickly got instability when try to simulate natural convection in 2D.

r/CFD - Thermal LBM convection instability
Evolution of instability of natural convection. Top row - mass, bottom row - energy (temperature).

Maybe somebody got any insights? where could problem be?

Code here