Two Dirac delta functions in an integral? For context, this is from a quantum mechanics lecture in which we were considering continuous eigenvalues of the position operator.
Starting with the position eigenvalue equation,
$$\hat{x}\,\phi(x_m, x)=x_m\phi(x_m,x)$$
where $x_m$ is the eigenvalue and $\phi(x_m, x)$ are the continuous eigenfunctions of the position operator $\hat x$. 
The professor wrote that $\phi(x_m, x)=\delta(x-x_m)$ and stated that this is because the eigenbasis is continuous.
But then he wrote the following
$$\begin{align}\int_{-\infty}^{\infty}\phi^*(x_m,x)\phi({x_m}',x)\,dx &=\int_{-\infty}^{\infty}\delta(x_m-x)\delta({x_m}'-x)\,dx \tag{1}\\ &=\color{red}{\delta(x_m-{x_m}')}\end{align}$$

I don't understand how the expression in $(1)$ can be equal to the expression in red. I do understand that $$\begin{align}\int_{-\infty}^{\infty}\delta(x_m-x)\phi({x_m}',x)\,dx &=\phi({x_m}',x_m)\tag{2}\\&=\color{#080}{\delta(x_m-{x_m}')}\end{align}$$
since the integral 'sifts' out the only value of $x$ where the argument of the Dirac delta function is zero (at $x_m$). I have applied this sifting property for one Dirac delta function (as in $(2)$).  
But I don't understand how this works when there are two Dirac deltas in the integrand, $(1)$.
By my logic, I think it should sift out each value, one at a time, so I think that $(1)$ should be 
$$\int_{-\infty}^{\infty}\delta(x_m-x)\delta({x_m}'-x)\,dx =\delta(x_m)+\delta({x_m}')$$
where the results of the integration are added, since the values $x_m$ and ${x_m}'$ are sifted out one after the other, depending on which of $x_m$ and ${x_m}'$ are larger (here I assumed ${x_m}'\gt x_m$).

Could someone please derive or explain why equation $(1)$ is true, or give me any hints to help me understand it.  
 A: The Dirac delta is a distribution, which means it acts on smooth functions. The Dirac delta itself is not a smooth function, which means that it cannot act on itself. The expression
$$
\int_{\mathbb R}\delta(x-x_1)\delta(x-x_2)\mathrm dx=\delta(x_1-x_2)\tag1
$$
is meaningless from a rigorous point of view. While it's true that it makes intuitive sense, making it rigorous requires a bit of care.
To this end, we shall regard the Dirac delta as the limit (in the sense of measures) of a mollifier:
$$
\delta_\epsilon(x):=\epsilon^{-1}\eta(x/\epsilon)\tag2
$$
where $\eta$ is an absolutely integrable function with unit integral.
With this,
$$
\lim_{\epsilon\to0^+}\int_{\mathbb R}f(x-y)\delta_\epsilon(x)\mathrm dx\equiv f(y)\tag3
$$
for good enough $f$. Note that this is just a convolution: $f\star\delta_\epsilon\to f$ as $\epsilon\to0^+$. This is a perfectly rigorous statement.
But note that $\delta_\epsilon$ is itself a smooth function, which means that we can apply this identity to $f=\delta_{\epsilon'}$:
$$
\lim_{\epsilon\to0^+}\int_{\mathbb R}\delta_{\epsilon'}(x-y)\delta_\epsilon(x)\mathrm dx\equiv \delta_{\epsilon'}(y)\tag4
$$
Finally, the formal limit $\epsilon'\to0^+$ yileds the identity in the OP (after the redefinition $y=x_2-x_1$ and the change of variables $x\to x-x_1$). Needless to say, this latter limit only makes sense in the sense of measures, which means you should regard it as a statement that both sides are equal when integrated over good-enough functions.
A: Distinguish two cases:


*

*if $x_m\ne x'_m$, the product of the deltas is always zero and so is the integral;

*if $x_m=x'_m$, you have a squared delta, which gives a divergent integral.
Hence, the value of the integral can be expressed as $\delta(x_m-x'_m)$.

For rigor, one should show that the integral of a squared delta is the "same infinity" as a single delta. In other terms
$$\int_{y\in\mathbb R}\int_{x\in\mathbb R}\delta^2(x)\,dx\,dy=1.$$
A: OP is asking:

Why is $ \int_{\mathbb{R}}\! \mathrm{d}y ~\delta(x-y) \delta(y-z)~=~\delta(x-z)~?$ 

Answer: Both the lhs. and the rhs. are informal notations 
$$ \iiint_{\mathbb{R}^3}\! \mathrm{d}x~\mathrm{d}y~\mathrm{d}z~ \delta(x-y) \delta(y-z) f(x,z)~=~u[f]~=~\iint_{\mathbb{R}^2}\! \mathrm{d}x~\mathrm{d}z~ \delta(x-z) f(x,z) $$
for the same distribution $u\in D^{\prime}(\mathbb{R}^2)$ defined as
$$u[f]~:=~ \int_{\mathbb{R}}\!\mathrm{d}y~f(y,y), \qquad f~\in~D(\mathbb{R}^2).$$
A: Intuitively, $\delta(x_m-x)$ is non-zero only when $x_m=x$. Likewise $\delta(x_m'-x)$ is non-zero only when $x_m'=x$. So the factor $\delta(x_m-x) \, \delta(x_m'-x)$ is non-zero only when $x_m=x=x_m'.$ After integration w.r.t. $x$ we still must have $x_m=x_m'$ which makes the result $\delta(x_m-x_m')$ quite natural.
More formally, recall the formula $\int_{-\infty}^{\infty} f(x) \, \delta(x_0-x) \, dx = f(x_0)$. Apply this to the given integral with $x_0=x_m'$ and $f(x) = \delta(x_m-x).$ Then the result $\delta(x_m-x_m')$ falls out.
A: $\delta(x)$ is never really well-defined by itself (at least not as a function). It's only defined when appearing in an integral (possibly, multiplied with another function). It is essentially syntactic sugar, and the “definition” is
$$ \int_\mathbb{R}\mathrm{d}x\ f(x)\cdot \delta(x-x_0) := f(x_0) $$
So, what the identity in question actually means is
$$
  \int_\mathbb{R}\mathrm{d}y\: f(y)\cdot \int_\mathbb{R}\mathrm{d}x\: \delta(y-x) \cdot \delta(x_m'-x)
      = f(x_m')
$$
for any function $f$.
Why would that be? Well, it's basically just a matter of changing the order of integration:
$$\begin{align}
  \int_\mathbb{R}\mathrm{d}y\: f(y)\cdot \int_\mathbb{R}\mathrm{d}x\: \delta(y-x) \cdot \delta(x_m'-x)
      =&   \int_\mathbb{R}\mathrm{d}y \int_\mathbb{R}\mathrm{d}x\: f(y)\cdot \delta(y-x) \cdot \delta(x_m'-x)
   \\ =&   \int_\mathbb{R}\mathrm{d}x \int_\mathbb{R}\mathrm{d}y\: f(y)\cdot \delta(y-x) \cdot \delta(x_m'-x)
   \\ =&   \int_\mathbb{R}\mathrm{d}x\: \delta(x_m'-x) \int_\mathbb{R}\mathrm{d}y\: f(y)\cdot \delta(y-x)
   \\ =&   \int_\mathbb{R}\mathrm{d}x\: \delta(x_m'-x) \int_\mathbb{R}\mathrm{d}y\: f(y)\cdot \delta(y-x)
   \\ =&   \int_\mathbb{R}\mathrm{d}x\: \delta(x_m'-x) f(x)
   \\ =&   \int_\mathbb{R}\mathrm{d}x\: \delta(x-x_m') f(x)
   \\ =&   f(x_m')
\end{align}$$
The above would be uncontroversial if $\delta$ were just a smooth, odd function with compact support. It's not in fact, but you can approximate it with such functions, so at least for continuous $f$ it's unproblematic.
A: Here's a derivation using calculus
$$
\int_{-\infty}^{\infty} \delta(x_m - x)\delta(x_m' - x)
$$
Integrating by parts
$$
= (\delta(x_m' - x) \int \delta(x_m - x')dx') |^\infty_{-\infty} - \int^\infty_{-\infty} \frac{d\delta(x_m' - x)}{dx} \int\delta(x_m-x')dx' dx
$$
Note that $\int \delta(x_m -x')dx' (x)$ is the heaviside step function, $H(x_m - x)$, which is 0 below $x_m$, $1/2$ at $x_m$, and $1$ above $x_m$. Thus we have
$$
= \lim_{x\rightarrow \infty} \delta(x_m' - x)H(x_m - x) - \lim_{x\rightarrow -\infty} \delta(x_m' - x)H(x_m - x) - 
\int^\infty_{-\infty} \frac{d\delta(x_m' - x)}{dx} H(x_m - x) dx
$$
The limits go to $0$, because the heaviside function takes on a finite value and the delta functions are zero as $x \rightarrow \infty$ and as $x \rightarrow -\infty$. Thus
$$
= - \int^\infty_{-\infty} \frac{d\delta(x_m' - x)}{dx} H(x_m - x) dx
$$
WLOG assume that $x_m \leq x_m'$. We note that the heaviside function is $0$ below $x_m$ and $1$ above $x_m$ (we can ignore the value $\frac{1}{2}$ at $x_m$ in integrating because that is a content zero set of values). Thus we have that
$$
= - \int^\infty_{x_m} \frac{d\delta(x_m' - x)}{dx} dx
$$
$$
= - (\delta(x_m' - x)|^\infty_{x_m}) = \delta(x_m' - x_m) = \delta(x_m - x_m')
$$
To address the assumption $x_m \leq x_m'$: we could have done integration by parts on the other delta function and gotten the same result with the assumption $x_m' \leq x_m$, so we cover all cases.
A: This might be the simplest way to think about it.
If we were working with discrete indices, $\sum_x\phi_{x_m,\,x}\delta_{x,\,x_m}=\phi_{x_m,\,x_m'}$ would be an equation in matrices, viz. $\phi\mathbb{I}=\phi$. The Dirac delta is an "identity matrix" on a vector space of uncountable dimension, where we integrate over $x$ instead of summing over it. The result you asked about is just $\mathbb{I}^2=\mathbb{I}$. In other words, the special case $\phi=\delta$ is legitimate; whether $\phi$ is a true function or a measure doesn't really matter.
A: If you take $u=x_m'-x$, then you have $\int_{-\infty}^{\infty}\delta(u+(x_m-x'_m))\delta(u)$.  Actually finding that rigorously is difficult, but if you use the heuristic that  $\int_{-\infty}^{\infty}f(u)\delta(u) = f(0)$, then taking $f(u)= \delta(u+(x_m-x'_m))$ and setting $u=0$ gets you $\delta(x_m-x'_m)$.
A: I want to refer you to this theory of divergent integrals.
Under this angle, long story short, the integral in question should be equal $\frac2\pi\delta(x_m-x'_m)$. There, $\delta(0)$ has a value with which we can manipulate.
$$\int_{-\infty}^\infty (\pi\delta(x))^2 dx=\int_{0^+}^\infty\frac2{x^2}dx=2\pi \delta(0)$$
Since this is different from the result of your lecturer, I suspect his answer did not predict any physical consequences that depend on that coefficient.
