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π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\|^{-β}\mathrm{d} x$, we obtain a multiple phase transition for the CLT for $U(t)$ from $β\in(0\,,1)$ to $β=1$ to $β\in(1\,,d\wedge 2)$.