When the variables $X, Y$ are independent, then the PDF of $Z = X + Y$ can be computed using convolutions: $$ f_Z(z) = \int_{-\infty}^{\infty} f_X(x)f_Y(z - x) dx $$
When the variables are dependent, apparently you can use $$ f_Z(z) = \int_{-\infty}^{\infty} f_{XY}(x, z - x) dx $$
I am wondering where the expression came from for the dependent case? It looks very similar to the independent case except you can't separate the joint distribution into marginals.
On
The characteristic function of $Z=X+Y$ is $$E[e^{iuZ}]=E[e^{iu(X+Y)}]=\int_{\mathbb{R}^2} e^{iu(x+y)}p_{X,Y}(x,y)dxdy$$ The pdf is found by inverse Fourier Transform. By exchanging integrals $$f_{X+Y}(z)=\frac{1}{2\pi}\int_{\mathbb{R}^2}\bigg(\int_{\mathbb{R}} e^{iu(x+y)}e^{-iuz}du\bigg)p_{X,Y}(x,y)dxdy=$$ $$=\int_\mathbb{R^2}\delta(z-(x+y))p_{X,Y}(x,y)dxdy$$ By symmetry of the Dirac delta $\delta(w)=\delta(-w)$ $$\int_\mathbb{R^2}\delta((x+y)-z)p_{X,Y}(x,y)dxdy=$$ $$(x+y)-z=0 \implies y=z-x $$ $$=\int_\mathbb{R}p_{X,Y}(x,z-x)dx$$
$$P(Z\leq z)=P(X+Y\leq z)=\int P(X+y\leq z|Y=y)p_Y(y)dy=$$ $$=\int_{\infty}^\infty\int^{z-y}_{-\infty}p_{X|Y}(x,y)p_y(y)dxdy=\int_{\infty}^\infty\int^{z-y}_{-\infty}p_{X,Y}(x,y)dxdy$$ $$f_Z(z)=\frac{d}{dz}P(Z\leq z)=\int_{-\infty}^\infty p_{X,Y}(z-y,y)dy$$ By using Leibniz integral rule.
On
From the cumulative distribution function: $$F_Z(z)= P(Z\le z ) = P(X+Y \le z) = P(Y \le z -X) = \int_{-\infty}^{\infty} \int_{-\infty}^{z-x} f_{XY}(x,y)\, dy\, dx $$
Now, $$f_Z(z) = \frac{d F_Z(z) }{ d z} = \int_{-\infty}^{\infty} \frac{d }{ d z} \left( \int_{-\infty}^{z-x} f_{XY}(x,y) dy\, \right) dx = \int_{-\infty}^{\infty} f_{XY}(x,z-x) dx $$
If $X$ and $Y$ are independent we can further write
$$f_Z(z) = \int_{-\infty}^{\infty} f_X(x) f_Y(z-x) dx $$
On
A pdf for a random vector $X:\Omega\to\Bbb R^n$ is any function $g:\Bbb R^n\to \Bbb R$ such that $E[h(X)]=\int_{\Bbb R^n}h(x)g(x)\,dx$ for, say, all bounded summable Borel functions $h:\Bbb R^n\to\Bbb R$.
So, if your random vector is $V=(X,Y)$, then, without too many ceremonies, you can use the change of variables $(t,s)= (x,y+x)$ and Fubini-Tonelli in the integral $$E[h(X+Y)]=\int_{\Bbb R^2} h(x+y)f_{X,Y}(x,y)\,dxdy=\int_{\Bbb R^2}h(s)f_{X,Y}(t,s-t)\,dsdt=\\=\int_{-\infty}^\infty f(s)\left(\int_{-\infty}^\infty f_{X,Y}(t,s-t)\,dt\right)\,ds$$
And there you have it.
There is a property for linear combinations that says if X, Y have joint pdf $f(x, y)$ and $Z=aX+bY+c$, then Z has pdf $$g(z)=\int_{-\infty}^\infty f\left(x, \frac{z-c-ax}b\right)\frac1{|b|}dx$$
For $Z=X+Y$, $$\begin{split}G(z)&=\int_{-\infty}^\infty\int_{-\infty}^{z-x}f(x,y)dydx\end{split}$$
Change of variables $w=x+y$
$$\begin{split}G(z)&=\int_{-\infty}^\infty\int_{-\infty}^{z}f(x,w-x)dwdx\\ &=\int_{-\infty}^z\int_{-\infty}^{\infty}f(x,w-x)dxdw\end{split}$$
Derivative with respect to $z$
$$\begin{split}g(z)&=\left[\frac{d}{dz}(z)\right]\int_{-\infty}^\infty f(x, z-x)dx-0\\ &=1\cdot\int_{-\infty}^\infty f(x, z-x)dx=\int_{-\infty}^\infty f(x, z-x)dx\end{split}$$