Why are residuals Normally distributed? When they assume the $\varepsilon_i$ have a normal distribution, with mean $0$ and variance $\sigma^2,$ then the estimated errors $\hat\varepsilon$ have a normal distribution with mean $0$ and variance $\operatorname{Var}(\hat\varepsilon).$
But what is the connection? Why, when errors are normally distributed, the estimated residuals also are normally distributed?

The connection is made through a specific procedure: namely, ordinary least squares.  When you suppose a response $Y$ is a random variable related to other variables $X$ in the form
$$Y = X\beta + \varepsilon$$
where $\varepsilon$ is assumed to have a joint Normal distribution (plus some other assumptions that don't matter here), then the least squares solution is
$$\hat\beta = (X^\prime X)^{-}X^\prime Y = (X^\prime X)^{-}X^\prime (X\beta + \epsilon) = \beta + J\varepsilon$$
for a matrix $J = (X^\prime X)^{-}X^\prime $ that is computed only from $X.$  The residuals, or estimated errors, are the differences
$$\hat\varepsilon = Y - \hat Y = (X\beta + \varepsilon) - X\hat\beta = \varepsilon - XJ\varepsilon = (\mathbb I - H)\varepsilon$$
where $H = XJ = X(X^\prime X)^{-}X^\prime$ is the "hat matrix."  This exhibits the residuals as linear combinations of $\varepsilon$ (with coefficients given by $\mathbb I - H$).  Linear combinations of jointly Normal variables are Normal, QED.
