Show $\mathrm{cov}(E_i,E_j)=\sigma^2p_{ij}$ where $E_i=Y_i-\bar y - \hat\beta(x_i-\bar x)$ The model is given by $E(Y_i)=\gamma+\beta (x_i-\bar x)$

Show $\mathrm{cov}(E_i,E_j)=\sigma^2p_{ij}$ where $$p_{ij}=1/n+\frac{(x_i-\bar x)(x_j-\bar x)}{S_{xx}}$$

I already showed that $E(E_i)=0$ and $\mathrm{cov}(Y_i,\bar Y)=\sigma^2/n$
I think I'm missing some easy properties that I can use.
I can see that somewhat related (not equal) $\sigma^2 p_{ij}=\mathrm{var}(\hat Y)$
I tried:
$$\mathrm{cov}(E_i,E_j)=\mathrm{cov}(Y_i-\bar y - \hat\beta(x_i-\bar x),Y_j-\bar y - \hat\beta(x_j-\bar x))$$
But then I don't know how to procceed because I have $3$ random variables. I can find $E(Y_i)=\gamma+\beta (x_i-\bar x),E(\bar y)=\gamma$ and $E(\hat\beta)=\beta$

It is not a complete answer, i.e., it seems that there is no way to avoid some tedious calculations, but may be it can help you:
1) $ cov(\hat{\beta}(x_i - \bar{x}), \hat{\beta}(x_j - \bar{x}) ) = \sigma^2\frac{(x_i-\bar{x})(x_j-\bar{x})}{\sum (x_i - \bar{x})^n}$.
2) $cov(\bar{Y}_n, \bar{Y}_n) = var(\bar{Y}_n)=\sigma^2/n $
3) If the error terms are i.i.d then $cov(Y_i, \bar{Y}_n) = cov(Y_i, Y_i/n)=\sigma^2/n $  
