How exactly is OLS derived from assumption of normally distributed residuals? Ordinary least square solution of linear regression can be derived from the assumption of normally distributed residuals:
$$
e_i=y_i-\hat{y_i}\\
e_i\sim N(0, \sigma^2)
$$
What I don't quite understand from reading a lot of derivations is how we go from
$MLE=\arg \max_{\theta}\sum{log(p_{\theta}^i)}$ to the sum of logarithms of PDFs of normal distribution. PDFs are called "probabilities" everywhere I read (e.g. https://machinelearningmastery.com/linear-regression-with-maximum-likelihood-estimation), which is incorrect. PDFs can be > 1. Only integrals of PDFs can be interpreted as probabilities of X falling into the integration interval.
Do I understand correctly, that the proper derivation should introduce a small interval $\epsilon$, so that we could first rewrite probabilities as:
$$
p_{\theta,\epsilon}^i=\int_{y_i-\epsilon}^{y_i+\epsilon}PDF(N(\hat{y^i}, \sigma^2))dy
$$
Then we could assume that with small enough $\epsilon$
$$
p_{\theta,\epsilon}^i\approx 2*\epsilon*[\frac 1 {\sqrt {2*\pi}*\sigma} e^{-(\hat{y_\theta^i} - y_i)^2/\sigma^2}]
$$
Now MLE can finally be approximated with
$$
\require{cancel}
MLE=\arg\max_{\theta}\sum log(2*\epsilon*[\frac 1 {\sqrt {2*\pi}*\sigma} e^{-(\hat{y_\theta^i} - y_i)^2/\sigma^2}])\\
=\arg\max_{\theta}\sum \cancel{log(2*\epsilon)} + \cancel{log(\frac 1 {\sqrt {2*\pi}*\sigma})} -(\hat{y_\theta^i} - y_i)^2/\sigma^2\\
=\arg\min_{\theta}\sum (\hat{y_\theta^i} - y_i)^2
$$
And the last expression is the OLS.

You seem to mix two issues: first, how to interpret maximum likelihood estimation and second the derivation of OLS as an ML estimator.
Regarding the first point: Leaving mathematical rigour aside and thinking about the likelyhood about the joint probabilities for all possible small intervals can be useful. Let $\mathbf{y}=\{y_1, y_2..., y_n\}$ the data we observe. For each datapoint, we want to choose an estimator which maximises the probability to observe data in a small neighbourhood of $y_1$ and in a small neighbourhood of $y_2$ and in a small neighbourhood of $y_3$ and so on...  We impose that these probabilites are compatible with the assumption of all $e_i=y_i-\hat{y}_i$ being drawn form a normally distributed random variable.
Regardingthe second point, i.e. estimation I agree with seanv507: Following wikipedia notation you are maximising the likelihood function which is the product of individual densities:
\begin{equation}
f_n(\mathbf{y,\theta})=\Pi_{k=1}^n f_k^{univar}(y_k;\theta)
\end{equation}
plugging in univariate normal densities, taking logs and simplifying along the lines you do gives you the desired result.
