  In this paper, we establish lower and upper Gaussian bounds for the
probability density of the mild solution to the stochastic heat equation with
multiplicative noise and in any space dimension. The driving perturbation is a
Gaussian noise which is white in time with some spatially homogeneous
covariance. These estimates are obtained using tools of the Malliavin calculus.
The most challenging part is the lower bound, which is obtained by adapting a
general method developed by Kohatsu-Higa to the underlying spatially
homogeneous Gaussian setting. Both lower and upper estimates have the same
form: a Gaussian density with a variance which is equal to that of the mild
solution of the corresponding linear equation with additive noise.
