In an exam, my professor gave the following exercise:
State and prove the mean value theorem for harmonic functions. Let $H$ be a harmonic function on $\mathbb{R}^n$. Show all of its dilations $H^\tau$ are harmonic. Here, $H^\tau(x)=H(\tau x)$ for all $\tau>0$ and all $x\in\mathbb{R}^n$. Let $P$ then be the function defined by:
$$P(x,y):=\frac{1-|x|^2|y|^2}{(1-2x\cdot y+|x|^2+|y|^2)^{\frac n2}},$$
for all $(x,y)\in\mathbb{R}^n\times\mathbb{R}^n$ for which the denominator is nonzero. Show $P$ is an extension of the Poisson kernel $P_0$ of the unit sphere in $\mathbb{R}^n$; in particular, show that:
$$P(x,y)=P_0(|x|y,x/|x|).$$
Deduce that for fixed $x$, the function $y\mapsto P(x,y)$ is hamonic.
My problem is the last part. Deducing from what was proved before seems to imply checking that $P_0(x,y)$ is harmonic in the first entry when the second one is fixed, since $P$ is a dilation of $P_0(\cdot,x/|x|)$. I first proved a general formula for $\Delta|x|^k$:
\begin{align*} \partial_j|x|^k={}&k|x|^{k-1}\partial_j|x|=k|x|^{k-2}x_j; \\ \partial_j^2|x|^k={}&k\partial_j(|x|^{k-2}x_j)=k[|x|^{k-2}+x_j\partial_j|x|^{k-2}]={} \\ {}={}&k[|x|^{k-2}+x_j(k-2)|x|^{k-4}x_j]; \\ \Delta|x|^k={}&\sum_{j=1}^n\partial_j^2|x|^k=k[n|x|^{k-2}+(k-2)|x|^{k-4}\sum_1^nx_j^2]={} \\ {}={}&k|x|^{k-2}[n+k-2]. \end{align*}
Then I remarked that (hopefully getting it right):
$$\Delta(uv)=u\Delta v+v\Delta u+\nabla u\cdot\nabla v.$$
Hence, the $x$-laplacian of the Poisson kernel would be:
\begin{align*} \Delta_x\left(\frac{1-|x|^2}{|x-y|^n}\right)={}&\frac{\Delta(1-|x|^2)}{|x-y|^n}+\nabla(1-|x|^2)\cdot\nabla(|x-y|^{-n})+\Delta_x|x-y|^{-n}\cdot(1-|x|^2)={} \\ {}={}&-\frac{2|x|^0[n+2-2]}{|x-y|^n}-2x\cdot(-n)|x-y|^{-n-1}\cdot\nabla|x-y|+{} \\ &{}+(1-|x|^2)(-n)|x-y|^{-n-2}[n-n-2]={} \\ {}={}&\frac{2n}{|x-y|^{n+2}}\left(-|x-y|^2+|x-y|x\cdot\frac{x-y}{|x-y|}+1-|x|^2\right)={} \\ {}={}&\frac{2n}{|x-n|^{n+2}}(1-|x-y|^2-x\cdot y), \end{align*}
which doesn't seem to be zero everywhere, not even where $x\neq y$. I then tried writing out averages of $P$ on balls, but went nowhere. So did I do something wrong up there? How do I do this exercise?
Update
As the answer shows, I found the mistake. But I am now wondering: is there a less brute-force-calculation-y way to do this?
I just found out the mistake:
$$\Delta(uv)=\Delta u\cdot v+u\Delta v+{\color{red}2}\nabla u\cdot\nabla v.$$
Redoing the calcs with this knowledge, and expanding the $|x-y|^2$ at the end, gave me:
\begin{align*} \Delta_x\left(\frac{1-|x|^2}{|x-y|^n}\right)={}&\frac{\Delta(1-|x|^2)}{|x-y|^n}+2\nabla(1-|x|^2)\cdot\nabla(|x-y|^{-n})+\Delta_x|x-y|^{-n}\cdot(1-|x|^2)={} \\ {}={}&-\frac{2|x|^0[n+2-2]}{|x-y|^n}-4x\cdot(-n)|x-y|^{-n-1}\cdot\nabla|x-y|+{} \\ &{}+(1-|x|^2)(-n)|x-y|^{-n-2}[n-n-2]={} \\ {}={}&\frac{2n}{|x-y|^{n+2}}\left(-|x-y|^2+2|x-y|x\cdot\frac{x-y}{|x-y|}+1-|x|^2\right)={} \\ {}={}&\frac{2n}{|x-y|^{n+2}}(1-|x-y|^2-2x\cdot y+|x|^2)=\frac{2n}{|x-y|^{n+2}}(1-|y|^2), \end{align*}
which is 0 for $y\in\mathbb{S}^{n-1}$. In $P$, we have $P_0$ with a modulus-1 second entry, $\frac{x}{|x|}$, which should be enough to conclude.