Sum of Random Variables Involving Dirac's Delta Studying Statistical Mechanics I came upon the following statement: Let $X$ and $Y$ be two random variables with joint density $f_{XY}$. The density of $Z=X+Y$  is given by:
$$f_Z(z)=\int d\textrm{x}d\textrm{y}f_{XY}(x,y)\delta(z-x-y)$$
where $f_{XY}(\cdot,\cdot)$ denotes the joint density of $X$ and $Y$ and $\delta(\cdot)$ is Dirac's Delta "function".
As a math student the use of Dirac's delta here caught me unexpectedly. Yet I tried to sketch a proof that tried to convince me this was true. I proceeded as follows:
$$f_Z(z)=\dfrac{\partial}{\partial z}P(X+Y\leq z)=\dfrac{\partial}{\partial z}\int_{x+y\leq z}d\textrm{x}d\textrm{y}f_{XY}(x,y)=\dfrac{\partial}{\partial z}\int_{-\infty}^\infty d\textrm{x}\int_{-\infty}^{z-x}d\textrm{y}f_{XY}(x,y)$$
I now made use of Heaviside's step function $\theta(\cdot)$, that is $\theta(x)=1$ if $x\geq 0$ and $\theta(x)=0$ if $x<0$ to rewrite the last integral as follows:
$$\dfrac{\partial}{\partial z}\int_{-\infty}^\infty d\textrm{x}\int_{-\infty}^{\infty}d\textrm{y}f_{XY}(x,y)\theta(z-x-y)=\int_{-\infty}^\infty d\textrm{x}\int_{-\infty}^{\infty}d\textrm{y}f_{XY}(x,y)\dfrac{\partial}{\partial z}\theta(z-x-y)$$
Lastly I used the fact that $\dfrac{d}{d\textrm{x}}\theta(x)=\delta(x)$ to write
$$\int_{-\infty}^\infty d\textrm{x}\int_{-\infty}^{\infty}d\textrm{y}f_{XY}(x,y)\dfrac{\partial}{\partial z}\theta(z-x-y)=\int_{-\infty}^\infty d\textrm{x}\int_{-\infty}^{\infty}d\textrm{y}f_{XY}(x,y)\delta(z-x-y) $$
which is esentially what I wanted to proof. I do still have concerns of the rigour in this procedure. I am not used to dealing with Heaviside or Dirac's Delta so I do not know if this reasoning is correct at all.
 A: You can do it from the definition of the 1D Dirac delta by changing variables in the integration (for each fixed $z$) to a coordinate system with $u=z-x-y$ as one of the coordinates and, say, $v=x$ as the other one. Then $\frac{\partial u}{\partial x}=-1,\frac{\partial u}{\partial y}=-1,\frac{\partial v}{\partial x}=1,\frac{\partial v}{\partial y}=0$ so the absolute Jacobian determinant of the forward transformation is $1$ and thus the absolute Jacobian determinant of the backward transformation is also $1$. Moreover $(x,y) \mapsto (z-x-y,x)$ maps onto the whole plane.
So then
$$\int_{-\infty}^\infty \int_{-\infty}^\infty f_{X,Y}(x,y) \delta(z-x-y) dx dy = \int_{-\infty}^\infty \int_{-\infty}^\infty f_{X,Y}(v,z-u-v) \delta(u) du dv \\ = \int_{-\infty}^\infty f_{X,Y}(v,z-v) dv$$
where in the second line I used the definition of $\delta(u)$. This effectively defines the delta function of a function of the independent variable(s) as "whatever makes change of variable make sense" (which is indeed how mathematicians define it, and it is indeed self-consistent).
Now note that if you integrate that in $z$ over whatever (measurable) subset of $\mathbb{R}$ you like, you get the integral of the joint density over the region where $x+y=z$, which is what you want to see.
Intuitively, the delta function is restricting the integration to the part of the space where $z-x-y=0$ which is what you want, but there is a question of whether maybe you need some additional coefficient to counterbalance the change of variable. This is what the Jacobian took care of for us here. In a different scenario (for a silly example, if you had written it with $\delta(2z-2x-2y)$ instead), the Jacobian would not necessarily have been $1$ and some coefficient would have appeared.
