  We consider a system of $d$ non-linear stochastic heat equations in spatial
dimension $k \geq 1$, whose solution is an $\R^d$-valued random field $u=
\{u(t\,,x),\, (t,x) \in \R_+ \times \R^k\}$. The $d$-dimensional driving noise
is white in time and with a spatially homogeneous covariance defined as a Riesz
kernel with exponent $\beta$, where $0<\beta<(2 \wedge k)$. The non-linearities
appear both as additive drift terms and as multipliers of the noise. Using
techniques of Malliavin calculus, we establish an upper bound on the two-point
density, with respect to Lebesgue measure, of the $\R^{2d}$-valued random
vector $(u(s,y),u(t,x))$, that, in particular, quantifies how this density
degenerates as $(s,y) \to (t,x)$. From this result, we deduce a lower bound on
hitting probabilities of the process $u$, in terms of Newtonian capacity. We
also establish an upper bound on hitting probabilities of the process in terms
of Hausdorff measure. These estimates make it possible to show that points are
polar when $d > \frac{4+2k}{2-\beta}$ and are not polar when $d <
\frac{4+2k}{2-\beta}$. %(if $\frac{4+2k}{2-\beta}$ is an integer, then the case
$d=\frac{4+2k}{2-\beta}$ is open). In the first case, we also show that the
Hausdorff dimension of the range of the process is $\frac{4+2k}{2-\beta}$ a.s.
