Step 1. Reduction of the integral
Let $I$ denote the integral in question. With the change of variable $v = \frac{\pi x^2}{2}$, we have
$$ I = \frac{1}{\pi} \int_{0}^{\infty} \left\{ (1 - 2 C(x) )^{2} + (1 - 2S(x))^{2} \right\}^{2} \, dv $$
where $x = \sqrt{2v / \pi}$ is understood as a function of $v$. By noting that
$$ 1-2 S(x) = \sqrt{\frac{2}{\pi}} \int_{v}^{\infty} \frac{\sin u}{\sqrt{u}} \, du \quad \text{and} \quad 1-2 C(x) = \sqrt{\frac{2}{\pi}} \int_{v}^{\infty} \frac{\cos u}{\sqrt{u}} \, du, $$
we can write $I$ as
$$ I = \frac{4}{\pi^{3}} \int_{0}^{\infty} \left| A(v) \right|^{4} \, dv \tag{1} $$
where $A(v)$ denotes the function defined by
$$ A(v) = \int_{v}^{\infty} \frac{e^{iu}}{\sqrt{u}} \, du. $$
Step 2. Simplification of $\left| A(v) \right|^2$.
Now we want to simplify $\left| A(v) \right|^2$. To this end, we note that for $\Re u > 0$,
$$ \frac{1}{\sqrt{u}}
= \frac{1}{\Gamma\left(\frac{1}{2}\right)} \frac{\Gamma\left(\frac{1}{2}\right)}{u^{1/2}}
= \frac{1}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-ux}}{\sqrt{x}} \, dx
= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} e^{-ux^{2}} \, dx \tag{2} $$
Using this identity,
\begin{align*}
A(v)
&= \frac{2}{\sqrt{\pi}} \int_{v}^{\infty} e^{iu} \int_{0}^{\infty} e^{-u x^2} \, dx du
= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \int_{v}^{\infty} e^{-(x^2-i)u} \, du dx \\
&= \frac{2 e^{iv}}{\sqrt{\pi}} \int_{0}^{\infty} e^{-v x^2} \int_{0}^{\infty} e^{-(x^2-i)u} \, du dx
= \frac{2 e^{iv}}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-v x^2}}{x^2-i} \, dx.
\end{align*}
Thus by the polar coordinate change $(x, y) \mapsto (r, \theta)$ followed by the substitutions $r^2 = s$ and $\tan \theta = t$, we obtain
\begin{align*}
\left| A(v) \right|^2
&= A(v) \overline{A(v)}
= \frac{4}{\pi} \int_{0}^{\infty} \int_{0}^{\infty} \frac{e^{-v (x^2+y^2)}}{(x^2-i)(y^2 + i)} \, dxdy \\
&= \frac{4}{\pi} \int_{0}^{\infty} \int_{0}^{\frac{\pi}{2}} \frac{r e^{-v r^2}}{(r^2 \cos^{2}\theta-i)(r^2 \sin^{2}\theta + i)} \, d\theta dr \\
&= \frac{2}{\pi} \int_{0}^{\infty} \int_{0}^{\frac{\pi}{2}} \frac{e^{-v s}}{(s \cos^{2}\theta-i)(s \sin^{2}\theta + i)} \, d\theta ds \\
&= \frac{2}{\pi} \int_{0}^{\infty} \frac{e^{-v s}}{s} \int_{0}^{\frac{\pi}{2}} \left( \frac{1}{s \cos^{2}\theta-i} + \frac{1}{s \sin^{2}\theta + i} \right) \, d\theta ds \\
&= \frac{2}{\pi} \int_{0}^{\infty} \frac{e^{-v s}}{s} \int_{0}^{\infty} \left( \frac{1}{s -i(t^2 + 1)} + \frac{1}{s t^2 + i (t^2 + 1)} \right) \, dt ds.
\end{align*}
Evaluation of the inner integral is easy, and we obtain
\begin{align*}
\left| A(v) \right|^2
&= 2 \int_{0}^{\infty} \frac{e^{-v s}}{s} \Re \left( \frac{i}{\sqrt{1 + is}} \right) \, ds.
\end{align*}
Applying $(2)$ again, we find that
\begin{align*}
\left| A(v) \right|^2
&= 2 \int_{0}^{\infty} \frac{e^{-v s}}{s} \Re \left( \frac{i}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-(1+is)u}}{\sqrt{u}} \, du \right) \, ds \\
&= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-v s}}{s} \int_{0}^{\infty} \frac{e^{-u} \sin (su)}{\sqrt{u}} \, du\, ds \\
&= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-u}}{\sqrt{u}} \int_{0}^{\infty} \frac{\sin (su)}{s} \, e^{-v s} \, ds\, du \\
&= \frac{2}{\sqrt{\pi}} \int_{0}^{\infty} \frac{e^{-u}}{\sqrt{u}} \arctan \left( \frac{u}{v} \right) \, du \\
&= \frac{4\sqrt{v}}{\sqrt{\pi}} \int_{0}^{\infty} e^{-vx^{2}} \arctan (x^2) \, dx \qquad (u = vx^2) \tag{3}
\end{align*}
Here, we exploited the identity
$$ \int_{0}^{\infty} \frac{\sin x}{x} e^{-sx} \, dx = \arctan \left(\frac{1}{s}\right), $$
which can be proved by differentiating both sides with respect to $s$.
Step 3. Evaluation of $I$.
Plugging $(3)$ to $(1)$ and applying the polar coordinate change, $I$ reduces to
\begin{align*}
I
&= \frac{64}{\pi^{4}} \int_{0}^{\infty} \int_{0}^{\infty} \int_{0}^{\infty} v e^{-v(x^{2}+y^{2})} \arctan (x^2) \arctan (x^2) \, dx dy dv \\
&= \frac{64}{\pi^{4}} \int_{0}^{\infty} \int_{0}^{\infty} \frac{\arctan (x^2) \arctan (x^2)}{(x^2 + y^2)^2} \, dx dy \\
&= \frac{64}{\pi^{4}} \int_{0}^{\frac{\pi}{2}} \int_{0}^{\infty} \frac{\arctan (r^2 \cos^2 \theta) \arctan (r^2 \sin^2 \theta)}{r^3} \, dr d\theta \\
&= \frac{32}{\pi^{4}} \int_{0}^{\frac{\pi}{2}} \int_{0}^{\infty} \frac{\arctan (s \cos^2 \theta) \arctan (s \sin^2 \theta)}{s^2} \, ds d\theta. \qquad (s = r^2) \tag{4}
\end{align*}
Now let us denote
$$ J(u, v) = \int_{0}^{\infty} \frac{\arctan (us) \arctan (vs)}{s^2} \, ds. $$
Then a simple calculation shows that
$$ \frac{\partial^{2} J}{\partial u \partial v} J(u, v) = \int_{0}^{\infty} \frac{ds}{(u^2 s^2 + 1)(v^2 s^2 + 1)} = \frac{\pi}{2(u+v)}. $$
Indeed, both the contour integration method or the partial fraction decomposition method work here. Integrating, we have
$$ J(u, v) = \frac{\pi}{2} \left\{ (u+v) \log(u+v) - u \log u - v \log v \right\}. $$
Plugging this to $(4)$, it follows that
\begin{align*}
I
&= -\frac{64}{\pi^{3}} \int_{0}^{\frac{\pi}{2}} \sin^2 \theta \log \sin \theta \, d\theta
= -\frac{16}{\pi^{3}} \frac{\partial \beta}{\partial z}\left( \frac{3}{2}, \frac{1}{2} \right)
\end{align*}
where $\beta(z, w)$ is the beta function, satisfying the following beta function identity
$$ \beta (z, w) = 2 \int_{0}^{\infty} \sin^{2z-1}\theta \cos^{2w-1} \theta \, d\theta = \frac{\Gamma(z)\Gamma(w)}{\Gamma(z+w)}. $$
Therefore we have
\begin{align*}
I
&= \frac{16}{\pi^{3}} \frac{\Gamma\left(\frac{3}{2}\right)\Gamma\left(\frac{1}{2}\right)}{\Gamma(2)} \left\{ \psi_{0} (2) - \psi_{0} \left(\tfrac{3}{2} \right) \right\}
= \frac{8}{\pi^2} \int_{0}^{1} \frac{\sqrt{x} - x}{1 - x} \, dx
= \frac{8 (2 \log 2 - 1)}{\pi^2},
\end{align*}
where $\psi_0 (z) = \dfrac{\Gamma'(z)}{\Gamma(z)}$ is the digamma function, satisfying the following identity
$$ \psi_{0}(z+1) = -\gamma + \int_{0}^{1} \frac{1 - x^{z}}{1 - x} \, dx. $$
I usually don't answer my own questions, but I literally just derived the solution and I think it's beautiful.
So we know from the link in the OP that:
$$\zeta(s,a)=\frac{1}{\Gamma(s)} \int_0^\infty \frac{z^{s-1} dz}{e^{a z} (1-e^{-z})}$$
Which makes our expression:
$$\zeta \left(\frac12, \frac{t}{2}\right)-\zeta \left(\frac12, \frac{t+1}{2}\right)=\frac{1}{\sqrt{\pi}} \int_0^\infty \frac{e^{- t z/2}(1-e^{-z/2}) dz}{(1-e^{-z}) \sqrt{z}}=$$
$$=\frac{2}{\sqrt{\pi}} \int_0^\infty \frac{e^{- t u^2/2}(1-e^{-u^2/2}) du}{1-e^{-u^2}}$$
Now we can see that it's very easy to take integral over $t$ (integration by parts):
$$\int_0^1 \sin \pi t~ e^{- t u^2/2} dt=\frac{4 \pi}{4 \pi^2 +u^4} (1+e^{-u^2/2})$$
Now we substitute this into the second integral to get (this is the beautiful part):
$$8 \sqrt{\pi} \int_0^\infty \frac{(1+e^{-u^2/2})(1-e^{-u^2/2}) du}{(1-e^{-u^2})(4 \pi^2 +u^4)}=8 \sqrt{\pi}\int_0^\infty \frac{du}{4 \pi^2 +u^4}$$
After a change of variables we have:
$$\frac{2 \sqrt{2}}{\pi} \int_0^\infty \frac{dv}{1 +v^4}=\frac{2 \sqrt{2}}{\pi} \frac{\pi}{2 \sqrt{2}}=1$$
Just as it was supposed to be.
God, Mathematics is simply perfect sometimes.
Appendix
Trying to prove in a simple way that the integral formula for $\zeta(1/2,a)$ equals the series definition.
$$ \int_0^\infty \frac{z^{-1/2} dz}{e^{a z} (1-e^{-z})}=2 \int_0^\infty \frac{e^{-a u^2}du}{ 1-e^{-u^2}}=2 \sum_{n=0}^\infty \int_0^\infty e^{-(a+n) u^2} du=$$
$$=2 \sum_{n=0}^\infty \frac{1}{\sqrt{a+n}} \int_0^\infty e^{-w^2} dw= \sqrt{\pi} \sum_{n=0}^\infty \frac{1}{\sqrt{a+n}} $$
Formally, this exactly fits the series definition, but it doesn't converge (it's alright, as the function for $s=1/2$ is defined by analytic continuation).
On the other hand, for the particular function in my case, the proof works well, as the series inside the integral becomes alternating due to a factor $(1-e^{-u^2/2})$ in the numerator, so everyting converges.
I would say that all of this constitutes a nice real proof of the Fresnel integrals, especially since the Poisson integral also has a few real proofs.
Though I'm definitely not claiming this is a new result. It was new for me, but a two second google search found me this paper https://www.jstor.org/stable/2320230, and I'm sure there's plenty more.
Best Answer
Let's consider the legality of doing an actual u-substitution, such as $z = \sqrt{i} x$. Not only must the integrand be rewritten, so must the limits of integration.
In the original definite integral you have $x$ going from $0$ to $\infty$. Of course this then gives a path of integration for $z$, but it's not sufficient to have just limits $0$ to $\infty$ in the complex plane to specify that path. So this would be a gray area where the limitations of your notation could let you down!
In the complex plane there are many paths from $0$ to $\infty$, even many straight such paths.