Consider the (physicist's) [Hermite polynomials][1] $H_n(x)$ which are divided (for normalization purposes) by
$$\sqrt{\sqrt{\pi} 2^n n!}.$$
They become orthogonal with respect to the weight function $e^{-x^2}$ such that $\int_{\mathbb{R}}H_n(x)^2\,e^{-x^2}\,dx=1$.
Define the family of functions
$$f_N(x,y)=e^{-\frac{N}2(x^2+y^2)}\sum_{\ell=0}^{N-1}H_{\ell}(x)H_{\ell}(y).$$
Using tools from Random Matrix Theory, the following when treated as the variance of the number of eigenvalues in an interval, it is proven positive. However, I would like to ask:

>**QUESTION.** Is there a direct proof of the below inequality that does not involve RMT?
$$\int_a^bf_N(x,x)\,dx-\int_a^b\int_a^bf_N(x,y)^2dxdy>0 \qquad \text{for all $N\geq1$}$$

**Remark.** It turns out that the above inequality is in fact valid for other orthogonal polynomials.

[1]: https://en.wikipedia.org/wiki/Hermite_polynomials