Let u be an harmonic function in $\mathbb{R}^n$ and g be a distribution function with compact support. I need to show that $u*g$ is an harmonic function.
Idea of proof
I start by considering a function $u:\mathbb{C}\rightarrow\mathbb{R}$ to be harmonic and bounded $\Rightarrow$ u is constant. Since $u$ is harmonic, we can use the mean value property. Let $a,b\in\mathbb{C}$, then we can show that $u(a)=u(b)$ by choosing a constant $R>0$ big enough so that: $u(j)=\frac{1}{\pi R^2}\int\int u(x,y)dxdy, \forall j\in\{a,b\}$. Then, $u(a)-u(b)=\frac{1}{\pi R^2}(\int\int_{D(a,R)\text{\ }D(b,R)}u(x,y)dxdy-\int\int_{D(b,R)\text{\ }D(a,R)}u(x,y)dxdy)$. Letting $R\rightarrow\infty$, $u(a)-u(b)=2M \times\{\text{smth goint to }0\}\Rightarrow u(a)=u(b)$, where $M$ is the upper bound for u.
Now, I know a result for the Poisson kernel via Laplace operator, but I don't think this is applicable in the case. Do you have any hints?
If $g$ is an honest function for which $u*g$ is defined as a function ($C_c^\infty$, for example) then you can differentiate the convolution under the integral sign to show that $\Delta(u*g) = 0$: the point is that you can let the Laplacian attack $u$ only, which (since $u$ is harmonic) annihilates the convolution. For more general $g$ you may need to demonstrate this identity in the sense of distributions, but the strategy is the same: pass the Laplacian from the test function to $u$ via integration by parts.