Understanding the math behind using maximum likelihood for linear regression
Note first that the first $=$ (equals) in $\frac{dl(\theta)}{d\theta} 
 = 0 = −\frac{1}{2\sigma^2}(0−2X^TY + X^TX \theta)$ should be interpreted as a "is set to", that is, we set $\frac{dl(\theta)}{d\theta} 
 = 0$. Given that (apparently) $\frac{dl(\theta)}{d\theta} = −\frac{1}{2\sigma^2}(0−2X^TY + X^TX \theta)$, $\frac{dl(\theta)}{d\theta} 
 = 0$ is equivalent to $0 = −\frac{1}{2\sigma^2}(0−2X^TY + X^TX \theta)$.

Now, let's apply some basic linear algebra:

\begin{align}
0 &amp;= −\frac{1}{2\sigma^2}(0−2X^TY + X^TX \theta) \iff \\
0 &amp;= −(0−2X^TY + X^TX \theta) \iff \\
0 &amp;= −0 + 2X^TY - X^TX \theta) \iff \\
0 &amp;= 2X^TY - X^TX \theta \iff \\
X^TX \theta &amp;= 2X^TY \iff \\
(X^TX)^{-1}(X^TX) \theta &amp;= (X^TX)^{-1}2X^TY \iff \\
\theta &amp;= (X^TX)^{-1}2X^TY
\end{align}

Now, you can ignore the $2$, because it is just a constant, and, when optimizing, this does not influence the result.

Note that using $\hat{\theta}$ instead of $\theta$ is just to indicate that what we will get is an "estimate" of the real $\theta$, because of round off errors during the computations, etc.