We consider a system of d non-linear stochastic heat equations in spatial dimension k ≥ 1, whose solution is an R d -valued random field u = {u(t , x), (t, x) ∈ R + × R k }. The d-dimensional driving noise is white in time and with a spatially homogeneous covariance defined as a Riesz kernel with exponent β, where 0 < β < (2 ∧ 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) → (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 >