Generalizing the trick for integrating $\int_{-\infty}^\infty e^{-x^2}\mathrm dx$?

I would like to answer to your question about Fresnel integral as I did this in my undergraduate studies. So, let us consider the integrals

$$S_1=\int_{-\infty}^\infty dx\sin(x^2) \qquad C_1=\int_{-\infty}^\infty dx\cos(x^2).$$

We want to apply the same technique used for Gauss integral in this case and consider the two-dimensional integrals

$$S_2=\int_{-\infty}^\infty\int_{-\infty}^\infty dxdy\sin(x^2+y^2) \qquad C_2=\int_{-\infty}^\infty\int_{-\infty}^\infty dxdy\cos(x^2+y^2).$$

If you go to polar coordinates, these integrals doe not converge. So, we introduce a convergence factor in the following way

$$S_2(\epsilon)=\int_{-\infty}^\infty\int_{-\infty}^\infty dxdy e^{-\epsilon(x^2+y^2)}\sin(x^2+y^2)$$ $$C_2(\epsilon)=\int_{-\infty}^\infty\int_{-\infty}^\infty dxdy e^{-\epsilon(x^2+y^2)}\cos(x^2+y^2).$$

Then one has, moving to polar coordinates,

$$S_2(\epsilon)=2\pi\int_0^\infty \rho d\rho e^{-\epsilon\rho^2}\sin(\rho^2) \qquad C_2(\epsilon)=2\pi\int_0^\infty \rho d\rho e^{-\epsilon\rho^2}\cos(\rho^2)$$

that is

$$S_2(\epsilon)=\pi\int_0^\infty dx e^{-\epsilon x}\sin(x) \qquad C_2(\epsilon)=\pi\int_0^\infty dx e^{-\epsilon x}\cos(x).$$

These integrals are well known and give

$$S_2(\epsilon)=\frac{\pi}{1+\epsilon^2} \qquad C_2(\epsilon)=\frac{\epsilon\pi}{1+\epsilon^2}$$

noting that integration variables are dummy. Now, in this case one can take the limit for $\epsilon\rightarrow 0$ producing

$$S_2=\int_{-\infty}^\infty\int_{-\infty}^\infty dxdy\sin(x^2+y^2)=\pi \qquad C_2=\int_{-\infty}^\infty\int_{-\infty}^\infty dxdy\cos(x^2+y^2)=0$$

By applying simple trigonometric formulas you will get back the value of the Fresnel integrals. But now, as a bonus, you have got the value of these integrals in two dimensions.

Robert Dawson wrote about the limitations of this method in this article in the American Mathematical Monthly.

The distance metric that we use is quadratic in nature. No matter what dimensional $\mathbb{R}^n$ we are working in, the distance between two points is given by $r^2=\sum(\Delta x_i)^2$. It's the ability to swap out $\sum(\Delta x_i)^2$ with $r^2$ that helps compute $\int_{-\infty}^{\infty}e^{-x^2}\,dx$. No higher dimensional space will make this particular exchange any nicer. So this helps explain why the property you mention ($f(x^2)f(y^2)=f(x^2+y^2)$) is important.

On the other hand, what happens here is an exchange between Cartesian coordinates for $\mathbb{R}^2$ and polar coordinates. So could we do something similar in $\mathbb{R}^3$ with spherical coordinates?

$$\left(\int_{\mathbb{R}}f(x)\,dx\right)^3=\iiint_{\mathbb{R}^3} f(x)f(y)f(z)\,dxdydz=\iiint_{\mathbb{R}^3} F(\rho,\theta,\varphi)\,d\rho d\theta d\varphi$$

If $F(\rho,\theta,\varphi)$ breaks up as $g(\rho)h(\theta)k(\varphi)$, then you might have a slick way to compute $\int_{\mathbb{R}}f(x)\,dx$.