Suppose we have a linear model $Y = X\beta + \varepsilon$, where $Y \in \mathbb{R}^n, \beta \in \mathbb{R}^p, X \in \mathbb{R}^{n \times p},$ and $\varepsilon \sim N_n(0,\sigma^2I_n)$ (multivariate normal). Are the residuals $e$ and fitted values $\hat{Y}$ uncorrelated?

When tackling this question, start by expressing both the residuals and fitted values in terms of the observed data, the model, and error terms. Remember that the residuals are the differences between observed and fitted values, and the fitted values are derived from the projection of Y onto the column space of X. Utilize your knowledge of linear algebra and properties of matrix operations to prove or disprove the correlation between residuals and fitted values.