Where do the Fresnel integrals come from?
A monochromatic scalar wave of wavelength $\lambda$ has wavenumber $k=\frac{2\pi}{\lambda},$ and propagation over a distance $r$ contributes the phase $e^{ikr}.$ Let $\xi$ be a coordinate in an aperture plane and $x$ a coordinate in an observation plane a distance $z$ away. Their separation is
\[r=\sqrt{z^2+(x-\xi)^2}.\]When $|x-\xi|\ll z$, the binomial expansion gives the Fresnel approximation
\[r =z\sqrt{1+\frac{(x-\xi)^2}{z^2}} \simeq z+\frac{(x-\xi)^2}{2z}.\]Therefore,
\[e^{ikr} \simeq e^{ikz} \exp\left[\frac{ik}{2z}(x-\xi)^2\right].\]Apart from an overall prefactor, the Huygens–Fresnel diffraction formula is then
\[U(x,z)\propto \int_{\text{aperture}} U_0(\xi) \exp\left[\frac{ik}{2z}(x-\xi)^2\right]\,\mathrm{d}\xi.\]For uniform illumination, observation at $x=0$, and an aperture extending from $\xi=0$ to $\xi=L$,
\[U(0,z)\propto \int_0^L \exp\left(\frac{i\pi\xi^2}{\lambda z}\right)\,\mathrm{d}\xi.\]Introduce the dimensionless variables
\[t=\sqrt{\frac{2}{\lambda z}}\,\xi, \qquad u=\sqrt{\frac{2}{\lambda z}}\,L.\]The diffraction field becomes
\[U(0,z)\propto \sqrt{\frac{\lambda z}{2}} \int_0^u e^{i\pi t^2/2}\,\mathrm{d}t.\]Using Euler’s formula,
\[e^{i\pi t^2/2} =\cos\left(\frac{\pi t^2}{2}\right) +i\sin\left(\frac{\pi t^2}{2}\right),\]we obtain
\[U(0,z)\propto \sqrt{\frac{\lambda z}{2}} \left[\mathcal{C}(u)+i\mathcal{S}(u)\right], \qquad \mathcal{C}(u) =\int_0^u\cos\left(\frac{\pi t^2}{2}\right)\,\mathrm{d}t, \qquad \mathcal{S}(u) =\int_0^u\sin\left(\frac{\pi t^2}{2}\right)\,\mathrm{d}t.\]These are the Fresnel integrals. Their parametric plot
\[u\longmapsto\bigl(\mathcal{C}(u),\mathcal{S}(u)\bigr)\]is called the Cornu spiral.
After rescaling once more, it is enough to evaluate
\[F=\int_0^\infty e^{ix^2}\,\mathrm{d}x\]and hence its real and imaginary parts,
\[C=\int_0^\infty\cos(x^2)\,\mathrm{d}x, \qquad S=\int_0^\infty\sin(x^2)\,\mathrm{d}x.\]These integrals are not absolutely convergent, because $|e^{ix^2}|=1$. They nevertheless converge as improper oscillatory integrals. Indeed, integration by parts gives, for $B>A>0$,
\[\int_A^B e^{ix^2}\,\mathrm{d}x =\int_A^B\frac{1}{2ix}\,\mathrm{d}\left(e^{ix^2}\right) =\left[\frac{e^{ix^2}}{2ix}\right]_A^B +\frac{1}{2i}\int_A^B\frac{e^{ix^2}}{x^2}\,\mathrm{d}x.\]The right-hand side tends uniformly to zero as $A\to\infty$, so the Cauchy criterion proves convergence.
1. Schwinger parameterization
Start with the definition of the Gamma function:
\[\Gamma(\nu) =\int_0^\infty t^{\nu-1}e^{-t}\,\mathrm{d}t, \qquad \operatorname{Re}\nu>0.\]For the moment, let $A>0$ and make the substitution $t=As$ and $\mathrm{d}t=A\,\mathrm{d}s$. Then
\[\Gamma(\nu) =\int_0^\infty (As)^{\nu-1}e^{-As}A\,\mathrm{d}s =A^\nu\int_0^\infty s^{\nu-1}e^{-As}\,\mathrm{d}s.\]Dividing by $A^\nu\Gamma(\nu)$ gives the Schwinger parameterization
\[\frac{1}{A^\nu} =\frac{1}{\Gamma(\nu)} \int_0^\infty s^{\nu-1}e^{-As}\,\mathrm{d}s, \qquad \operatorname{Re}A>0,\quad \operatorname{Re}\nu>0.\]The derivation above used real $A>0$. With the principal branch of $A^\nu$, both sides are analytic in the right half-plane, so the identity extends to every complex $A$ with $\operatorname{Re}A>0$.
For the Fresnel integral, a small positive regulator is essential. Define
\[F_\varepsilon =\int_0^\infty e^{-(\varepsilon-i)x^2}\,\mathrm{d}x, \qquad \varepsilon>0.\]Unlike the original integral, $F_\varepsilon$ is absolutely convergent. Set $s=x^2$, so that $x=\sqrt{s}$ and $\mathrm{d}x=\frac12s^{-1/2}\,\mathrm{d}s$. Then
\[F_\varepsilon =\frac12\int_0^\infty s^{-1/2}e^{-(\varepsilon-i)s}\,\mathrm{d}s.\]Apply the Schwinger identity with
\[\nu=\frac12, \qquad A=\varepsilon-i.\]Since $\operatorname{Re}(\varepsilon-i)=\varepsilon>0$ and $\Gamma(\frac12)=\sqrt{\pi}$, we obtain
\[F_\varepsilon=\frac12\Gamma\left(\frac12\right)(\varepsilon-i)^{-1/2}=\frac{\sqrt{\pi}}{2}(\varepsilon-i)^{-1/2}.\]We use the principal square root, with argument in $(-\pi,\pi)$. Writing
\[\varepsilon-i =\sqrt{1+\varepsilon^2}\, e^{-i\arctan(1/\varepsilon)} \quad \Rightarrow \quad (\varepsilon-i)^{-1/2} =(1+\varepsilon^2)^{-1/4} e^{\frac{i}{2}\arctan(1/\varepsilon)}.\]As $\varepsilon\to0^+$,
\[\lim_{\varepsilon\to0^+}F_\varepsilon =\frac{\sqrt{\pi}}{2}e^{i\pi/4} =\frac{\sqrt{\pi}}{2\sqrt{2}}(1+i).\]It remains to connect this Abel-regularized limit to the original improper integral. For a fixed $R$,
\[\lim_{\varepsilon\to0^+} \int_0^R e^{-(\varepsilon-i)x^2}\,\mathrm{d}x =\int_0^R e^{ix^2}\,\mathrm{d}x\]by dominated convergence. The tails can be controlled uniformly for small $\varepsilon$ by integration by parts:
\[\int_R^\infty e^{-(\varepsilon-i)x^2}\,\mathrm{d}x =\int_R^\infty \frac{-1}{2(\varepsilon-i)x}\, \mathrm{d}\left(e^{-(\varepsilon-i)x^2}\right) =\frac{e^{-(\varepsilon-i)R^2}}{2(\varepsilon-i)R} -\frac{1}{2(\varepsilon-i)} \int_R^\infty\frac{e^{-(\varepsilon-i)x^2}}{x^2}\,\mathrm{d}x.\]For $0<\varepsilon\leq1$, its absolute value is at most $1/R$. Thus we may first choose $R$ large and then let $\varepsilon\to0^+$. This proves that the regulated limit equals the ordinary Fresnel integral:
\[F=\frac{\sqrt{\pi}}{2}e^{i\pi/4}.\]Taking real and imaginary parts yields
\[\boxed{ \int_0^\infty\cos(x^2)\,\mathrm{d}x =\int_0^\infty\sin(x^2)\,\mathrm{d}x =\sqrt{\frac{\pi}{8}} }\]2. Complex-analysis method
The same result follows geometrically by rotating the contour through an angle of $\pi/4$. Consider the entire function
\[f(z)=e^{iz^2}\]on the positively oriented boundary of the sector with vertices $0$, $R$, and $Re^{i\pi/4}$. Denote the circular arc by
\[\gamma_R(\theta)=Re^{i\theta}, \qquad 0\leq\theta\leq\frac{\pi}{4}.\]
Since $f$ is entire, Cauchy’s integral theorem gives
\[\int_0^R e^{ix^2}\,\mathrm{d}x +\int_{\gamma_R}e^{iz^2}\,\mathrm{d}z -\int_0^R e^{i(re^{i\pi/4})^2}e^{i\pi/4}\,\mathrm{d}r =0.\]On the diagonal ray, $(re^{i\pi/4})^2=ir^2$, so
\[e^{i(re^{i\pi/4})^2}=e^{-r^2}.\]Consequently,
\[\int_0^R e^{ix^2}\,\mathrm{d}x =e^{i\pi/4}\int_0^R e^{-r^2}\,\mathrm{d}r -\int_{\gamma_R}e^{iz^2}\,\mathrm{d}z.\]We now show that the arc contribution vanishes. On $z=Re^{i\theta}$,
\[\left|e^{iz^2}\right| =e^{-R^2\sin(2\theta)}.\]Because $\sin(2\theta)\geq4\theta/\pi$ for $0\leq\theta\leq\pi/4$,
\[\left|\int_{\gamma_R}e^{iz^2}\,\mathrm{d}z\right| \leq R\int_0^{\pi/4}e^{-R^2\sin(2\theta)}\,\mathrm{d}\theta \leq R\int_0^{\pi/4}e^{-4R^2\theta/\pi}\,\mathrm{d}\theta =\frac{\pi}{4R}\left(1-e^{-R^2}\right) \longrightarrow0.\]Letting $R\to\infty$ and using the Gaussian integral
\[\int_0^\infty e^{-r^2}\,\mathrm{d}r=\frac{\sqrt{\pi}}{2}\]therefore gives
\[F=e^{i\pi/4}\frac{\sqrt{\pi}}{2},\]in agreement with the Schwinger-parameter calculation.
The contour rotation also explains the phase $e^{i\pi/4}$. Along the real axis, $e^{iz^2}$ oscillates without decaying. Rotating $z=re^{i\theta}$ changes its magnitude to $e^{-r^2\sin(2\theta)}$; at $\theta=\pi/4$, the oscillatory exponential becomes the decaying Gaussian $e^{-r^2}$, while the line element contributes $\mathrm{d}z=e^{i\pi/4}\,\mathrm{d}r$.
3. Common normalization and a scaled formula
Many texts define the Fresnel functions by
\[\mathcal{C}(u)=\int_0^u\cos\left(\frac{\pi t^2}{2}\right)\,\mathrm{d}t, \qquad \mathcal{S}(u)=\int_0^u\sin\left(\frac{\pi t^2}{2}\right)\,\mathrm{d}t.\]Using $x=\sqrt{\pi/2}\,t$ in the result above gives
\[\mathcal{C}(\infty)=\mathcal{S}(\infty)=\frac12.\]More generally, for $a>0$, the substitution $u=\sqrt{a}\,x$ gives
\[\int_0^\infty e^{iax^2}\,\mathrm{d}x =\frac{\sqrt{\pi}}{2\sqrt{a}}e^{i\pi/4},\]and hence
\[\boxed{ \int_0^\infty\cos(ax^2)\,\mathrm{d}x =\int_0^\infty\sin(ax^2)\,\mathrm{d}x =\sqrt{\frac{\pi}{8a}} }\]