Help clarify the implication of linearity in an Ordinary Least Squares (OLS) Regression

If the dependent variable is linearly related to the independent variables, there should be no systematic relationship between the residuals and the fitted values. In other words, the model should capture all the systematic variance present in the data, leaving nothing but random noise.

I don't understand the implication in bold. I don't see how the conclusion (no systematic relationship between the residuals and the fitted values) follows from the hypothesis (the dependent variable is linearly related to the independent variables). Can someone please explain it in the context of a simple linear regression?

This is an abuse of notation, but you can always write $Y \,\vert\, X = \mathbb{E}\left[Y\,\vert\, X\right] + \epsilon \,\vert\, X$ where $\mathbb{E}\left[\epsilon \,\vert\, X\right]$ is necessarily zero. (If it were not , you could take expectations on both sides and get $\mathbb{E}\left[Y\,\vert\, X\right] \neq \mathbb{E}\left[Y\,\vert\, X\right]$.)

In words: conditional on X, Y equals its conditional mean plus some mean-zero residual (whose distribution might nonetheless depend on X, e.g. the variance of epsilon might change with X).

I interpret "the dependent variable is linearly related to the independent variables" in your quote as meaning $\mathbb{E}\left[Y\,\vert\, X\right] = \beta^{T}X$, i.e. the mean of $Y$ given $X$ is linear in $X$.

If that assumption is true, the model residuals $\hat{\epsilon} \equiv Y - \hat{\beta}^{T}X$, where $\hat{\beta}$ is your estimate of $\beta$, should have (approximately) zero mean for all values of $X$, and hence for all values of $\hat{\beta}^{T}X$. In other words, you're looking to see whether $\mathbb{E}\left[Y - \beta^{T}X \,\vert\, X\right] = 0$ for all $X$, using the estimated coefficients since you don't know the true $\beta$.

Perhaps the point is made most clearly if I write: $\mathbb{E}\left[Y - \mathbb{E}\left[Y\,\vert\, X\right] \,\vert\, X\right] = 0$ always holds. Linearity of the conditional mean implies $\mathbb{E}\left[Y - \beta^{T}X \,\vert\, X\right] = 0$ for all values of $X$, which implies there is "no systematic relationship between the residuals and the fitted values" -- if you interpret "systematic relationship" as meaning a change in conditional mean.

An example where the linearity assumption is violated: suppose $Y \,\vert\, X = \sin X + \epsilon \,\vert\, X$ where $\epsilon \,\vert\, X$ has mean zero:

n <- 10^3
df <- data.frame(x=runif(n, 1, 10))
df$mean.y.given.x <- sin(df$x)
df$y <- df$mean.y.given.x + rnorm(n)
model <- lm(y ~ x, data=df)
plot(df\$x, residuals(model))  # Looks like the mean of the residuals varies with X
plot(predict(model, newdata=df), residuals(model))  # Exact same issue