Adding a variable with no relationship with $Y$ increases $R^2$ I'm confused by something I've found by adding a variable with no relationship with a DV, using multiple regression with four predictors and one DV (Y). If I regress $Y$ onto $X_1,~ X_2,$ and $X_3$, the multiple $R$ is less than if I add a 4th predictor with no relationship with $Y$. I didn't think this was possible. I've done this via a simulation and also more manually, with each shown below. What's even more confusing is that if I specify the 4th variable to have a correlation of $.2$ with the DV, the $R^2$ is less than if the 4th variable has a correlation of $0$ with the DV. How is this possible?
### via simulation ###

library(MASS)
library(psych)

rx12 = .2
rx13 = .25
rx14 = .3
rx23 = .35
rx24 = .3
rx34 = .4

rx1y = .15
rx2y = .25
rx3y = .2
rx4y = 0

corr_matrix  <- matrix(c(1, rx12, rx13, rx14, rx1y, 
                         rx12, 1, rx23, rx24, rx2y,
                         rx13, rx23, 1, rx34, rx3y,
                         rx14, rx24, rx34, 1, rx4y,
                         rx1y, rx2y, rx3y, rx4y, 1), nrow=5)
corr_matrix #this shows the correlation is zero#

set.seed(33)
data = as.data.frame (mvrnorm(n=1000, mu=c(.0, .0, .0, .0, 0), Sigma=corr_matrix, empirical=TRUE)) 
psych::corr.test(data)$r #this shows the correlation is zero#

summary(lm(V5 ~ V1 + V2 + V3, data=data)) #R^2 = .0833
summary(lm(V5 ~ V1 + V2 + V3 + V4, data=data)) #R^2 = .1044

### matrix multiplication with all 4 variables ###
corr_matrix_x  <- matrix(c(1, rx12, rx13, rx14,
                           rx12, 1, rx23, rx24, 
                           rx13, rx23, 1, rx34,
                           rx14, rx24, rx34, 1), nrow=4)
corr_matrix_y  <- matrix(c(rx1y, rx2y, rx3y, rx4y), nrow=4)
corr_matrix_y #this shows the correlation is zero#

x_inverse <- solve(corr_matrix_x)
betas <- as.matrix(x_inverse %*% corr_matrix_y)
t(betas) %*% corr_matrix_y   #R^2 = .1044

### 3 variables ###
corr_matrix_x  <- matrix(c(1, rx12, rx13, 
                           rx12, 1, rx23, 
                           rx13, rx23, 1), nrow=3)
corr_matrix_y  <- matrix(c(rx1y, rx2y, rx3y), nrow=3)

x_inverse <- solve(corr_matrix_x)
betas <- as.matrix(x_inverse %*% corr_matrix_y)
t(betas) %*% corr_matrix_y #R^2 = .0833

Follow-up:
My friend visually examined the issue of R^2 being larger when the X4 correlation with Y is zero or small (~ <.2), and then only increasing at correlations ~>.2. Image below:


You should consider the so-called Adjusted R-squared that imposes a penalty for considering additional independent variables to a model. The Un-adjusted  R-squared can never fall when a new independent variable is added to a regression equation: this is because Sum of Squared Residuals (SSR) never goes up (and usually falls) as more independent variables are added (assuming we use the same set of observations). Adjusted R-squared can be calculated as follows:
$$
\overline{R}^2=1-\frac{(1-R^2)(n-1)}{n-k-1}
$$
where $R^2$ is the Un-adjusted R-squared, $k$ is the number of independent variables and $n$ the sample size.
