[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*}[/step]

[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: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]