What does it mean by "maximum likelihood estimation (MLE) problem is unbounded"? I saw the following statement in this paper and wanted to understand its meaning:


$\DeclareMathOperator{\diag}{diag}$
$\DeclareMathOperator{\tr}{tr}$
A great question. Most standard multivariate analysis texts only treated the $N > m$ case to give the MLE of $\mu$ and $X$ (or more frequently, $\Sigma = X^{-1}$).  The $N \leq m$ case seems always overlooked (or merely qualitatively mentioned).
In short, "the MLE problem is unbounded" means that the objective function
\begin{align}
l(\mu, X) := -\log\det(X) + \frac{1}{N}\sum_{i = 1}^N(\hat{\xi}_i - \mu)^TX(\hat{\xi}_i - \mu), \quad \mu \in \mathbb{R}^m, X \in S_+^m 
\end{align}
does not have a fixed lower bound (which depends on $\hat{\xi}_1, \ldots, \hat{\xi}_N$ only) when $m \geq N$ (by contrast, when $N > m$, almost every multivariate analysis textbook will show you $l(\mu, X)$ is bounded below).  To show that $l(\mu, X)$ can be arbitrarily small, first introduce the following notations:
\begin{align}
& Z = \begin{bmatrix}
\hat{\xi}_1^T \\
\vdots \\
\hat{\xi}_N^T
\end{bmatrix} \in \mathbb{R}^{N \times m}, \;
e = \begin{bmatrix}
1 \\
\vdots \\
1
\end{bmatrix} \in \mathbb{R}^{N \times 1}, \;
\bar{\xi} = N^{-1}\sum_{i = 1}^N\hat{\xi}_i, \\
& A = \sum_{i = 1}^N(\hat{\xi}_i - \bar{\xi})(\hat{\xi}_i - \bar{\xi})^T = Z^T(I_{(N)} - N^{-1}ee^T)Z := Z^TPZ \in \mathbb{R}^{m \times m}.
\end{align}
It can be seen that $P$ is idempotent, whence $A = (PZ)^TPZ$ is positive semi-definite, and
\begin{align}
r:= \operatorname{rank}(A) = \operatorname{rank}(PZ) \leq \operatorname{rank}(P) =  N - 1 < m.
\end{align}
Therefore, the spectral decomposition of $A$ can be written as
\begin{align}
A = O\diag(\Lambda, 0_{(m - r)})O^T,
\end{align}
where $O$ is an order $m$ orthogonal matrix, $\Lambda = \diag(\lambda_1, \ldots, \lambda_r)$, $\lambda_i (1 \leq i \leq r)$ are positive eigenvalues of $A$.
Using these notations, $l(\mu, X)$ can be rewritten as:
\begin{align}
l(\mu, X) = -\log\det(X)+\frac{1}{N}\tr(XA) + (\bar{\xi} - \mu)^TX(\bar{\xi} - \mu).
\end{align}
Let $X^* = O\diag(\Lambda^{-1}, kI_{(m - r)})O^T$ (where $k > 0$ is to be determined), $\mu^* = \bar{\xi}$, it then follows that
\begin{align}
& l(\mu^*, X^*) = \sum_{i = 1}^r\log\lambda_i - (m - r)\log k+ \frac{r}{N},
\end{align}
whence $l(\mu^*, X^*) \to -\infty$ as $k \to \infty$, i.e., $l$ cannot be bounded from below.
If the above matrix operations are too complicated to understand, you can get some sense by considering the case $m = N = 1$ (univariate normal distribution with only one observation), for which the objective function is
\begin{align}
l(\mu, X) = -\log X + (\hat{\xi}_1 - \mu)^2X.
\end{align}
Therefore $l(\hat{\xi}_1, X) = -\log X \to -\infty$ as $X \to +\infty$, i.e., the objective function is unbounded.
