On page 599 of this book, the author states (without proving) that for random samples $(X_1, Y_1)$, ..., $(X_n, Y_n)$ from a bivariate normal distribution, the sample correlation coefficient \begin{align} r &= \dfrac{\sum_{i=1}^n (X_i - \bar X )(Y_i - \bar Y)}{\sqrt{\sum_{i=1}^n (X_i - \bar X )^2 \sum_{i=1}^n(Y_i - \bar Y)^2}} \\ & = \dfrac{S_{xy}}{\sqrt{S_{xx} S_{yy}}} \end{align} is also the maximum-likelihood estimator of $\rho$, the correlation coefficient.
However I don't know how to prove it. For bivariate distribution $(X, Y)$, we have pdf $$f(x,y) = \frac{1}{2\pi \sqrt{1-\rho^2} \sigma_x \sigma_y} \exp(-\frac{1}{2(1-\rho^2)}[(\frac{x-\mu_x}{\sigma_x})^2 - 2 \rho (\frac{x-\mu_x}{\sigma_x})(\frac{y-\mu_y}{\sigma_y})+(\frac{y-\mu_y}{\sigma_y})^2]).$$
The log-likelihood would yield $$\ln L = -n \ln (2\pi \sqrt{1-\rho^2} \sigma_x \sigma_y) + \sum -\frac{1}{2(1-\rho^2)}[...]$$ and taking the partial derivative w.r.t. $\rho$ and setting it to be zero doesn't seem to take me anywhere.
I have found some good resources for special case when $\mu_x = \mu_y = 0$ and $\mathbb{E}[X_i^2] = \mathbb{E}[Y_i^2] = 1$, but I am wondering how to approach the general problem, or if someone can direct me to a great resource I would really appreciate it.