Omitted Variable Bias (OVB) and multicollinearity In a linear regression model, the reason we control for variables is to prevent the omitted variable bias (OVB). That is, suppose we are trying to fit the model
$$
Y = \beta_{0} + \beta_{1}X_{1} + \varepsilon
$$
however, there is another variable $X_{2}$ that is correlated with $X_{1}$ and influences $Y$, then the estimate for $\beta_{1}$, which we can call $\hat{\beta_{1}}$, will be biased. The way to remove this bias is to add variable $X_{2}$ in our model, i.e. fit the model
$$
Y = \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + \varepsilon
$$
However, my question is why this doesn't violate the assumption that covariates in a multiple linear regression model shouldn't be highly correlated with one another (no multicollinearity)? What if $X_{1}$ and $X_{2}$ are multicollinear? In that case, omitting $X_{2}$ might lead to OVB, but including $X_{2}$ might lead to multicollinearity. What is done in cases like this?
 A: This is a good question. The confusion stems from the "assumption" of no multicollinearity. From the Wikipedia page on multicollinearity:

Note that in statements of the assumptions underlying regression analyses such as ordinary least squares, the phrase "no multicollinearity" usually refers to the absence of perfect multicollinearity, which is an exact (non-stochastic) linear relation among the predictors. In such case, the data matrix $X$ has less than full rank, and therefore the moment matrix $X^TX$ cannot be inverted. Under these circumstances, for a general linear model $y = X\beta + \epsilon$ , the ordinary least squares estimator $\hat\beta_{OLS} = (X^TX)^{-1} X^T y $ does not exist.

Multicollinearity in the sense that you describe will inflate the variance of the OLS estimator, but unless you include $X_2$ in the regression, the OLS estimator is biased. In short, if you have to worry about OVB, you should not be worrying about multicollinearity. Why would we want a more precise but biased estimator?
At more length, I am not sure that multicollinearity (or variance inflation) is at all meaningful to consider when we are concerned with OVB. Assume
$$ Y = 5X_1 + X_2 + \epsilon $$
$$ X_1 = -0.1X_2 + u $$
If $\text{Cov}(X_2, u) = 0$, the correlation between $X_1$ and $X_2$ is
$$ \rho = \frac{\sigma_{x_1x_2}}{\sigma_{x_1}\sigma_{x_2}} = \frac{-0.1\sigma_{x_2}}{\sqrt{0.01\sigma_{x_2}^2 + \sigma_u^2}} $$
If we let $\sigma_{x_2} = \sigma_{x_1}$, then $\rho \approx -0.1$ (which is a case where we would not worry about multicollinearity). Simulating in R, we see that an OLS regression of $Y$ on $X_1$ controlling for $X_2$ is unbiased. However, the bias that we get by excluding $X_2$ is pretty small.
iter <- 10000 # NUMBER OF ITERATIONS
n <- 100 # NUMBER OF OBSERVATIONS PER SAMPLE
sigma_e = sigma_u = sigma_x2 = 5
mu_e = mu_u = mu_x2 = 0
res0 = res1 = list() # LISTS FOR SAVING RESULTS

for(i in 1:iter) {
  
  #print(i)
  
  x2 <- rnorm(n, mu_x2, sigma_x2)
  u <- rnorm(n, mu_u, sigma_u)
  e <- rnorm(n, mu_e, sigma_e)
  
  x1 <- -0.1*x2 + u
  y <- 5*x1 + x2 + e
  
  res0[[i]] <-  lm(y ~ x1 + x2)$coef
  res1[[i]] <-  lm(y ~ x1)$coef
  
}

res0 <- as.data.frame(do.call("rbind", res0))
res1 <- as.data.frame(do.call("rbind", res1))


If we increase the variance of $X_2$ so that $\rho \approx -0.95$
sigma_x2 <- 150

and repeat the simulation we see that this does not affect the precision of the estimator for $X_1$ (but the precision for $X_2$ increases). However, the bias is now pretty big, which means that there is a big difference between the association between $X_1$ and and $Y$, where other factors (that is, $X_2$) are not held constant, and the effect of $X_1$ on $Y$ ceteris paribus. As long as there is some variation in $X_1$ that does not depend on $X_2$ (i.e., $\sigma_u^2 > 0$), we can retrieve this effect by OLS; the precision of the estimator will depend on the size of $\sigma_u^2$ compared to $\sigma_\epsilon^2$.

We can illustrate the effect of variance inflation by simulating with and without correlation between $X_1$ and $X_2$ and regressing $Y$ on $X_1$ and $X_2$ for both the correlated and uncorrelated case.
install.packages("mvtnorm")
library(mvtnorm)

sigma_x2 <- 5 # RESET STANDARD DEVIATION FOR X2
res0 = res1 = list()

Sigma <- matrix(c(sigma_x1^2, sigma_x1*sigma_x2*-0.95, 0,
                  sigma_x1*sigma_x2*-0.95, sigma_x2^2, 0,
                  0, 0, sigma_e^2), ncol = 3)

Sigma0 <- matrix(c(sigma_x1^2, 0, 0,
                   0, sigma_x2^2, 0,
                   0, 0, sigma_e^2), ncol = 3)


for(i in 1:iter) {
  
  print(i)
  
  tmp <- rmvnorm(n, mean = c(mu_x1, mu_x2, mu_e), sigma = Sigma0)
  x1 <- tmp[,1]
  x2 <- tmp[,2]
  e <- tmp[,3]
  
  y <- 5*x1 + x2 + e
  res0[[i]] <-  lm(y ~ x1 + x2)$coef
  
  tmp <- rmvnorm(n, mean = c(mu_x1, mu_x2, mu_e), sigma = Sigma)
  x1 <- tmp[,1]
  x2 <- tmp[,2]
  e <- tmp[,3]
  
  y <- 5*x1 + x2 + e
  res1[[i]] <-  lm(y ~ x1 + x2)$coef
  
}

res0 <- as.data.frame(do.call("rbind", res0))
res1 <- as.data.frame(do.call("rbind", res1))


This shows that the precision of the estimator would be better if $X_1$ and $X_2$ were uncorrelated, but if they are not, there is nothing we can do about it. It seems about as valuable as knowing that if our sample size were greater, then the precision would be better.
I can think of one example in which we could potentially care about both OVB and multicollinearity. Say that $X_2$ is a theoretical construct and you are unsure about how to measure it. You could use $X_{2A}$, $X_{2B}$, and/or $X_{2C}$. In this case, you might choose to just include one of theses measures of $X_2$ rather than all of them to avoid too much multicollinearity. However, if you are primarily interested in the effect of $X_1$ this is not a major concern.
