Find harmonic function s.t. $\int_{\mathbb R^d}|D^2 u|^2<\infty $ Find all harmonic functions $u$ defined in $\mathbb R^d$ ($d\geq 2$) s.t. 
$$\int_{\mathbb R^d}|D^2 u|^2<\infty .$$
I recall that $$|D^2 u|^2=\sum_{i=1}^d\sum_{j=1}^d \left(\frac{\partial^2 u}{\partial x_i\partial x_j}\right)^2.$$
My attempts
My idea was to set $h_i=\frac{\partial u}{\partial x_i}$ and remark that $$|D^2 u|^2=\sum_{i=1}^d|\nabla h_i|^2=\sum_{i=1}^d\nabla h_i\cdot \nabla h_i=\sum_{i=1}^d \text{div}(h_i\nabla h_i)-\sum_{i=1}^n (h_i \Delta h_i).$$
Since $u$ is harmonic, then $u\in \mathcal C^\infty $ and thus $\Delta h_i$ is well defined, but I can't conclude. I wanted to use divergence theorem, but since $\mathbb R^d$ is not bounded, it's not possible.
A second idea was to sue the fact that $u$ is harmonic, and thus
$$|D^2u|^2=2\sum_{i=1}^{n-2}\sum_{j=1}^{n-1}\left(\frac{\partial ^2 u}{\partial x_i\partial x_j}\right)^2,$$
and get $$\int_{\mathbb R^d}|D^2 u|^2<\infty \iff 2\sum_{i=1}^{n-2}\sum_{j=i+1}^{n-1}\int_{\mathbb R^d}\left(\frac{\partial ^2u}{\partial x_i\partial x_j}\right),$$
but still, I can't continue.
 A: Suppose that $u$ is a harmonic function on $\mathbb{R}^n$.  Then by the mean-value property we know that for each $x \in \mathbb{R}^n$ and $r>0$ we have that
$$
u(x) = \frac{1}{\omega_n r^{n}} \int_{ B(x,r)}  u(y) dy 
$$
where $\omega_n$ is the volume of the unit ball in $\mathbb{R}^n$.  From this and Cauchy-Schwarz we get the estimate
$$
|u(x)| \le \frac{1}{\omega_n r^n} \int_{B(x,r)} |u(y)| dy \le \frac{1}{\omega_n r^n} \left( \int_{B(x,r)} |u(y)|^2 dy \right)^{1/2} |B(x,r)|^{1/2} \\
=\frac{1}{\sqrt{\omega_n} r^{n/2}} \left( \int_{B(x,r)} |u(y)|^2 dy \right)^{1/2}
$$
Now, if $u$ is such that 
$$
\int_{\mathbb{R}^n} | u(y)|^2 dy = M^2 < \infty,
$$
then we in turn find that 
$$
|u(x)| \le \frac{1}{\sqrt{\omega_n} r^{n/2}} \left( \int_{B(x,r)} |u(y)|^2 dy \right)^{1/2} \le \frac{1}{\sqrt{\omega_n} r^{n/2}} \left( \int_{\mathbb{R}^n} |u(y)|^2 dy \right)^{1/2}  = \frac{M}{\sqrt{\omega_n} r^{n/2}}.
$$
Then for any fixed $x \in \mathbb{R}^n$ we can send $r \to \infty$ to get that 
$$
|u(x)| \le \lim_{r \to \infty}\frac{M}{\sqrt{\omega_n} r^{n/2}} =0,
$$
and hence we find that $u(x) =0$ for all $x$.  Thus the only harmonic function that is square-integrable over all of $\mathbb{R}^n$ is $0$.
Having established this, we can now answer your question.  Suppose now that $u$ is harmonic on $\mathbb{R}^n$ and $D^2 u$ is square-integrable.  Since $u$ is harmonic, it's smooth, and so we can apply arbitrarily many derivatives to the PDE $\Delta u =0$.  In particular we find that if $|\alpha|=2$ then
$$
0 = \partial^\alpha \Delta u = \Delta \partial^\alpha u,
$$
which means that all of the second-order partial derivatives are also harmonic.  Since
$$
\int_{\mathbb{R}^n} |\partial^\alpha u |^2 \le  \int_{\mathbb{R}^n} |D^2 u |^2 < \infty
$$
we find that $\partial^\alpha u$ is a harmonic function that is square-integrable.  By the above analysis we then know that $\partial^\alpha u =0$ for all multi-indices with $|\alpha |=2$.  In other words, we know that $D^2 u =0$, and so basic calculus tells us that $u$ is linear, i.e. $u(x) = a + b \cdot x$ for some $a \in \mathbb{R}$ and $b\in \mathbb{R}^n$.  It turns out that all such functions are trivially harmonic.
Thus a harmonic functions has $D^2 u$ square integrable if and only if it's linear.
