As you said, the inequality holds for the orthonormal polynomials $P_{k}$, $k\geq0$, of any positive measure $\mu$, with support $K$, a compact subset of $\mathbb{C}$. It follows from a basic result about the Christoffel-Darboux kernel $K_{N}(w,z)$ and the Christoffel function $\lambda_{N}(z)$, associated to $\mu$, respectively defined by
$$
K_{N}(w,z)=\sum_{k=0}^{N-1}P_{k}(w)\overline{P_{k}(z)},\qquad
\lambda_{N}(z)=\min_{\deg P\leq N-1,~P(z)=1}\int_{K}|P(w)|^{2}d\mu(w).$$
Theorem :
One has
$$
\lambda_{N}(z)=\frac{1}{\sum_{k=0}^{N-1}|P_{k}(z)|^{2}}=\frac{1}{K_{N}(z,z)},\qquad z\in\mathbb{C},$$
and the minimum is attained (uniquely) when
$$
P(w)=\Pi_{N,z}(w)
:=K_{N}(w,z)/K_{N}(z,z).
$$
For completeness, here is the proof :
let $P(w)=\sum_{k=0}^{N-1}\alpha_{k}P_{k}(w)$ be any polynomial of degree $\leq N-1$. Then
$$
\lambda_{N}(z)^{-1}=\max_{\alpha_{k}}
\frac{|\sum_{k=0}^{N-1}\alpha_{k}P_{k}(z)|^{2}}
{\sum_{k=0}^{N-1}|\alpha_{k}|^{2}}\leq
\sum_{k=0}^{N-1}|P_{k}(z)|^{2},
$$
and the Cauchy-Schwarz inequality becomes an equality when
$\alpha_{k}=\overline{P_{k}(z)}$, $k=0,\ldots,n$, which also gives the expression for the minimizing polynomial $\Pi_{N,z}(w)$.
Making use of the above result, we have, for any subset $K_{0}$ of $K$,
\begin{align*}
\int_{K_{0}}\int_{K_{0}} |K_{N}(w,z)|^{2}d\mu(z)d\mu(w)
& =\int_{K_{0}}K_{N}(z,z)^{2}\int_{K_{0}}|\Pi_{N,z}(w)|^{2}d\mu(w)d\mu(z)
\\[5pt]
& \leq\int_{K_{0}}K_{N}(z,z)^{2}\lambda_{N}(z)d\mu(z)
=\int_{K_{0}}K_{N}(z,z)d\mu(z).
\end{align*}
By the way, I think there is a missing normalization of the variables $x$ and $y$ in the Hermite polynomials $H_{l}$, in your definition of $f_{N}(x,y)$.
