Is square of Delta function defined somewhere? 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
 A: There are whole theories in microlocal analysis that deal with the issues here, I believe. Some heuristics are that the "singular support" of a distribution controls what it can be multiplied by in a naive sense (distributions with a disjoint singular support). So squaring the delta function is the first bad case - whatever the singular support means, it must be the set containing 0 for the delta function. Need more heuristics.
One insight is that one dimension may be too few to show the real picture. "Microlocal" tends to mean localising in (co)tangential directions, and one dimension offers only two. Hyperfunctions in the case of one dimension make something of this by considering the real line as the boundary of the upper half complex plane. I.e up is not the same as down. Boundary values of functions holomorphic in the upper half plane have a candidate for the delta function analogue: take 1/z. No problem squaring that. More of a problem saying what this analogy means that is worth anything. Mikio Sato did that. Now I shall be quiet, because this is probably already wrong enough.
A: 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.
A: The extent to which multiplication of distributions is defined was examined by Richards & Youn and some of the results are in their short and fairly elementary joint book on distributions.
One can multiply something fairly exotic like the third derivative of the delta function by a very well-behaved function; that much everybody knows.  But I think they had a result that as one factor becomes progressively less well-behaved the other must become more well-behaved in order to make multiplication possible.  I don't recall the details.  But I'm pretty sure theirs is not the last word on the subject.
A: I'd like to point out that several of the concepts mentioned here are explained on the nLab:
multiplication of distributions, 
while several are missing and the parts on microlocal analysis and hyperfunctions could use some help .
A: The theory of distributions and operations on them are generally only useful in so far as they extend the operations on smooth functions.  If you look in Hörmander, there is a criterion in terms of wavefront sets which is very useful (mentioned by others), and you'll also notice that the wavefront sets of $\delta$ and $\delta$ collide.  The reason you can't square the delta-function is that when you approximate it by smooth functions, there is no unique limit.  If you wanted to restrict to a smaller space of test functions, you would clearly have to consider test functions which vanish at the origin in some way.  But do you have a particular purpose in mind for this question?
EDIT:  Sorry -- this was supposed to be a comment, not an answer.
A: Denis Serre's answer is just perfect. Let me add a couple of examples of distributions that can be squared:
(1) With $H$ the Heaviside function, define $\operatorname{Log}(x+i0)=\ln(\vert x\vert)+i\pi H(-x) $
and
$$T_1=\frac{1}{x+i0}=\frac{d}{dx}(\operatorname{Log}(x+i0)) = \operatorname{pv}\frac 1{x}-i\pi \delta_0(x).$$
It is easy to see that
$
WF T_1=[0]\times (0,+\infty),
$
so that $WF T_1+WF T_1$ does not meet $0.$
Then there is no difficulty to define $T^2$ say as
$$
\langle T^2,\phi\rangle=\lim_{\varepsilon\rightarrow 0_+} \int\frac{\varphi(x) \, dx}{(x+i\varepsilon)^2}.
$$
(2) Let us consider a smooth hypersurface $\Sigma$
of $\mathbf R^d$ defined by the equation $f(x)=0$ with a smooth $f$ such that $df\not=0$ at $f=0$ and let $\delta_\Sigma$ be the Euclidean measure on $\Sigma$. Then
$$
T_2=pv\frac{1}{f}-i\delta_\Sigma
$$
can be squared. The reason is the same than for the previous example, since
$
WF T_2$
is the positive conormal of $\Sigma$. A point $(x,\xi)\in WF T_2$ iff
$$
x\in \Sigma\quad \xi =\lambda df(x) \text{ with $\lambda >0$}.
$$
Then of course, if $(x,\xi_j)$, $j=1,2$ are both in $WF T_2$ then
$$
\xi_1+\xi_2\not=0.
$$
A: $\delta_0$ vanishes identically on the space of test functions you've defined.  So it's not surprising that its square is well-defined:   $0\cdot 0 = 0$.
I suspect you'll have a much harder time defining $\delta_0^2$ on test functions which don't vanish at $0$.
A: I've seen the idea of it used in image processing for denoising; the total variation energy
$E_{TV}(f)  = \displaystyle\int\left(|\nabla f| +(f-u)^2\right)$
is generally used instead of Tikhonov regularization
$E_{Tikhonov}(f)  = \displaystyle\int\left(|\nabla f|^2 +(f-u)^2\right)$
as the latter never has a discontinuous solution (since the integral would be infinite).
I don't remember how rigorously this idea was developed -  "Mathematical Problems in Image Processing" by Aubert and Kornprobst was the textbook I used at the time, but there are probably some more recent references in the field.
A: Check out the papers by Accardi and Boukas
[added by S. Carnahan: The relevant part of their first ArXiv paper is that they regularize powers of $\delta$, not by setting $\delta^n(x) = c_n \delta(x)$ for some real $c_n$ (which they find to work poorly for their purposes), but by using two variables and setting $\delta(t-s)^n = \delta(s)\delta(t-s)$ for all $n \geq 2$.  This definition allows them to define certain representations of Lie algebras by Fock space methods.  As far as I can tell, this does not yield a workable definition of $\delta^2$ for analysis on the real line.]
