Suppose that ${u(t,, x)}_{t >0, x inmathbb{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 Dalangs condition. Let $boldsymbol{p}_t(x):=(2pi 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): xin mathbb{R}^d}$ is stationary using Feynman-Kacs 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 $betain(0,,1)$ to $beta=1$ to $betain(1,,dwedge 2)$.