Looking at other question, I see several with $P(X^2+Y^2<1)$ but none with $P(X^2+Y^2<t)$.
One question is how to use one of these other solutions with t instead of 1.
A fellow classmate suggested letting $Z=X^2+Y^2$ then we can use $Z$~ χ$^2$ distribution with 2 degrees of freedom. I am not sure how/why I can use this. His method looks a lot shorter than the ones posted for the other solutions. But I want to understand more on how and why.
Let $Z$ be a standard normal random variable and define $W = Z^2$. Then for $w \ge 0$ we have $$F_W(w) = \Pr[Z^2 \le w] = \Pr[-\sqrt{w} \le Z \le \sqrt{w}] = 2 \Pr[Z \le \sqrt{w}] = 2 F_Z(\sqrt{w}).$$ It follows by the chain rule that $$f_W(w) = 2 \cdot \frac{1}{2\sqrt{w}} f_Z(\sqrt{w}) = \frac{f_Z(\sqrt{w})}{\sqrt{w}} = \frac{e^{-w/2}}{\sqrt{2\pi w}}, \quad w > 0.$$ The moment-generating function of $W$ is therefore $$M_W(t) = \operatorname{E}[e^{tW}] = \int_{w = 0}^\infty e^{tw} \frac{e^{-w/2}}{\sqrt{2\pi w}} \, dw = \frac{1}{\sqrt{2\pi}} \int_{w=0}^\infty \frac{e^{-(1/2 - t)w}}{\sqrt{w}} \, dw = \frac{1}{\sqrt{\pi(1-2t)}} \int_{u=0}^\infty \frac{e^{-u}}{\sqrt{u}} \, du.$$ We can now perform a substitution $u = v^2$, $du = 2v \, dv$, to obtain $$M_W(t) = \frac{1}{\sqrt{\pi(1-2t)}} \int_{v=0}^\infty 2e^{-v^2} \, dv = \frac{1}{\sqrt{1-2t}}.$$ We could have also obtained this result by recognizing that $$W \sim \operatorname{Gamma}(\alpha = 1/2, \theta = 2)$$ where $\alpha$ is the shape and $\theta$ is the scale parameter, and since the MGF of a gamma distribution is $M_W(t) = (1-\theta t)^{-\alpha}$, we obtain the result above. Finally, suppose $Z_1, \ldots, Z_n$ are iid standard normal variables. Then $$X_n = \sum_{i=1}^n Z_i^2$$ will have an MGF that equals the product of the individual MGFs of $Z_i^2$; that is to say, $$M_{X_n}(t) = \left(M_W(t)\right)^n = (1 - 2t)^{-n/2}.$$ But as we have seen, this simply tells us $$X_n \sim \operatorname{Gamma}(\alpha = n/2, \theta = 2),$$ with PDF $$f_{X_n}(x) = \frac{x^{n-1/2} e^{-x/2}}{2^{n/2} \Gamma(n/2)}, \quad x > 0.$$ This is precisely the $\chi^2$ density with $n$ degrees of freedom. When $n = 2$ as in your case, we easily obtain $$f_{X_2}(x) = \frac{1}{2} e^{-x/2}, \quad x > 0;$$ i.e., the sum of the squares of two iid standard normal random variables is exponential with mean $2$.