Let lowercase bold letters denote vectors. The linear model is:
$$ y_i = \mathbf{x}_i' \mathbf{b} + \epsilon_i $$
In matrix notation for all observations $i=1, \ldots, n$:
$$\mathbf{y} = X \mathbf{b} + \boldsymbol{\epsilon} $$
The OLS estimator for $\mathbf{b}$ is:
$$\hat{\mathbf{b}} = (X'X)^{-1}X'\mathbf{y} $$
Substituting:
$$\hat{\mathbf{b}} = (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) $$
Take the variance (conditional on X)$:
\begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= \operatorname{Var}\left( (X'X)^{-1}X' \left( X \mathbf{b} + \boldsymbol{\epsilon}\right) \mid X\right) \\ &= \operatorname{Var}\left( (X'X)^{-1}X' \boldsymbol{\epsilon} \mid X\right) \\ &= (X'X)^{-1}X' \operatorname{Var}\left( \boldsymbol{\epsilon} \mid X\right) X(X'X)^{-1} \end{align*}
We assumed that that $\epsilon_i$ were IID. Hence for all $i$ we have $\operatorname{Var}(\epsilon_i) = \sigma^2$, and for the whole vector $\boldsymbol{\epsilon})$ we have $\operatorname{Var}(\boldsymbol{\epsilon}) = \sigma^2 I$ where $I$ is the identity matrix.
Continuing: \begin{align*} \operatorname{Var}\left( \hat{\mathbf{b}} \mid X\right) &= (X'X)^{-1}X' \sigma^2 IX(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1}X'X(X'X)^{-1} \\ &= \sigma^2 (X'X)^{-1} \end{align*}
$\sigma^2$ is a scalar and can be moved wherever by the commutative property of multiplication. $(X'X)^{-1}X'X = I$ by definition of an inverse matrix.