Exercise :
Consider the Two-Way Anova model $Y_{ij} = \mu + a_i + b_j + e_{ij}$ with $i = 1, \dots, p$ and $j=1,\dots,q$. Show that : $$\text{cov}(\hat{e_{ij}},\hat{e_{i\ell}}) = -\sigma^2\left(\frac{1}{q} - \frac{1}{pq}\right), \quad j \neq \ell$$
Attempt/Thoughts :
So, this is part of an exercise involving the proofs of all potential covariances of errors. I have managed to prove the most straightforward case, which is :
$$\text{cov}(\hat{e_{ij}}, \hat{e_{ij}}) = V(\hat{e_{ij}}) = \sigma^2\left[1-\left(\frac{1}{p} + \frac{1}{q} - \frac{1}{pq}\right)\right]$$
But in this case, it's rather simple, since first of all you have the variance and secondly
$$\hat{e_{ij}} = Y_{ij} - \hat{Y_{ij}} $$ which makes all the work after that just calculations of Variances and some covariances.
But for this particular exercise listed, it will be :
$$\hat{e_{ij}} = Y_{ij} - \hat{Y_{ij}} \quad \text{and} \quad \hat{e_{i\ell}} = Y_{i\ell} - \hat{Y_{i\ell}}$$
So, that yields :
\begin{align*} \text{cov}(\hat{e_{ij}},\hat{e_{i\ell}})&=\text{cov}(Y_{ij} - \hat{Y_{ij}},Y_{i\ell} - \hat{Y_{i\ell}})\\ &= \mathbb{E}[(Y_{ij} - \hat{Y_{ij}})\cdot(Y_{i\ell} - \hat{Y_{i\ell}})] - \mathbb{E}[Y_{ij} - \hat{Y_{ij}}]\mathbb{E}[Y_{i\ell} - \hat{Y_{i\ell}}] \end{align*}
This expression seems much complex though. How would one proceed after that ?
Or if my approach thus far is wrong, how would one generally find the covariance asked ?