  Suppose that $\{u(t\,, x)\}_{t >0, x \in\mathbb{R}^d}$ is the solution to a
$d$-dimensional parabolic Anderson model with delta initial condition and
driven by a Gaussian noise that is white in time and has a spatially
homogeneous covariance given by a nonnegative-definite measure $f$ which
satisfies Dalang's condition. Let $\boldsymbol{p}_t(x):=(2\pi
t)^{-d/2}\exp\{-\|x\|^2/(2t)\}$ denote the standard Gaussian heat kernel on
$\mathbb{R}^d$. We prove that for all $t>0$, the process $U(t):=\{u(t\,,
x)/\boldsymbol{p}_t(x): x\in \mathbb{R}^d\}$ is stationary using Feynman-Kac's
formula, and is ergodic under the additional condition $\hat{f}\{0\}=0$, where
$\hat{f}$ is the Fourier transform of $f$. Moreover, using Malliavin-Stein
method, we investigate various central limit theorems for $U(t)$ based on the
quantitative analysis of $f$. In particular, when $f$ is given by Riesz kernel,
i.e., $f(\mathrm{d} x) = \|x\|^{-\beta}\mathrm{d} x$, we obtain a multiple
phase transition for the CLT for $U(t)$ from $\beta\in(0\,,1)$ to $\beta=1$ to
$\beta\in(1\,,d\wedge 2)$.
