Suppose $y \sim N(\mu, \Sigma)$, where $y \in \mathbb{R}^n$. Let $X \in \mathbb{R}^{n \times p}$ denote a full rank design matrix. By ordinary least squares, the residuals are $$\hat{e} = (I - X(X^TX)^{-1}X^T)y$$ Let $u_1 \in \mathbb{R}^p$ denote a unit vector in $\mathbb{R}^p$. Let $\hat{\beta}_1 = u_1^T (X^TX)^{-1}X^Ty$ denote the OLS estimate for the first covariate in $X$. Suppose I'm interested in the condition distribution of $$\hat{\beta}_1 | \hat{e} = e_0$$
My question is, once I condition on a specific set of the residuals $\hat{e} = e_0$, is there still randomness left in $\hat{\beta}_1$? I'm confused about this because $\hat{\beta}_1$ is a function of $y$, so any variability in $\hat{\beta}_1$ stems from the fact that $y$ is random. Since $\hat{e}$ is a function of $y$, does conditioning on a specific value of $\hat{e}$ imply that there's no more variability left in $\hat{\beta}_1$? In other words, is it possible for different $y$ vectors to give rise to the same residuals, $e_0$?