Don't quite understand an answer for finding a pdf from a joint pdf Here's the problem

I have the solution

However, I don't understand their choice of integration bounds. I understand that the function g doesn't exist at zero, but I miss why they choose the bounds of [y,1]. I can see that it works, but I don't understand the why.
 A: As pointed out in the comments, the answer to your problem is already explicitly included in the solution you provided. That said, this answer is meant to elaborate on why the answer is what it is.
In these sorts of problems, where you are performing transformations on multiple ranfom variables, it is of great help to draw pictures for the support of the density function before and after the transformation. In this specific problem recall that we are performing the transformation $(Y,Z)=(X_1,X_1X_2)$, which gives the inverse transformation $(X_1,X_2)=(Z,Y/Z)$. Now, the original density $f_{X_1X_2}(x_1,x_2)$ has support on the region $0<x_1,x_2<1$. Taking the four corners of this region and applying are transformation to them gives us a way to draw the support of $g_{YZ}(y,z)$. We have,
$$
\left(
\begin{array}{c}
 (x_1,x_2) \to (y,z)\\
 (0,0) \to (0,0) \\
 (1,0) \to (0,1) \\
 (0,1) \to (0,0) \\
 (1,1) \to (1,1) \\
\end{array}
\right),
$$
which allows us to draw the following diagram.

Notice that for any fixed $y\in[0,1]$ the support of $g$ with respect to $z$ is bounded below by $y$ and bounded above by $1$. Hence, to get the density of $Y$ we integrate joint density $g_{YZ}(y,z)$ over $z$ to get
$$
f_Y(y)=\int_y^1g_{YZ}(y,z)\,\mathrm dz,
$$
which gives the result in the solution you posted.
