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 $+$.
$\newcommand{\+}{^{\dagger}}
\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle}
\newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace}
\newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack}
\newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,}
\newcommand{\dd}{{\rm d}}
\newcommand{\down}{\downarrow}
\newcommand{\ds}[1]{\displaystyle{#1}}
\newcommand{\expo}[1]{\,{\rm e}^{#1}\,}
\newcommand{\fermi}{\,{\rm f}}
\newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,}
\newcommand{\half}{{1 \over 2}}
\newcommand{\ic}{{\rm i}}
\newcommand{\iff}{\Longleftrightarrow}
\newcommand{\imp}{\Longrightarrow}
\newcommand{\isdiv}{\,\left.\right\vert\,}
\newcommand{\ket}[1]{\left\vert #1\right\rangle}
\newcommand{\ol}[1]{\overline{#1}}
\newcommand{\pars}[1]{\left(\, #1 \,\right)}
\newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\newcommand{\pp}{{\cal P}}
\newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,}
\newcommand{\sech}{\,{\rm sech}}
\newcommand{\sgn}{\,{\rm sgn}}
\newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}}
\newcommand{\ul}[1]{\underline{#1}}
\newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}
\newcommand{\wt}[1]{\widetilde{#1}}$
\begin{align}
\color{#c00000}{\int_{0}^{1}{\dd x \over \pars{x^{2} - x^{3}}^{1/3}}}&
=\int_{\infty}^{1}{-\dd x/x^{2} \over \pars{1/x^{2} - 1/x^{3}}^{1/3}}
=\int_{1}^{\infty}{\dd x \over x\pars{x - 1}^{1/3}}
=\color{#c00000}{\int_{0}^{\infty}{x^{-1/3}\,\dd x \over x + 1}}
\end{align}
Use the following contour with
$\ds{z^{-1/3} = \verts{z}^{-1/3}\expo{-\ic\phi\pars{z}/3}\,,\qquad
0 < \phi\pars{z} < 2\pi}$:
\begin{align}
\color{#c00000}{\int_{0}^{\infty}{x^{-1/3}\,\dd x \over x + 1}}
&=2\pi\ic\expo{-\pi\ic/3}
-\int_{\infty}^{0}{x^{-1/3}\expo{-2\pi\ic/3}\,\dd x \over x + 1}
=2\pi\ic\expo{-\pi\ic/3}
+\expo{-2\pi\ic/3}\color{#c00000}{\int_{0}^{\infty}{x^{-1/3}\,\dd x \over x + 1}}
\end{align}
\begin{align}
\color{#c00000}{\int_{0}^{\infty}{x^{-1/3}\,\dd x \over x + 1}}&=
2\pi\ic\,{\expo{-\pi\ic/3} \over 1 - \expo{-2\pi\ic/3}}
=
2\pi\ic\,{1 \over \expo{\pi\ic/3} - \expo{-\pi\ic/3}}
={2\pi\ic \over 2\ic\sin\pars{\pi/3}}={\pi \over \sin\pars{\pi/3}}
\\[3mm]&={\pi \over \root{3}/2} = {2\root{3} \over 3}\,\pi
\end{align}
$$
\color{#66f}{\large\int_{0}^{1}{\dd x \over \pars{x^{2} - x^{3}}^{1/3}}
={2\root{3} \over 3}\,\pi} \approx 3.6276
$$
Best Answer
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$.