I would greatly appreciate any references were they study the stochastic equation in higher dimensions: $u_{t}=\Delta u+f$ in great detail, especially in dimension 2.

In Hairer's Spde notes , he mentions that the solutions will not be function-valued, and I am curious what the particular regularity will be with ,say, $u_{0}\equiv 0$.

Q: Over which space they will be distributions?

For $u_{0}\equiv 0$ the covariance for $u(x,t)$ is

$$E[u(x,t)u(y,s)]=2^{-(n+1)}\int_{[|s-t|,s+t]}l^{-n/2}exp(-|x-y|^{2}/4l)dl=:C(t,s,x,y)$$

and one idea is to find a Gelfand-triple so that using Gross results, we can define a Gaussian process with this covariance. To obtain some space of measures s.t. $<u,\rho>$ is a Gaussian process with covariance

$$E[<u,\rho><u,\mu>]=\int \int C(t,s,x,y)d\rho(x,t) d\mu(y,s). $$

This is in spirit of the mathematical construction for white noise and Gaussian free field.

Update: You can find more information as described in the answer below in "Stochastic PDEs, Regularity Structures, and Interacting Particle Systems" Theorem 2.8.