If $f$ is a probability density of a random variable $X$, then we have
$$
\Pr(X\in A) = \int\limits_A f(x)\, dx.
$$
Now suppose $X$ and $Y$ are independent random variables with respective densities $f$ and $g.$ Then
\begin{align}
\Pr(X\in A\ \&\ Y\in B) & = \Pr(X\in A)\Pr(Y\in B) \\[10pt]
& = \int\limits_A f(x)\,dx \cdot \int\limits_B g(x)\, dx \\[10pt]
& = \int\limits_A f(x)\,dx \cdot \int\limits_B g(y)\, dy
\end{align}
because a bound variable can be freely re-named as long as the new name is not already taken. And you in an integral $\int_A \cdot\cdots \, dx$ with respect to $x,$, anything not depending on $x$ is a "constant", so you can write $\int_A f(x)\,dx \cdot 5 = \int_A f(x)\cdot5\, dx,$ etc., and accordingly we have
\begin{align}
& = \int\limits_A f(x) \left( \int\limits_B g(y)\, dy \right) \, dx \\[10pt]
\text{and then } & = \int\limits_A \left( \int\limits_B f(x)g(y) \, dy \right) \, dx \tag 1 \\[10pt]
& = \iint\limits_{A\times B} f(x) g(y) \, d(x,y) \text{ by Tonelli's theorem}. \tag 2
\end{align}
Now, what does this mean? The "double integral"
$$
\iint\limits_{A\times B} \cdots \, d(x,y)
$$
is an integral with respect to the $2$-dimensional measure in the plane, in which the measure of a region is its area. The "iterated integral"
$$
\int\limits_A \left( \int\limits_B \cdots\, dy \right) \, dx
$$
involves two integrals with respect to $1$-dimensional measure in the line, in which the measure of an integral is its length. Tonelli's theorem says that as long as the function getting integrated is everywhere non-negative, then the values of the iterated integral in line $(1)$ and the double integral in line $(2)$ are the same (regardless of whether the value is a finite number or $+\infty$). (There is also the related Fubini's theorem, which says that if the double integral of the absolute value of the function is finite, then iterated integral and the double integral are equal.)
So now we have
$$
\Pr(X\in A\ \&\ Y\in B) = \iint\limits_{A\times B} f(x) g(y) \, d(x,y).
$$
Therefore the function $(x,y) \mapsto f(x) g(y)$ behaves like the density function of the random pair $(X,Y),$ at least as far as product sets like $A\times B$ are concerned. Does it still behave like the density of $(X,Y)$ if instead we integrate over, for example a disk in the plane? Here one would need to show that if it works for things like $A\times B,$ then it also works for things like a disk and for all other "measurable sets" in the plane. That's a fairly hairy argument, but if you've got that, then we're done.