It is well-known that one cannot, in general, multiply distributions. In particular, it is not clear what does it mean to define a square delta distribution. However, as discussed here, expressions like: \begin{equation} \int dz \delta(x-z)\delta(y-z) = \delta(x-y) \tag{1}\label{1} \end{equation} have mathematical meaning.
When I learned distributions, I remember that the Dirac delta was defined as follows. Take a smooth function $f$ on $\Omega \subseteq \mathbb{R}^{n}$ and fix $x \in \Omega$. Then the Dirac delta distribution $\delta_{x}$ is the distribution defined by: $$\langle \delta_{x},f\rangle := f(x) \tag{2}\label{2}$$ The right hand side of (\ref{2}) is a number, so it does not make sense to apply a new Dirac delta to this object. However, we can think of $f$ itself as a function $\Omega \ni x \mapsto g(x) :=\langle\delta_{x},f\rangle \equiv f(x)$, so it would make sense to think of: $$\langle \delta_{y},g\rangle = f(y)\tag{3}\label{3}$$
Relation (\ref{3}) seems to be well-defined as, although it differs from applying the Dirac delta again in (\ref{2}), the difference seems sutile.
I want to know: is (\ref{3}) what is really going on in (\ref{1})? I think the integral representation of the Dirac delta ends up making manipulations easier, but I wanted to give meaning to (\ref{1}) by using the Dirac delta as usually is done in distribution theory courses (as far as I know, the integral representation of the delta is more common in the physics literature).
