This approach can (and likely will) produce values that are impossible, and there are other methods like logistic regression that avoid this issue. That’s one reason why OLS might not be preferred for such a situation.
However, the math works out fine. Nothing will stop you from fitting an OLS model with your parameters estimated the usual way with $\hat\beta=(X^TX)^{-1}X^Ty$. If such a model works best for what you’re doing, as proponents of linear probability models believe, then go for it! Nothing about the derivation of $\hat\beta=(X^TX)^{-1}X^Ty$ makes a distribution assumption about the error term.
For instance, in the following simulation, the linear probability model has tight confidence intervals around the true parameters (intercept is $0$, slope is $1$), and that even happens with a small sample size of just ten.
set.seed(2023)
R <- 1000
N <- 10
intercepts <- slopes <- rep(NA, R)
for (i in 1:R){
x <- runif(N, 0, 1)
y <- rbinom(N, 1, x)
L <- lm(y ~ x)
intercepts[i] <- summary(L)$coef[1, 1]
slopes[i] <- summary(L)$coef[2, 1]
}
t.test(intercepts)$conf.int[c(1, 2)] # Confidence interval (-0.01600467, 0.01664349)
t.test(slopes, mu = 1)$conf.int[c(1, 2)] # (0.9723355, 1.0277178)
If the true relationship is $\mathbb E\left[Y\vert X\right] = \beta_0 + \beta_1X$, as is the case here, then the linear model captures this relationship.
As with Frank Harrell in the comments, I have my doubts about how well you will get a linear probability model to fit most data, but if it works...