Let $f\in\mathbb{Z}[X_1,\ldots,X_n]$ be a Diophantine equation which, for the purposes of this question, I will assume is homogeneous and nonsingular on $\mathbb{R}^n\setminus\{0\}$ (so that $\nabla f\not=0$). Supposing that it has infinitely many primitive integer zeros, we can posit that they are smoothly distributed in an asymptotic sense. Writing $V(R)\subseteq R^n$ for the set of primitive solutions to $f(x)=0$ in a ring $R$, the integer solutions $V(\mathbb{Z})$ clearly lie on the manifold $V(\mathbb{R})$. So, I am looking for a density $\rho\colon V(\mathbb{R})\to\mathbb{R}$ with $$\vert V(\mathbb{Z})\cap U\vert\sim\int_{V(\mathbb{R})\cap U}\rho(x)\,d\sigma(x),\qquad\qquad{\rm(1)}$$ for subsets $U\subseteq\mathbb{R}^n$, where $d\sigma$ is the standard surface integral on $V(\mathbb{R})$. This should hold asymptotically as $U$ is scaled up, and for reasonably regular regions $U$.
My question is regarding a simple (but incorrect -- see below) heuristic argument for calculating $\rho$. Choosing positive integer $N$ and real $a\gg N$ then, for large regions $U$, the set of $x\in U$ with $\vert f(x)\vert < 2a$ has volume about $2a\int_{V(\mathbb{R})\cap U}\Vert\nabla f\Vert^{-1}\,d\sigma$, so should contain about that number of integer points. The probability of a random $x\in\mathbb{Z}^n$ being relatively prime to $N$ and satisfying $f(x)=0$ (mod $N$) is $N^{-n}\vert V(\mathbb{Z}/N\mathbb{Z})\vert$. Conditional on $\vert f(x)\vert < 2a$ and $f(x)=0$ (mod $N$), it seems reasonable to suppose that $f(x)=0$ with probability $N/(2a)$. Multiplying these terms together and taking the limit as $N$ increases to include all prime-powers as factors, we get the following expression for $\rho$. \begin{align} &\rho(x)=\Vert\nabla f(x)\Vert^{-1}\prod_p c_p,\qquad\qquad{\rm(2)}\\ &c_p=\lim_{r\to\infty}p^{-r(n-1)}\left\vert V(\mathbb{Z}/p^r\mathbb{Z})\right\vert. \end{align} The product is taken over all primes $p$. This seems like a very neat expression, and can be seen that it gives the correct result for linear equations. However, it is not correct in general. Just looking at quadratic forms for $f$, the expression given by (2) is wrong. I do not have any good feeling as to where exactly this heuristic goes astray, and if it is possible to fix it. Maybe this approach and the reason that it does not quite work is well known. This is not an area in which I am any kind of expert, so maybe others on MathOverflow would be able to help?
For example, consider $f=x^2+y^2-z^2$, so that we are looking for primitive Pythagorean triples. Euclid's parameterization $(x,y,z)=(a^2-b^2,2ab,a^2+b^2)$ can be used to show that $\rho=\sqrt{2}\pi^{-2}\vert z\vert^{-1}$. However, on $V(\mathbb{R})$ we have $\Vert\nabla f\Vert = 2\sqrt{2}\vert z\vert$ and you can calculate $c_2=1$ and $c_p=1-p^{-2}$ for odd prime $p$. Using (2) would lead to $\rho=2\sqrt{2}\pi^{-2}\vert z\vert^{-1}$, which is out by exactly a factor of 2. If we look at Pythagorean quadruples $f=w^2+x^2+y^2-z^2$ instead, then we can calculate $c_p=(1-p^{-1})(1+2p^{-1}1_{\{p\equiv1{\rm\ mod\ }4\}}+p^{-2})$ for odd primes $p$, so the product in (2) is not unconditionally convergent.
Is there a known or, even, just conjectural expression for the asymptotic density $\rho$? And, is it possible to explain precisely how the heuristic used to derive (2) fails?