$R^2$ of regression of $\epsilon$ on $X$ (residuals on independent variables) Assume you have a linear model.
What would happen if one regressed the OLS residuals $\epsilon = y − X \hat\beta$ on $X$?
Would it yield by construction an $R^2$ of zero? If so, why?
 A: Yes, it will, always, regardless whether we have specified our model correctly or not.
Before explaining why, let's first see some examples. I don't know what statistical software you use, but consider the following example using R.
## We consider a true model: a quadratic polynomial on [0, 1]
set.seed(0); x <- sort(runif(50))  ## 50 sampling points 
y <- rnorm(50, mean = 0.1 + 0.2 * x + 0.5 * x ^ 2, sd = 0.1)  ## observations

We first fit our model correctly.
## fit the model correctly
m0 <- lm(y ~ poly(x, 2, raw = TRUE))
r0 <- residuals(m0)  ## take out residuals
mr0 <- lm(r0 ~ poly(x, 2, raw = TRUE))  ## fit residuals to the same RHS
coef(mr0)  ## extract coefficients
# (Intercept) poly(x, 2, raw = TRUE)1 poly(x, 2, raw = TRUE)2
#1.204349e-18           -4.702044e-18            4.373216e-18

On the first inspection, we see that all coefficients are estimated to be 0. This means the fitted values of this residual model will be 0. Thus, if we compute multiple $R^2$, we get 0:
drop(crossprod(fitted(mr0))) / drop(crossprod(r0))
# [1] 0

Now, we fit a wrong model, by just fitting a line.
## fit the model incorrectly
m1 <- lm(y ~ x)
r1 <- residuals(m1)  ## take out residuals
mr1 <- lm(r1 ~ x)  ## fit residuals to the same RHS
coef(mr1)  ## extract coefficients
#  (Intercept)             x 
# 4.140599e-18 -1.201517e-17 

Again, all coefficients of the residual model are estimated to be 0. So once again, we will have zero fitted values and thus zero multiple $R^2$.

So why? Because the least square fitting of $y = X\beta + e$ ensures that residuals $r = y - X\hat{\beta}$ are orthogonal to design matrix $X$, i.e., $r'X = X'r = 0$. This is pretty straightforward to prove using that $\hat{\beta} = (X'X)^{-1}X'y$, and I am sure there are questions on this so I will not repeat. But it implies that if we regress $r$ onto $X$, we get zero estimate: $(X'X)^{-1}X'r = (X'X)^{-1}0 = 0$, hence zero fitted values, therefore zero regression sum of squares, and in the end zero $R^2$.
A: In my opinion the answer to this question realies heavily on the definition of $R^2$. Many books for instance define $R^2$ as $$R^2 = \frac{y'(P-Q)y}{y'(I-Q)y}$$
where $P=X(X'X)^{-1}X'$ and $Q=\frac{1}{n}1_n1_n'$. Now plugging in the residual vector instead of $y$ yields
$$R^2 = \frac{-\hat\epsilon'Q\hat\epsilon}{\hat\epsilon'(I-Q)\hat\epsilon}$$
which is negative unless the model contains an intercept (in the answer above using R an intercept is included. Note that R automatically caluclates the uncentered $R^2$ when not including a constant).
