Why a regression of OLS residuals on regressors, yields a $R^2$ of 0? In a regression of residuals
on X, the OLS estimator would be $\hat{\gamma} = (X'X)^{−1}X'\hat{\varepsilon}$ which is a vector of zeros ($\hat{\gamma} = 0$)
since $X'\hat{\varepsilon} = 0$. Could anyone tell me how did we arrive to the above conclusion about the OLS estimator?

$$
\hat\beta_{ols} = \left(X^TX\right)^{-1}X^Ty\\\Downarrow\\
\left(X^TX\right)\hat\beta_{ols} = X^Ty\\\Downarrow\\
\left(X^TX\right)\hat\beta_{ols} - X^Ty = 0\\\Downarrow\\
X^T\left(X\hat\beta_{ols} - y\right) = 0\\\Downarrow\\
X^T\hat\epsilon = 0
$$
(When you go through the derivation of the OLS estimator, $\hat\beta_{ols}$, it is the second line that implies the first. Let's set that aside.)
Since $X^T\hat\epsilon = 0$, then $\left(X^TX\right)^{-1}X^T\hat\epsilon = \left(X^TX\right)^{-1}0 = 0$.
To get into why $R^2=0$ for this regression on the residuals, note that this means all of the coefficients (including the intercept) are zero. Also note that the mean of the residuals is zero. As $R^2$ compares the performance of your model (which always predicts zero) to the performance of a model that always predicts the mean of the outcome (which will be the residuals in the case, so a mean of zero), the two models have identical performance.
Take the $y_i$ to be observed residuals from the original model, $\hat y_i$ to be the predictions made by regressing on the residuals from the original model, and $\bar y$ to be the mean of the residuals (which will be zero).
$$
R^2 = 1-\left(
\dfrac{
\underset{i = 1}{\overset{N}{\sum}}\left(
y_i -\hat y_i
\right)^2
}{
\underset{i = 1}{\overset{N}{\sum}}\left(
y_i -\bar y
\right)^2
}
\right)
= 1-\left(
\dfrac{
\underset{i = 1}{\overset{N}{\sum}}\left(
y_i -0
\right)^2
}{
\underset{i = 1}{\overset{N}{\sum}}\left(
y_i -0
\right)^2
}
\right) = 1-1 = 0
$$
This is reflected in an R simulation.
set.seed(2023)
N <- 10                 # sample size
X <- rnorm(N)           # simulate a feature
Y <- X + rnorm(N)       # simulate a response
L <- lm(Y ~ X)          # fit a regression
y <- resid(L)           # extract the residuals of that regression
new_L <- lm(y ~ X)      # regress those residuals on the original feature
y_hat <- predict(new_L) # make predictions from that second regression "new_L"

1 - (sum((y - y_hat)^2))/(sum((y - mean(y))^2)) # This is R^2. I get 0.
```

