First, the Thomas-Fermi screening is a semiclassical static theory which assumes that the total potential $\phi(\mathbf{r})$ varies slowly in the scale of the Fermi length $l_{\text{F}}$, the chemical potential $\mu$ is constant and that $T$ is low. In principle, it does not rely on linear response theory.
The condition of slowly varying potential is a general condition of validity of semiclassical models. Physically, if the particle [electron] is represented by a wave packet, what is tellying us is that all the waves in the wavepacket will see the same potential and the particle will suffer [or enjoy!] a force as if it was point-like ["classical"] because such potentials gives rise to ordinary forces in the equation of motion describing the evolution of the position and wavevector of the packet. The wavepacket must have a well-defined wavevector on scale of the Brillouin zone [thus $\Delta k \simeq k_{\text{F}}$] and therefore can be spread in the real space over many primitive cells.
Mathematically, the assumption that your potential is a slowly varying function of the position implies that the theory is not valid for $|\mathbf{q}| \gg k_{\text{F}}$ [and therefore for $|\mathbf{r}| \ll l_{\text{F}}$].
On the other hand, the static Lindhard dielectric function is a fully quantum treatment of the problem and it is valid for all the ranges of $\mathbf{q}$. It includes, in the limit $\mathbf{q} \rightarrow 0$, the linearized Thomas-Fermi dielectric function. It only assumes linear response, that is, the induced density of charge is proportional to the total potential $\phi(\mathbf{r})$.
Note also that the Lindhard treatment is far more general than the Thomas-Fermi in the sense that it can describe both dynamic and static screening.
These Feynman diagrams can be summed by solving the Dyson-Schwinger equation
$$
G = G_0 + G_0\Sigma G
$$
This is a self-consistency equation for $G$. Now write $G_0$ and $G$ in terms of single particle wave functions,
$$ G(x,x';\omega)=\sum_j \phi_j(x)\phi^*_j(x')\left[ \frac{\Theta(E_j-E_F)}{\omega-E_j+i\epsilon} +\frac{\Theta(E_F-E_j)}{\omega-E_j-i\epsilon} \right].
$$
Then the Dyson-Schwinger equation becomes a coupled set of equations for the eigenfunctions $\phi_j$ and the eigenvalues $E_j$. These are the standard Hartee-Fock equations. This is explained in some detail in many text books,
for example Negele and Orland, or Fetter and Walecka.
Best Answer
1) The Thomas-Fermi method is a useful way to estimate the "screening cloud" surrounding electrons in a metal. A quantification of the screening is the inverse dielectric function of the material:
\begin{equation} \frac{\phi_{ext}(\textbf{q})}{\epsilon(\textbf{q})} = \phi_{total}(\textbf{q}) \end{equation}
Slowly-varying in this context means that the wavelength of the perturbation considered is long, i.e. it does not vary much over the lattice spacing or over the Fermi wavelength ($\textbf{q}<<\textbf{k}_F$).
As the wavelength of the perturbation gets shorter and shorter, the perturbation would start to "see" more of the graininess in the electron gas/liquid. Therefore, with a shorter wavelength, it would be difficult to define a "screening cloud" as the perturbation would be varying inside the "cloud".
2) As for the difference between the Thomas-Fermi method and the Hartree-Fock method, the simplest answer is that they seek to calculate different quantities. The Hartree-Fock method is concerned with calculating the many-body wavefunction and the ground state energies by including the effects of electron-electron interactions. This can lead to predictions of experimental quantities such as the low-temperature electronic specific heat as well as the spin susceptibility.
On the other hand, the Thomas-Fermi method seeks to calculate the dielectric function (the screening parameter) -- a quantity that can be probed by e.g. inelastic electron scattering.