The entire discussion below, as well as Selberg's result, is under the assumption of the Riemann Hypothesis.
A proof is sketched in Soundararajan's expository notes "The distribution of prime numbers", see pages 9-10 in the arXiv version. The explanation there is lucid and there is nothing I can add.
Selberg's 1943 proof can be found in his collected papers, specifically it is paper no. 12 in his "Collected Works. Volume I" from 1989.
In page 10 of his paper (page 169 of the volume), you'll see that he essentially proves the following:
$$\int_{1}^{\infty} y^{-20/\log T}\left|\psi\left(y+\frac{y}{T}\right) -\psi(y)-\frac{y}{T}\right|^2 \frac{dy}{y^2}\ll \frac{\log^2 T}{T}$$
for all large enough $T$. If in Selberg's result you only integrate from $x$ to $2x$, and rename $y$ to $t$, you obtain the following direct consequence:
$$(\star)\,\int_{x}^{2x} \left|\psi\left(t + \frac{t}{T}\right) - \psi(t) - \frac{t}{T}\right|^2 \ll x \frac{x}{T} \log^2 x$$
holds as long as $\log x \asymp \log T$.
In the result you are after, $h$ plays the role of $t/T$ as well as of $x/T$. However, $t/T$ varies as $t$ varies, while $h$ does not, so setting $T:=t/h$ does not recover the result you desire.
Fortunately, Saffari and Vaughan, in Lemma 6 of "On the fractional parts of $\{x/n\}$ and related sequences. II" (1977), found a simple way to recover the result that you state as a direct consequence of Selberg's bound. Namely, Lemma 6 proves the bound you are after, using the following elementary general bound implicit in the proof of their lemma (see page 25):
$$\int_{X}^{2X} |F(x+H)-F(x)|^2 dx \ll \sup_{\theta \in \left[ \frac{H}{3X},\frac{3X}{H}\right]} \int_{X}^{3X} |F(u+\theta u)-F(u)|^2 du$$
for any $1\le H\le X$ and any square-integrable $F$. This immediately implies the result you ask about given Selberg's result.
Two remarks:
- At least as I wrote it, combining Selberg with Saffari-Vaughan recovers the result you desire only in the range $h\le x^{1-\varepsilon}$ for some $\varepsilon>0$. However, your result stays true for larger $h$, as observed by Saffari-Vaughan in their Lemma 5 (see the lemma and the comparison with Selberg immediately after its statement). The general and correct form of the result you seek, as stated in Lemma 6 of Saffari-Vaughan, is
$$(\star\star)\,\int_{1}^{x} \left|\psi(t+ h)-\psi(t)-h\right|^2 dt \ll x \left(1+\log^2(x/h)\right)$$
for $0< h \le x$, which they deduce from the following estimate (proved in Lemma 5): For $T\ge 1$ and $x\ge 4$, one has
$$\int_{x}^{2x} \left|\psi\left(t + \frac{t}{T}\right) - \psi(t) - \frac{t}{T}\right|^2 \ll x \frac{x}{T} \log^2 (2T)$$
without the restriction $\log T \asymp \log x$ I imposed when deducing $(\star)$ above.
- Shaving a log from the right-hand side of $(\star\star)$, even for some range of $h$, is in fact a big open problem, which is expected to be true in view of Montgomery's famous Pair Correlation Conjecture, which in particular implies that
$$\int_{1}^{x} \left|\psi(t+ h)-\psi(t)-h\right|^2 dt \sim x \log(x/h)$$
holds for $1\le h\le x^{1-\varepsilon}$ and $x\to \infty$. For more on this see Goldston and Montgomery's "Pair correlation of zeros and primes in short intervals" (1987). Section 9 of Goldston's lecture notes "Notes on Pair Correlation of Zeros and Prime Numbers" is also relevant.
