The pdf of a bivariate normal is
$$f(x_1,x_2) = \frac{1}{2\pi \det|\Sigma|^{1/2}} e^{-\frac{1}{2}\big((x-\mu)^T \Sigma^{-1}(x-\mu)\big)}$$
Which can be written as:
$$f(x_1,x_2) = \propto \exp\bigg[-\frac{1}{2}\bigg({(\frac{x_1-\mu_1}{\sqrt{\sigma_{11}}})^2}+ (\frac{x_2-\mu_2}{\sqrt{\sigma_{22}}})^2 -2\rho_{12}(\frac{x_1-\mu_1}{\sqrt{\sigma_{11}}})(\frac{x_2-\mu_2}{\sqrt{\sigma_{22}}})\bigg)\bigg]$$
From my lecture notes the contours are all ellipses because the exponent can be factorized into the formula for an ellipse.
$$\frac{(x-h)^2}{a^2} + \frac{(y-k)^2}{b^2} = c^2$$
I dont understand how the exponent can be factorized into the above form if the random variables are not independent. What happens to the term in $x_1x_2$
You're right. The contours form ellipses all right, but with the correlation coefficient term, they are rotated with respect to the axes, which the generic ellipse formula you list doesn't cover. You would need to perform a linear transformation to rotate axes to get it into the form shown.