Hi,
It seems that, as you are saying, the paper fails to give a hint on how to compute grad(phi) in Eq. (8). Therefore, you don’t know how to translate the inter-particle potential phi into a force term to be used in the LB scheme. Although I know nothing more than you about the original intent of the authors of this paper, I do have a suggestion which is based on other multi-phase models.
The following is a general formula for the numerical approximation to the gradient of a scalar field phi, with an accuracy that is consistent with the overall accuracy of LB schemes (second-order in space):
grad(\phi) = 1/c_s^2 \sum_i t_i c_i \phi(x+c_i),
where c_s^2 is the lattice constant corresponding to the square-speed-of-sound (usually 1/3), t_i are the lattice weights (4/9, 1/9 and 1/36 for D2Q9) and c_i the lattice velocities ( (-1,1) is for example one of the nine lattice velocities in D2Q9), and x is the point is space where grad(phi) is evaluated. All units are lattice units (dx=dt=1). It is easy to show that the above formula is right, by Taylor-expanding phi around the position x:
\phi(x+c_i) = \phi + c_{i\beta}\partial_\beta \phi + 1/2 \partial_{i\beta}\partial_{i\gamma} c_{i\beta} c_{i\gamma} \phi .
Plug this approximation into the above formula for grad(\phi), and use the lattice symmetries (section 2.1.1 on page 15 in my thesis). You will find that many terms cancel, and the only one left is
\sum_i t_i c_{i\alpha} \phi(x+c_i) = c_s^2\partial_\alpha \phi .
Of course, this is just one possible choice for the numerical evaluation of a gradient. Although I have no quantitative argument to support it against, say, a simple finite difference scheme along each of the main axes, I would say that we know from numerical experiment that it is a good choice for multi-phase fluids. Among others, it is the approach chosen to evaluate grad(phi) in the widely used Shan/Chen multi-phase model. As another qualitative argument, one can argue that this approach is “more isotropic” than straight finite differences, because it also takes into account diagonal directions. It is therefore likely to do a better job for modeling the shape of an interface between phases.