Harmonic function and interation Let $f:\mathbb{R}^k \to \mathbb{R}$ be a bounded harmonic function (i.e. $f$ is of class $C^2$ and $\sum_{r=1}^k\frac{\partial^2f}{\partial y^2_r}=0$) and $(x_1,...,x_k) \in \mathbb{R}^k.$
Prove that $$f(x_1,...,x_k)=\frac{1}{(\sqrt{2\pi})^k}\int_{\mathbb{R}^k}f(x_1+y_1,...,x_k+y_k)e^{-\frac{1}{2}\sum_{r=1}^ky_r^2}dy_1...dy_k.$$
We note that $$\int_{\mathbb{R}^k}f(x_1+y_1,...,x_k+y_k)e^{-\frac{1}{2}\sum_{r=1}^ky_r^2}dy_1...dy_k=\int_{\mathbb{R}^k}f(y_1,...,y_k)e^{-\frac{1}{2}\sum_{r=1}^k(y_r-x_r)^2}dy_1..dy_k,$$
how to use the fact that $f$ is harmonic ?
 A: Note: I've used $n$ instead of $k$ - sorry this is a force of habit.
Since $u$ is harmonic it satisfies the mean value property: $$u(x) =\frac 1{ n \omega_n r^{n-1}}\int_{\partial B_r(x)}u(y)\,dS \qquad \text{for all }r>0$$ where $\omega_n$ is the volume of the $n$-ball with radius 1. Hence,$$n \omega_n r^{n-1}u(x)=\int_{\partial B_r(x)}u(y)\,dS \qquad \text{for all }r>0.$$  Multiply through by $e^{-\frac 12r^2}$ on both side then integrate from $0$ to $\infty$. On the left hand side you get $$n \omega_n \int_0^\infty u(x) r^{n-1} e^{-\frac 12r^2} \, dr =n \omega_n u(x)\int_0^\infty  r^{n-1} e^{-\frac 12r^2} \, dr = n \omega_n u(x) \frac{\Gamma(n/2)}{2^{1+n/2}}. $$  It is known that $\omega_n = \frac{\pi^{n/2}}{\Gamma(n/2+1)} $. Hence the LHS is: $$nu(x) \frac{\pi^{n/2}}{\Gamma(n/2+1)} \cdot\frac{\Gamma(n/2)}{2^{1+n/2}}=u(x) (2\pi)^{n/2} $$ using properties of the gamma function. On the right hand side you get $$ \int_0^\infty\int_{\partial B_r(x)}u(y)e^{-\frac12 r^2}\,dSdr = \int_{\mathbb R^n }u(y)e^{-\frac12\vert y-x \vert^2}\,dy$$ via polar coordinates. Thus, $$ u(x) = \frac 1 {(2\pi)^{n/2}}\int_{\mathbb R^n }u(y)e^{-\frac12\vert y-x \vert^2}\,dy.$$
