Grid resolution dependence in the rayleighTaylor2D example

We have modified the rayleighTaylor2D exapmle in order to be able to parametrize the case by defining Reynolds and Atwood numbers, since they are the two numbers driving the phenomenon.

We have followed carefully the physical-to-lattice units conversion instructions provided by Timm Krüger here and we are running 4 cases at Re=256,2048 for ny=512,1024 and nx=ny/4 in order to observe one mode.

//Aspect ratio of the domain is 4
nx = ny/4;

//Transfer from physical to lattice quantities
T dx = W/nx;                  		 //W = 1
T uMaxLB=0.2;                               //low lattice maximum velocity to keep low lattice Mach number (LBM regime) 
T uMax=sqrt(W*g);                        //He,Chen,Zhang99, v/sqrt(Wg) = 0.2
T dt = dx*uMaxLB/uMax;
T nu_LB = nx*uMaxLB/Reynolds;                       	//I don't know which approach is the best one, both are buggy anyway
T force_LB = pow(Reynolds*nu_LB/nx/sqrt(nx),2);
nu = nu_LB*dx*dx/dt;
omega = 2./(6.*nu_LB+1.);
g2 = pow(Reynolds*nu/W/sqrt(W),2);

and then instead of applying -force_LB we apply -2*A/(A+1)*force to model the influence of the Atwood number. (but that's not important for this question).

However our simulations seem to be very diffusive (more for Re=256 and less for Re=2048, but still). The diffusion is reduced when increasing the resolution of our mesh from ny 512 to 1024. This seems to be in contradiction with a paper we found online Do you have any idea of why this is happening or how to avoid it, so that we can reproduce the known contours of the RTI?