Don't understand the integral over the square of the Dirac delta function

Well, the Dirac delta function $\delta(x)$ is a distribution, also known as a generalized function.

One can e.g. represent $\delta(x)$ as a limit of a rectangular peak with unit area, width $\epsilon$, and height $1/\epsilon$; i.e.

$$\tag{1} \delta(x) ~=~ \lim_{\epsilon\to 0^+}\delta_{\epsilon}(x), $$ $$\tag{2} \delta_{\epsilon}(x)~:=~\frac{1}{\epsilon} \theta(\frac{\epsilon}{2}-|x|) ~=~\left\{ \begin{array}{ccc} \frac{1}{\epsilon}&\text{for}& |x|<\frac{\epsilon}{2}, \\ \frac{1}{2\epsilon}&\text{for}& |x|=\frac{\epsilon}{2}, \\ 0&\text{for} & |x|>\frac{\epsilon}{2}, \end{array} \right. $$

where $\theta$ denotes the Heaviside step function with $\theta(0)=\frac{1}{2}$.

The product $\delta(x)^2$ of the two Dirac delta distributions does strictly speaking not$^1$ make mathematical sense, but for physical purposes, let us try to evaluate the integral of the square of the regularized delta function

$$\tag{3} \int_{\mathbb{R}}\! dx ~\delta_{\epsilon}(x)^2 ~=~\epsilon\cdot\frac{1}{\epsilon}\cdot\frac{1}{\epsilon} ~=~\frac{1}{\epsilon} ~\to~ \infty \quad \text{for} \quad \epsilon~\to~ 0^+. $$

The limit is infinite as Griffiths claims.

It should be stressed that in the conventional mathematical theory of distributions, the eq. (2.95) is a priori only defined if $f$ is a smooth test-function. In particular, it is not mathematically rigorous to use eq. (2.95) (with $f$ substituted with a distribution) to justify the meaning of the integral of the square of the Dirac delta distribution. Needless to say that if one blindly inserts distributions in formulas for smooth functions, it is easy to arrive at all kinds of contradictions! For instance,

$$ \frac{1}{3}~=~ \left[\frac{\theta(x)^3}{3}\right]^{x=\infty}_{x=-\infty}~=~\int_{\mathbb{R}} \!dx \frac{d}{dx} \frac{\theta(x)^3}{3} $$ $$\tag{4} ~=~\int_{\mathbb{R}} \!dx ~ \theta(x)^2\delta(x) ~\stackrel{(2.95)}=~ \theta(0)^2~=~\frac{1}{4}.\qquad \text{(Wrong!)} $$

--

$^1$ We ignore Colombeau theory. See also this mathoverflow post.

You need nothing more than your understanding of $$ \int_{-\infty}^\infty f(x)\delta(x-a)dx=f(a) $$ Just treat one of the delta functions as $f(x)\equiv\delta(x-\lambda)$ in your problem. So it would be something like this: $$ \int\delta(x-\lambda)\delta(x-\lambda)dx=\int f(x)\delta(x-\lambda)dx=f(\lambda)=\delta(\lambda-\lambda) $$ So there you go.

Suppose I want to show $$\int \delta(x-a)\delta(x-b)\; dx = \delta(a-b) $$ To do that , I need to show $$\int g(a)\int \delta(x-a)\delta(x-b) \;dx \;da = \int g(a)\delta(a-b)\; da$$ for any function $g(a)$. \begin{align}\textrm{LHS}& = \int \int g(a) \delta(x-a)\;da \ \delta(x-b) \;dx\\ &=\int g(x)\delta(x-b)\;dx \\&=g(b) \end{align} But $\textrm{RHS}$ clearly $=g(b)$ too.

The result follows putting $a=b=\lambda$