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](http://www.hairer.org/notes/SPDEs.pdf) , 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. 



Regularity for Stochastic heat equation with additive noise in d=2
I hope this is an answer to your question. In any case it is too long for a comment.


With which space of distributions you are working depends on where you are solving the equation. Two natural choices are $[0,T]\times \mathbb{R}^d$ or $[0,T]\times \mathbb{T}^d$, the latter being a torus (i.e. periodic boundary conditions). Let us call $\mathcal{X}$ any of the above two, then white noise has a modification that is a Schwartz distribution:
$$
\xi \in \mathcal{D}(\mathcal{X}).
$$
Meaning that it is a distribution both in space and in time. Indeed the noise does not see the difference, and it can be proven that it lies in a space of negative Holder regularity. 
$$
\xi \in C^{-(d+1)/2-\epsilon} (\mathcal{X}), 
$$
for any $\epsilon >0.$ Now the rule of thumb is that solving the heat equation (i.e. convolution with the heat operator) makes you gain $+2$ in parabolic regularity, where parabolic stands for the fact that time regularity should be counted as double of the space regularity. In this sense $\xi$ should be of parabolic regularity:
$$
\xi \in C^{-(d+2)/2-\epsilon} (\mathcal{X})
$$
so by Schauder estimates (which I believe you can find in some of M. Hairers papers) integrating against the heat kernel (which we denote wich $\mathcal{I}$) gives you parabolic regularity:
$$
u = \mathcal{I}(\xi) \in C^{-d/2 +1-\epsilon} (\mathcal{X}).
$$
so that in dimension $1$ this is still function valued, while in dimension $2$ or greater it is distribution valued.


