[proofplan]
We compute the law of $X+Y$ on an arbitrary Borel set $B\subset\mathbb R$. Independence identifies the joint law of $(X,Y)$ with the product of the marginal laws, hence with the density $(x,y)\mapsto f_X(x)f_Y(y)$ on $\mathbb R^2$. Tonelli's theorem and [translation invariance](/theorems/4911) of one-dimensional [Lebesgue measure](/page/Lebesgue%20Measure) then rewrite the integral over the set $\{(x,y):x+y\in B\}$ as an integral over $B$ of the convolution function. Since this holds for every Borel set $B$, that convolution is a density for $X+Y$, after modifying it on a null set if necessary.
[/proofplan]
[step:Identify the joint law of $(X,Y)$ from independence]
Let $\mu_{(X,Y)}$ denote the law of the random vector
\begin{align*}
(X,Y):(\Omega,\mathcal F)\to(\mathbb R^2,\mathcal B(\mathbb R^2)).
\end{align*}
Thus, for every $A\in\mathcal B(\mathbb R^2)$,
\begin{align*}
\mu_{(X,Y)}(A)=\mathbb P((X,Y)\in A).
\end{align*}
Let $\mathcal L^2$ denote two-dimensional Lebesgue measure on $\mathbb R^2$.
Define the measure $\nu$ on $(\mathbb R^2,\mathcal B(\mathbb R^2))$ by
\begin{align*}
\nu(A)=\int_A f_X(x)f_Y(y)\,d\mathcal L^2(x,y).
\end{align*}
The integrand $(x,y)\mapsto f_X(x)f_Y(y)$ is Lebesgue measurable and non-negative because $f_X$ and $f_Y$ are Lebesgue-measurable non-negative densities.
We show that $\mu_{(X,Y)}=\nu$. Let $A,C\in\mathcal B(\mathbb R)$. Since $X$ and $Y$ are independent, the events $\{X\in A\}$ and $\{Y\in C\}$ are independent, so
\begin{align*}
\mu_{(X,Y)}(A\times C)=\mathbb P(X\in A,Y\in C)=\mathbb P(X\in A)\mathbb P(Y\in C).
\end{align*}
Because $f_X$ and $f_Y$ are densities,
\begin{align*}
\mathbb P(X\in A)\mathbb P(Y\in C)=\left(\int_A f_X(x)\,d\mathcal L^1(x)\right)\left(\int_C f_Y(y)\,d\mathcal L^1(y)\right).
\end{align*}
By Tonelli's theorem for non-negative functions on the product [measure space](/page/Measure%20Space) $(\mathbb R,\mathcal B(\mathbb R),\mathcal L^1)\times(\mathbb R,\mathcal B(\mathbb R),\mathcal L^1)$,
\begin{align*}
\left(\int_A f_X(x)\,d\mathcal L^1(x)\right)\left(\int_C f_Y(y)\,d\mathcal L^1(y)\right)=\int_{A\times C} f_X(x)f_Y(y)\,d\mathcal L^2(x,y)=\nu(A\times C).
\end{align*}
Thus $\mu_{(X,Y)}$ and $\nu$ agree on all measurable rectangles $A\times C$. These rectangles form a $\pi$-system generating $\mathcal B(\mathbb R^2)$, and both $\mu_{(X,Y)}$ and $\nu$ are probability measures. Hence the uniqueness theorem for measures generated by a $\pi$-system gives $\mu_{(X,Y)}=\nu$ on $\mathcal B(\mathbb R^2)$.
[/step]
[step:Rewrite the event $\{X+Y\in B\}$ as an integral over a sublevel region]
Let $B\in\mathcal B(\mathbb R)$ be arbitrary. Define the addition map
\begin{align*}
S:\mathbb R^2\to\mathbb R,\qquad S(x,y)=x+y.
\end{align*}
The map $S$ is continuous, hence Borel measurable, so $S^{-1}(B)\in\mathcal B(\mathbb R^2)$.
Since $X+Y=S\circ(X,Y)$, the event $\{X+Y\in B\}$ equals $\{(X,Y)\in S^{-1}(B)\}$. Using the joint density identified in the previous step, we obtain
\begin{align*}
\mathbb P(X+Y\in B)=\int_{S^{-1}(B)} f_X(x)f_Y(y)\,d\mathcal L^2(x,y).
\end{align*}
For any set $E$, let $\mathbb 1_E$ denote its indicator function, so $\mathbb 1_E(u)=1$ for $u\in E$ and $\mathbb 1_E(u)=0$ for $u\notin E$. Equivalently, using the indicator function $\mathbb 1_{S^{-1}(B)}:\mathbb R^2\to\{0,1\}$,
\begin{align*}
\mathbb P(X+Y\in B)=\int_{\mathbb R^2}\mathbb 1_B(x+y)f_X(x)f_Y(y)\,d\mathcal L^2(x,y).
\end{align*}
[/step]
[step:Use Tonelli and translation invariance to isolate the convolution kernel]
The function
\begin{align*}
h:\mathbb R^2\to[0,\infty),\qquad h(x,y)=\mathbb 1_B(x+y)f_X(x)f_Y(y)
\end{align*}
is Lebesgue measurable and non-negative. Therefore Tonelli's theorem applies and gives
\begin{align*}
\mathbb P(X+Y\in B)=\int_{\mathbb R} f_X(x)\left(\int_{\mathbb R}\mathbb 1_B(x+y)f_Y(y)\,d\mathcal L^1(y)\right)\,d\mathcal L^1(x).
\end{align*}
Fix $x\in\mathbb R$. Define the translation map
\begin{align*}
\tau_x:\mathbb R\to\mathbb R,\qquad \tau_x(y)=x+y.
\end{align*}
Translation invariance of $\mathcal L^1$ says that the pushforward measure $(\tau_x)_\#\mathcal L^1$ equals $\mathcal L^1$. Applying the change of variables $z=\tau_x(y)=x+y$, or equivalently $y=z-x$, gives
\begin{align*}
\int_{\mathbb R}\mathbb 1_B(x+y)f_Y(y)\,d\mathcal L^1(y)=\int_{\mathbb R}\mathbb 1_B(z)f_Y(z-x)\,d\mathcal L^1(z).
\end{align*}
Substituting this into the previous display gives
\begin{align*}
\mathbb P(X+Y\in B)=\int_{\mathbb R}\int_{\mathbb R}\mathbb 1_B(z)f_X(x)f_Y(z-x)\,d\mathcal L^1(z)\,d\mathcal L^1(x).
\end{align*}
The integrand $(x,z)\mapsto \mathbb 1_B(z)f_X(x)f_Y(z-x)$ is non-negative and Lebesgue measurable, so Tonelli's theorem applies again:
\begin{align*}
\mathbb P(X+Y\in B)=\int_{\mathbb R}\mathbb 1_B(z)\left(\int_{\mathbb R}f_X(x)f_Y(z-x)\,d\mathcal L^1(x)\right)\,d\mathcal L^1(z).
\end{align*}
Since $\mathbb 1_B(z)$ restricts the outer integral to $B$, this is
\begin{align*}
\mathbb P(X+Y\in B)=\int_B\left(\int_{\mathbb R}f_X(x)f_Y(z-x)\,d\mathcal L^1(x)\right)\,d\mathcal L^1(z).
\end{align*}
[guided]
We now perform the central computation. Throughout this computation, $\mathcal L^2$ denotes two-dimensional Lebesgue measure on $\mathbb R^2$, and for any set $E$, $\mathbb 1_E$ denotes the indicator function of $E$. The goal is to convert the probability of the event $\{X+Y\in B\}$ into an integral over the target variable $z$, because a density for $X+Y$ must satisfy a formula of the form
\begin{align*}
\mathbb P(X+Y\in B)=\int_B f_{X+Y}(z)\,d\mathcal L^1(z).
\end{align*}
From the previous step, we already know that
\begin{align*}
\mathbb P(X+Y\in B)=\int_{\mathbb R^2}\mathbb 1_B(x+y)f_X(x)f_Y(y)\,d\mathcal L^2(x,y).
\end{align*}
Define the measurable map
\begin{align*}
h:\mathbb R^2\to[0,\infty),\qquad h(x,y)=\mathbb 1_B(x+y)f_X(x)f_Y(y).
\end{align*}
This function is Lebesgue measurable and non-negative, so Tonelli's theorem applies without requiring prior integrability. Thus we may integrate first in $y$ and then in $x$:
\begin{align*}
\mathbb P(X+Y\in B)=\int_{\mathbb R} f_X(x)\left(\int_{\mathbb R}\mathbb 1_B(x+y)f_Y(y)\,d\mathcal L^1(y)\right)\,d\mathcal L^1(x).
\end{align*}
For fixed $x\in\mathbb R$, the inner integral is still written in terms of $y$, but the event depends on the sum $x+y$. To expose the variable of the sum, define
\begin{align*}
\tau_x:\mathbb R\to\mathbb R,\qquad \tau_x(y)=x+y.
\end{align*}
This translation map sends $y$ to $z=x+y$. Translation invariance of one-dimensional Lebesgue measure says that $(\tau_x)_\#\mathcal L^1=\mathcal L^1$. Therefore the substitution $z=x+y$, with inverse relation $y=z-x$, gives
\begin{align*}
\int_{\mathbb R}\mathbb 1_B(x+y)f_Y(y)\,d\mathcal L^1(y)=\int_{\mathbb R}\mathbb 1_B(z)f_Y(z-x)\,d\mathcal L^1(z).
\end{align*}
Substituting this identity into the iterated integral yields
\begin{align*}
\mathbb P(X+Y\in B)=\int_{\mathbb R}\int_{\mathbb R}\mathbb 1_B(z)f_X(x)f_Y(z-x)\,d\mathcal L^1(z)\,d\mathcal L^1(x).
\end{align*}
We want the outer integral to be over $z$, because $z$ is the value of the [random variable](/page/Random%20Variable) $X+Y$. The integrand
\begin{align*}
(x,z)\mapsto \mathbb 1_B(z)f_X(x)f_Y(z-x)
\end{align*}
is non-negative and Lebesgue measurable, so Tonelli's theorem applies a second time and permits the order of integration to be exchanged:
\begin{align*}
\mathbb P(X+Y\in B)=\int_{\mathbb R}\mathbb 1_B(z)\left(\int_{\mathbb R}f_X(x)f_Y(z-x)\,d\mathcal L^1(x)\right)\,d\mathcal L^1(z).
\end{align*}
Finally, multiplication by $\mathbb 1_B(z)$ restricts the outer integral to $B$, so
\begin{align*}
\mathbb P(X+Y\in B)=\int_B\left(\int_{\mathbb R}f_X(x)f_Y(z-x)\,d\mathcal L^1(x)\right)\,d\mathcal L^1(z).
\end{align*}
This is exactly the density identity, with the inner integral as the candidate density at the point $z$.
[/guided]
[/step]
[step:Choose a finite version of the convolution and conclude absolute continuity]
Define the measurable map $g:\mathbb R\to[0,\infty]$ by
\begin{align*}
g(z)=\int_{\mathbb R}f_X(x)f_Y(z-x)\,d\mathcal L^1(x).
\end{align*}
The measurability of $g$ follows from Tonelli's theorem applied to the non-negative [measurable function](/page/Measurable%20Function) $(x,z)\mapsto f_X(x)f_Y(z-x)$.
Taking $B=\mathbb R$ in the formula from the previous step gives
\begin{align*}
\int_{\mathbb R}g(z)\,d\mathcal L^1(z)=\mathbb P(X+Y\in\mathbb R)=1.
\end{align*}
Hence $g<\infty$ for $\mathcal L^1$-a.e. $z\in\mathbb R$. Let
\begin{align*}
N=\{z\in\mathbb R:g(z)=\infty\}.
\end{align*}
Then $\mathcal L^1(N)=0$. Define the finite-valued measurable map $f_{X+Y}:\mathbb R\to[0,\infty)$ by setting $f_{X+Y}(z)=g(z)$ for $z\in\mathbb R\setminus N$ and $f_{X+Y}(z)=0$ for $z\in N$.
For every $B\in\mathcal B(\mathbb R)$, changing the integrand on the $\mathcal L^1$-null set $N$ does not change the integral, so
\begin{align*}
\mathbb P(X+Y\in B)=\int_B g(z)\,d\mathcal L^1(z)=\int_B f_{X+Y}(z)\,d\mathcal L^1(z).
\end{align*}
Therefore the law of $X+Y$ is absolutely continuous with respect to $\mathcal L^1$, and $f_{X+Y}$ is a Lebesgue density for $X+Y$. Since $f_{X+Y}=g$ $\mathcal L^1$-a.e., this density is given, up to $\mathcal L^1$-a.e. equality, by
\begin{align*}
f_{X+Y}(z)=\int_{\mathbb R} f_X(x)f_Y(z-x)\,d\mathcal L^1(x).
\end{align*}
Thus $X+Y$ is a [continuous random variable](/page/Continuous%20Random%20Variable) with the asserted convolution density.
[/step]