Consider the (physicist's) Hermite polynomials $H_n(x)$ which are divided by $$\sqrt{\sqrt{\pi} 2^n n!}$$ for the purpose of normalization.

These are orthogonal with respect to the weight function $e^{-x^2}$ and $\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.