This is explained nicely in arXiv:1012.0461 [cond-mat.stat-mech] (section 2.1) in the case where there are no thermal effects. I leave out the thermal effects in the following. Since you know your paper better than me you should be able to generalise this.
The projection comes from the incompressibility of the fluid. The Navier-Stokes equations read
∂tvα+vβ∂βvα=ν∂β∂βvα−∂αP+fα,∂βvβ=0.
ν is the kinematic viscosity, ∂t=∂/∂t, ∂β=∂/∂xβ and upper indices indicate the components of the corresponding vector. Repeated indices are summed over.
The pressure can be eliminated from these equations by taking the divergence of Navier-Stokes equations,
∂α(∂tvα+vβ∂βvα=ν∂β∂βvα−∂αP+fα).
Taking incompressibility into account, this leads to
∂α∂αP=−∂αvβ∂βvα+∂αfα, or P=−1∇2(∂γvβ∂βvγ−∂γfγ).
Then Navier-Stokes equations are written as
∂tvα+(δα,γ−∂α1∇2∂γ)vβ∂βvγ=ν∂β∂βvα+(δα,γ−∂α1∇2∂γ)fγ,
after inserting the second expression for the pressure. Finally, since the velocity field is incompressible,
(δα,γ−∂α1∇2∂γ)vγ=vα,
we can define new fields as
ˆvα=(δα,γ−∂α1∇2∂γ)vγ, and ˆfα=(δα,γ−∂α1∇2∂γ)fγ,
in terms of which Navier-Stokes equations become
∂tˆvα+ˆvβ∂βˆvα=ν∂β∂βˆvα+ˆfα,∂βvβ=0.
The pressure term has been completely eliminated. The price to pay is that the velocity field is implicitly constrained to be incompressible and that the force is replaced by ˆf.
If you compute the correlator of ˆf in Fourier space you get,
⟨ˆfα(p)ˆfβ(q)⟩=∫x,yeipx+qy(δα,γ1−∂α1∇2∂γ1)(δβ,γ2−∂β1∇2∂γ2)⟨fγ1(x)fγ2(y)⟩
=(δα,γ1−pαpγ1p2)(δβ,γ2−qβqγ2q2)∫x,yeipx+qy⟨fγ1(x)fγ2(y)⟩
Finally, inserting
∫x,yeipx+qy⟨fγ1(x)fγ2(y)⟩=δγ1,γ2δ(p+q)D0(p),
and contracting the indices leads to your expression for the correlator of ˆf.