In a Simple or Multiple Linear Regression Model, we are given the assumption that the error terms, $\epsilon$, are distributed like so:
$ε ∼ N(0, I_nσ^2).$
Can anybody explain where the mathematical logic comes from to derive this assumption, or what the intuition is behind it?
Some of the reasons for this choice, I would say, include
Intuitively reasonable: the Normal distribution is very nice, in that it is mean zero, symmetric (i.e. there is no bias direction in the error), has a memorable bell-shape, etc. All these properties make it a very "plausible" assumption for how errors would be distributed.
Mathematically convenient: this is essentially the backbone for a lot of statistical inference, including hypothesis tests, models, etc. For example, this allows you to deduce the distributions of different test statistics, sums, etc.
Of course, it is not the only possibility; after all, it is only an assumption. For example, generalisations such as GLMs allow for other error distributions of the response variable.