everyone. I am wondering whether anyone knows if the square of Dirac Delta function is defined somewhere. 

In the beginning, this question might look strange. But by restricting the space of the test functions, I think it is still possible. For example, in order to make sense of $\delta_0^2$, we can think that it is the limit of $\frac{e^{-x^2/t}}{2\pi t}$ as $t\rightarrow 0_+$. Now choose the test function $f(x)=x^2$. It is clear that
$$
\int_{-\infty}^{\infty} x^2 \frac{e^{-x^2/t}}{2\pi t} d x = \frac{1}{2\sqrt{\pi t}} \int_{-\infty}^{\infty} x^2 \frac{e^{-x^2/t}}{\sqrt{\pi t}} d x = \frac{1}{2\sqrt{\pi t}} \cdot \frac{t}{2} = \frac{\sqrt{t}}{4\sqrt{\pi}}\;.
$$
Then let $t$ tend to $0$, we get $\langle\delta_0^2,f\rangle=0$ in this case. So we can restrict, for example, all test functions tend to 0 at the speed no less than $x^2$.

I don't want to invent the whole stuff if it already exists. Otherwise, I might take care of the every details. Thank you in advance for any hints.

EDIT:
Here are some references that I found to be useful.


Mikusiński, J. On the square of the Dirac delta-distribution. (Russian summary) 
Bull. Acad. Polon. Sci. Sér. Sci. Math. Astronom. Phys. 14 1966 511–513. 
44.40 (46.40) 
Ta Ngoc Tri, The Colombeau theory of generalized functions Master thesis, 2005


When L. Schwartz "invented" distributions (actually, he only invented the mathematical theory as a part of functional analysis, because distributions were already used by physicists), he proved incidentally that it is impossible to define a product in such a way that distributions form an algebra with acceptable topological properties. What is possible is to define the product of distributions when their wave front sets do not meet. This is why $fT$ makes sense if $T$ is a distribution and $f$ is $C^\infty$, for instance, because the front set of $f$ is void. But you can also multiply that way genuine distributions; for instance in $\mathbb R^2$, 
$$(1)\qquad\delta_{x=0}=\delta_{x_1=0}\delta_{x_2=0}.$$ 

J.-F. Colombeau invented in the 70's an algebra of generalized functions, which has something to do with distributions. But each distribution has infinitely many representatives in the algebra, and you have to play with the equality and a "weak equality" (or "association"). I don't know of an example where this tool solved an open problem. In Colombeau's algebra, the square of $\delta_0$ makes sense, but is highly non unique.

Edit (May 2020). I'd like to share the following generalization of identity (1) above, which I found in developing my theory of Divergence-free positive symmetric tensor. In ${\mathbb R}^d$, consider the one-dimensional Lebesgue measure ${\cal L}_j$ along the $j$-th axis, for $1\le j\le d$. Then
$$({\cal L}_1\cdots{\cal L}_d)^{\frac1{d-1}}=\delta_{x=0}.$$
There are a lot of reasons why this equality makes sense and is valid. For instance, if you approach ${\cal L}_j$ by $(2\epsilon)^{1-d}dx|_{K_j(\epsilon)}$ where $K_j(t)=(-\epsilon,\epsilon)^{d-1}\times {\mathbb R}\vec e_j$, then the left-hand side equals $(2\epsilon)^{-d}dx|_{(-\epsilon,\epsilon)^d}$, which approaches the Dirac at the origin. There is an analogous identity when the orthogonal axes are replaced by an arbitrary list of $d$ axes; then the right-hand sides is $C\delta$, where the constant $C$ is computed by solving a case of Minkowski's Problem.
