Two Dirac delta functions in an integral?

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.

$\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.

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.