In THIS ANSWER, I used straightforward complex analysis to show that
$$\gamma =2\int_0^\infty \frac{\cos(x^2)-\cos(x)}{x}\,dx \tag 1$$
where $\gamma =-\int_0^\infty \log(x) e^{-x}\,dx$ is the Euler-Mascheroni Constant.
The key in the derivation of $(1)$ was to transform the cosine terms into real exponential ones.
To date, I have been unable to use strictly real analysis, without appealing to tabulated results of special functions (e.g., use of the $\text{Cin}(x)$ and $\text{Ci}(x)$ functions), to prove $(1)$.
I have tried introducing a parameter and using "Feynman's Trick to augment the integral into something manageable. Or somewhat equivalently, rewriting the integral in $(1)$ as a double integral and proceeding by exploiting Fubini-Tonelli.
QUESTION: What are ways to prove $(1)$ without relying on complex analysis and without simply appealing to tabulated relationships of special functions. For example, stating that the $Ci(x)$ function is defined as $\text{Ci}(x)\equiv -\int_x^\infty\frac{\cos(t)}{t}\,dt=\gamma +\log(x) +\int_0^x \frac{\cos(t)-1}{t}\,dt$ is unsatisfactory unless one proves the latter equality.
It turns out that we have the following observation:
The reasoning is surprisingly simple: First, define $g(x) = (f(x) - f(0)\mathbf{1}_{(0,1)}(x))/x$ and notice that
$$ \int_{\epsilon}^{\infty} \frac{f(x)}{x} \, dx = -f(0)\log\epsilon + \int_{\epsilon}^{\infty} g(x) \, dx. $$
Assuming that the LHS of $\text{(1)}$ exists for all $\epsilon > 0$ and that $f$ behaves nice near $x = 0$, this implies $\text{(1)}$. Next, notice that $c(f) = \mathcal{L}g(0)$ and that $-(\mathcal{L}g(s))' = \mathcal{L}f(s) - f(0) (1-e^{-s})/s$. Therefore
\begin{align*} c(f) &= \lim_{R\to\infty} \int_{0}^{R} (-\mathcal{L}g(s))' \, ds \\ &= \lim_{R\to\infty} \left( \int_{0}^{R} \mathcal{L}f(s) \, ds - f(0) (1 - e^{-R})\log R + f(0) \int_{0}^{R} e^{-s}\log s \, ds \right) \\ &= \lim_{R\to\infty} \left( \int_{0}^{R} \mathcal{L}f(s) \, ds - f(0) \log R \right) - f(0)\gamma. \end{align*}
At this moment this is just a heuristic computation. For a broad class of functions for which the LHS of $\text{(1)}$ exists, however, this computation can be made rigorous. This is particularly true for our function $f(x) = \cos x$. Now plugging $\mathcal{L}f(s) = \frac{s}{s^2+1}$ shows that $c(f) = -\gamma$ and thus
$$ \int_{\epsilon}^{\infty} \frac{\cos x}{x} \, dx = -\log\epsilon - \gamma + o(1). $$
Plugging this asymptotics, we have
$$ \int_{\epsilon}^{\infty} \frac{\cos(x^2) - \cos x}{x} \, dx = \frac{1}{2}\int_{\epsilon^2}^{\infty} \frac{\cos x}{x} \, dx - \int_{\epsilon}^{\infty} \frac{\cos x}{x} \, dx = \frac{1}{2}\gamma + o(1) $$
and the identity follows by letting $\epsilon \to 0^+$.
Here, the constant $c(f)$ can be thought as a regularized value of the divergent integral $\int_{0}^{\infty} \frac{f(x)}{x} \, dx$. This has the following nice properties (whenever they exist)
Together with some known values, we can easily compute other types of integrals. For instance, using the fact that $c(\cos x) = -\gamma$ and $c(e^{-x}) = -\gamma$, we have
\begin{align*} \int_{0}^{\infty} \frac{\cos (x^p) - \exp(-x^q)}{x} \, dx &= c\left\{ \cos (x^p) - \exp(-x^q) \right\} \\ &= \frac{1}{p}c(\cos x) - \frac{1}{q}c(e^{-x}) = \gamma\left( \frac{1}{q} - \frac{1}{p}\right) \end{align*}
for $p, q > 0$.