I am struggling to understand the final step of the below derivation from my machine learning textbook (Section 3.6 of Machine Learning: A Bayesian and Optimization Perspective, 2e).
Let $y_n=\theta+\eta_n$ where $\theta$ is some constant and $\eta_n$ are $N$ noise samples from a normal distribution with mean $0$ and variance $\sigma_{\eta}^2$, and $\bar{y}$ represent the sample mean. The variance of the estimator $\bar{y}$ is derived as follows:
\begin{align} \sigma^2_{\bar{y}} & = \mathbb{E}\left[(\bar{y}-\theta)^2\right] = \mathbb{E}\left[\frac{1}{N^2}\left(\sum_{n=1}^N (y_n-\theta)\right)^2\right] \newline & = \frac{1}{N^2}\mathbb{E}\left[\left(\sum_{n=1}^N \eta_n\right)^2\right] = \frac{1}{N^2}\mathbb{E}\left[\sum_{i=1}^N \eta_i \sum_{j=1}^N \eta_j \right] \newline &= \frac{1}{N^2}\sum_{i=1}^N \sum_{j=1}^N \mathbb{E}[\eta_i \eta_j] = \frac{\sigma_{\eta}^2}{N} \end{align}
I follow up until $$\frac{1}{N^2}\sum_{i=1}^N \sum_{j=1}^N \mathbb{E}[\eta_i \eta_j] = \frac{\sigma_{\eta}^2}{N}$$ By what definition of variance does the above hold?
I know that $\textrm{Var}[x]=\mathbb{E}\left[\left(x-\mathbb{E}[x]\right)^2\right]=\mathbb{E}[x^2]-\mathbb{E}[x]^2$, and $\mathbb{E}[f(x)]=\lim_{N\rightarrow\infty} \frac{1}{N}\sum_{i=1}^N f(x_i)$, but neither of these seem to help here.
If $i \ne j$, if $\eta_i$ and $\eta_j$ are independent, then $E[\eta_i\eta_j]=E[\eta_i]E[\eta_j]=0\cdot 0=0$.
If $i=j$, then we have $e[\eta_i\eta_j]=E[\eta_i^2]=\sigma^2$.
Hence $$\frac{1}{N^2}\sum_{i,j}E[\eta_i\eta_j]=\frac1{N^2}\sum_i E[\eta^2]=\frac{N\sigma^2}{N^2}=\frac{\sigma^2}{N}.$$