Using whuber's method, we'll reject if $\ell(\theta_0; \vec x) - \ell(\theta_{MLE}; \vec x) \le k$ for some constant $k$, where $\ell(\theta; \vec x) = \ln L(\theta ; \vec x) = \ln \left(\theta^n \prod_{i=1}^n x_i^\theta \right) = n \ln \theta + \theta \sum_{i=1}^n \ln x_i$.
We see that $\ell(\theta_0; \vec x) - \ell(\theta_{MLE}; \vec x) = n (\ln \theta_0 - \ln \theta_{MLE}) + (\theta_0 - \theta_{MLE}) \sum_{i=1}^n \ln x_i$.
Using some generic constant $k$ at each step, we reject if
\begin{align*}
n (\ln \theta_0 - \ln \theta_{MLE}) + (\theta_0 - \theta_{MLE}) \sum_{i=1}^n \ln x_i \le k \\
-n \ln \theta_{MLE} + (\theta_0 - \theta_{MLE}) \left( - \frac{n}{\theta_{MLE}} \right) \le k \\
-n \ln \theta_{MLE} - \frac{n\theta_0}{\theta_{MLE}} \le k \\
\ln \theta_{MLE} + \frac{\theta_0}{\theta_{MLE}} \ge k
\end{align*}
Now, consider the function $f(z) = \ln(z) + \frac{\theta_0}{z}$. This has a minimum at $z = \theta_0$ and is concave up and $\lim \limits_{z \to 0^+} f(z) = \lim \limits_{z \to \infty} f(z) = \infty$.
So $f(z) \ge k$ means $z$ is either sufficiently small or sufficiently large.
That is, $ - \frac{n}{\sum_{i=1}^n \ln(x_i)}$ is either sufficiently small or sufficiently large; i.e. $- \sum_{i=1}^n \ln x_i \le c_1 $ or $ - \sum_{i=1}^n \ln x_i \ge c_2$. We want to choose $c_1$ and $c_2$ so that the test is size $\alpha$ under $H_0$.
Note that under $H_0$, $- \sum_{i=1}^n \ln X_i \sim \Gamma(n,\theta_0) \implies -2\theta_0 \sum_{i=1}^n \ln X_i \sim \Gamma(n,\frac{1}{2} ) = \chi^2(2n).$
So one LRT of size $\alpha$ is to reject $H_0$ in favor of $H_1$ if $- 2 \theta_0 \sum_{i=1}^n \ln X_i \le \chi^2_{1-\alpha/2}(2n)$ or if $- 2 \theta_0\sum_{i=1}^n \ln X_i \ge \chi^2_{\alpha/2}(2n)$.
The problem is in the matrix differentiation. As the covariance matrix is symmetric, we have
$
\frac{\partial l}{\partial \Sigma}=-\Sigma^{-1}+\frac{diag(\Sigma^{-1})}{2}+\Sigma^{-1}(x-\mu)(x-\mu)'\Sigma^{-1}-\frac{diag(\Sigma^{-1}(x-\mu)(x-\mu)'\Sigma^{-1})}{2}
$
where $l$ is the log-likelihood function.
Best Answer
If one includes the notational dependency on $n$: $$ \begin{align*} \ell_n\left(\theta\right) & = \ell_n\left(\widehat{\theta}_n\right)+\frac{\partial\ell_n\left(\theta\right)}{\partial\theta}\Bigr|_{\theta=\widehat{\theta}_n}\left(\theta-\widehat{\theta}_n\right)+ o_n(|\widehat{\theta}_n - \theta|^2) \\ \end{align*} $$ we see that the puzzling point is the $n$-dependency of the $o$.
A rigorous way to get an approximate with a $n$-independent $o$: $$ \begin{align*} \ell_n\left(\theta\right) & = \ell_n\left(\widehat{\theta}_n\right)+\frac{\partial\ell_n\left(\theta\right)}{\partial\theta}\Bigr|_{\theta=\widehat{\theta}_n}\left(\theta-\widehat{\theta}_n\right)+ o(|\widehat{\theta}_n - \theta|^2) \\ \end{align*} $$ is Taylor-Lagrange's inequality: if you are able to majorate $\ell_n'' \leq M$ uniformly in $n$ (on an appropriate interval) then you get the uniform $o$ by Taylor-Lagrange's inequality.