I figured out a rigorous proof by myself.
If $t = 0$, then it is an integration of a real-valued function, and clearly, $\int_0^\infty e^{-x} dx = 1$.
If $t \neq 0$, without losing of generality, assume $t > 0$. Consider the contour below:
In the picture, $n$ is a positive number that will be sent to $\infty$, the top line passes the origin and the point $(-1, t)$. And we set $f(z) = e^z, z \in \mathbb{C}$. By Cauchy's integration theorem,
\begin{align}
0 = \int_\Gamma f(z) dz = \int_{\Gamma_1} e^z dz + \int_{\Gamma_2} e^z dz + \int_{\Gamma_3} e^z dz
\end{align}
Let's denote the angle between $\Gamma_1$ and the real axis by $\theta_0$.
Clearly, $\int_{\Gamma_3}e^z dz = \int_{-n}^0 e^x dx = 1 - e^{-n}$.
On $\Gamma_1$, $z$ has the representation $z = (it - 1)x, x \in (0, n/|1 - it|)$, thus
\begin{align*}
\int_{\Gamma_1}e^z dz = \int_0^{n/\sqrt{1 + t^2}}e^{(it - 1)x}(it - 1) dx =
(it - 1) \int_0^{n/\sqrt{1 + t^2}} e^{(it - 1)x} dx
\end{align*}
To get the desired result, it remains to show $\int_{\Gamma_2} e^z dz \to 0$ as $n \to \infty$. Let $z = ne^{i\theta}$ with $\theta \in (\theta_0, \pi)$. It follows that
\begin{align*}
& \left|\int_{\Gamma_2} e^z dz\right|
= \left|\int_{\theta_0}^\pi e^{ne^{i\theta}}nie^{i\theta} d\theta\right| \\
\leq & \int_{\theta_0}^\pi e^{n\cos\theta}n d\theta \\
\leq & ne^{n\cos{\theta_0}}(\pi - \theta_0) \to 0
\end{align*}
as $n \to \infty$. Here we used the fact that $\pi/2 < \theta_0 < \pi$.
Best Answer
Aside from some trigonometric substitutions and identities, we will need the following identity, which can be shown using integration by parts twice: $$ \int_0^{\infty}\cos(\alpha t)e^{-\lambda t}\,\mathrm{d}t=\frac{\lambda}{\alpha^2+\lambda^2}\tag{1} $$ We will also use the standard arctangent integral: $$ \int_0^\infty\frac{\mathrm{d}t}{a^2+t^2}=\frac\pi{2a}\tag{2} $$ Now $$ \begin{align} &\left(\int_0^\infty\color{#C00000}{\sin}(x^2) e^{-\lambda x^2}\,\mathrm{d}x\right)^2\\ &=\int_0^\infty\int_0^\infty \color{#C00000}{\sin}(x^2)\color{#C00000}{\sin}(y^2) e^{-\lambda(x^2+y^2)}\,\mathrm{d}y\,\mathrm{d}x\tag{3.1}\\ &=\frac12\int_0^\infty\int_0^\infty \left(\cos(x^2-y^2) \color{#FF0000}{-}\cos(x^2+y^2)\right) e^{-\lambda(x^2+y^2)}\,\mathrm{d}y\,\mathrm{d}x\tag{3.2}\\ &=\frac12\int_0^{\pi/2}\int_0^\infty \left(\cos(r^2\cos(2\phi)) \color{#FF0000}{-}\cos(r^2)\right)e^{-\lambda r^2} \,r\,\mathrm{d}r\,\mathrm{d}\phi\tag{3.3}\\ &=\frac14\int_0^{\pi/2}\int_0^\infty \left(\cos(s\cos(2\phi)) \color{#FF0000}{-}\cos(s)\right) e^{-\lambda s} \,\mathrm{d}s\,\mathrm{d}\phi\tag{3.4}\\ &=\frac14\int_0^{\pi/2}\left(\frac{\lambda}{\cos^2(2\phi)+\lambda^2} \color{#FF0000}{-}\frac{\lambda}{1+\lambda^2}\right)\,\mathrm{d}\phi\tag{3.5}\\ &=\frac12\int_0^{\pi/4} \frac{\lambda}{\cos^2(2\phi)+\lambda^2}\,\mathrm{d}\phi \color{#FF0000}{-}\frac{\lambda\pi/8}{1+\lambda^2}\tag{3.6}\\ &=\frac14\int_0^{\pi/4} \frac{\lambda\,\mathrm{d}\tan(2\phi)} {1+\lambda^2+\lambda^2\tan^2(2\phi)} \color{#FF0000}{-}\frac{\lambda\pi/8}{1+\lambda^2}\tag{3.7}\\ &=\frac14\int_0^\infty\frac{\mathrm{d}t}{1+\lambda^2+t^2} \color{#FF0000}{-}\frac{\lambda\pi/8}{1+\lambda^2}\tag{3.8}\\ &=\frac{\pi/8}{\sqrt{1+\lambda^2}} \color{#FF0000}{-}\frac{\lambda\pi/8}{1+\lambda^2}\tag{3.9} \end{align} $$
$(3.1)$ change the square of the integral into a double integral
$(3.2)$ use $2\color{#C00000}{\sin}(A)\color{#C00000}{\sin}(B)=\cos(A-B)\color{#FF0000}{-}\cos(A+B)$
$(3.3)$ convert to polar coordinates
$(3.4)$ substitute $s=r^2$
$(3.5)$ apply $(1)$
$(3.6)$ pull out the constant and apply symmetry to reduce the domain of integration
$(3.7)$ multiply numerator and denominator by $\sec^2(2\phi)$
$(3.8)$ substitute $t=\lambda\tan(2\phi)$
$(3.9)$ apply $(2)$
Finally, take the square root of both sides of $(3)$ and let $\lambda\to0^+$ to get $$ \int_0^\infty\color{#C00000}{\sin}(x^2)\,\mathrm{d}x=\sqrt{\frac\pi8}\tag{4} $$
Addendum
I just noticed that the same proof works for $$ \int_0^\infty\cos(x^2)\,\mathrm{d}x=\sqrt{\frac\pi8}\tag{5} $$ if each red $\color{#C00000}{\sin}$ is changed to $\cos$ and each red $\color{#FF0000}{-}$ sign is changed to $+$.