As I mentioned in my comment and as Royi pointed out, you should check your formula on the CRLB. Additionally, since you are using additive white Gaussian noise, the method of maximum likelihood (ML) converges to the non-linear least squares (NLS) solution. Furthermore, since you appear to be estimating a single frequency, I will also show that the NLS solution for a single frequency is the peak of the periodogram.
So, firstly defining the NLS problem. The NLS problem is defined as
\begin{equation}
f(\omega,\alpha,\phi) = \sum_{n=0}^{N-1}\left\lvert y[n]-\sum_{m=0}^{M-1}\alpha_{m}e^{j\left(\omega_{m}n+\phi_{m}\right)}\right\rvert^{2}
\end{equation}
We assume our signal is of the form
\begin{equation}
y = \psi(\gamma) + e
\end{equation}
where $e$ is white Gaussian noise, $\gamma$ is some unknown parameter, and $\psi(\gamma)$ is a deterministic function of $\gamma$. The pdf is then
\begin{equation}
p\left(y\,|\,\psi(\gamma),\sigma^{2}\right) = \frac{1}{(2\pi)^{N}\left(\sigma^{2}\right)^{N}}e^{-\frac{\left\lvert y-\psi(\gamma)\right\rvert^{2}}{2\sigma^{2}}}
\end{equation}
Taking the log we get
\begin{align}
\log\left(p\left(y\,|\,\psi(\gamma),\sigma^{2}\right)\right) &= -N\log(2\pi) - N\log(\sigma^{2})-\frac{\left\lvert y-\psi(\gamma)\right\rvert^{2}}{2\sigma^{2}} \\ -\log\left(p\left(y\,|\,\psi(\gamma),\sigma^{2}\right)\right) &= N\log(2\pi) + N\log(\sigma^{2})+\frac{\left\lvert y-\psi(\gamma)\right\rvert^{2}}{2\sigma^{2}}
\end{align}
We can choose to minimize the negative log-likelihood to maximize the log-likelihood. Taking the derivative with respect to $\sigma^{2}$ gives
\begin{equation}
\frac{\partial}{\partial\sigma^{2}}\left[-\log\left(p\left(y\,|\,\psi(\gamma),\sigma^{2}\right)\right)\right] = \frac{N}{\sigma^{2}} - \frac{\left\lvert y-\psi(\gamma)\right\rvert^{2}}{2\sigma^{4}} = 0
\end{equation}
This gives
\begin{equation}
\hat{\sigma}^{2} = \frac{\left\lvert y-\psi(\gamma)\right\rvert^{2}}{N}
\end{equation}
Plugging this back in and gives
\begin{equation}
-\log\left(p\left(y\,|\,\psi(\gamma),\sigma^{2}\right)\right) = N\log(2\pi) + N\log\left(\frac{\left\lvert y-\psi(\gamma)\right\rvert^{2}}{N}\right)+N
\end{equation}
Since we wish to minimize the negative log-likelihood, and $\log()$ is a monotonically increasing function, we find that
\begin{equation}
\hat{\gamma} = \text{arg }\min_{\gamma}\left\lvert y-\psi(\gamma)\right\rvert^{2}
\end{equation}
which is precisely the NLS problem!
Now, let's reformulate the NLS problem into vector form. Let the following be true
\begin{align}
\beta_{m} &= \alpha_{m}e^{j\phi_{m}} \\ \beta &= \begin{bmatrix}\beta_{0} & \cdots & \beta_{M-1} \end{bmatrix} \\ \underline{y} &= \begin{bmatrix} y[0] & \cdots & y[N-1] \end{bmatrix} \\ \mathbf{B} &= \begin{bmatrix}1 & \cdots & 1 \\ \vdots & & \vdots \\ e^{j\omega_{0}(N-1)} & \cdots & e^{j\omega_{M-1}(N-1)} \end{bmatrix}
\end{align}
where $\mathbf{B}$ is a Vandermonde matrix. It can be shown 1 under certain assumptions that the NLS estimates of the frequencies are
\begin{equation}
\hat{\omega} = \text{arg }\max_{\omega}\left[\underline{y}^{H}\mathbf{B}\left(\mathbf{B}^{H}\mathbf{B}\right)^{-1}\mathbf{B}^{H}\underline{y}\right]
\end{equation}
In the case that there's only one true frequency, the matrix $\mathbf{B}$ becomes a steering vector to that frequency
\begin{equation}
\underline{a}(\omega) = \begin{bmatrix} 1 & e^{j\omega} & \cdots & e^{j\omega(N-1)} \end{bmatrix}^{T}
\end{equation}
The equation for $\hat{\omega}$ then becomes
\begin{align}
\hat{\omega} &= \text{arg }\max_{\omega} = \frac{\underline{y}^{H}\underline{a}(\omega)\underline{a}^{H}(\omega)\underline{y}}{\underline{a}^{H}(\omega)\underline{a}(\omega)} \\ &= \text{arg }\max_{\omega}\frac{\left\lvert\underline{a}^{H}(\omega)\underline{y}\right\rvert^{2}}{N}
\end{align}
which is the periodogram at frequency $\omega$ (see Computing modern spectral estimation techniques with FFTs for more details on why this is the case)!
There is a very in depth analysis in 1 on the bias and variance of the periodogram. Since the periodogram is asymptotically unbiased, the MSE will be the variance, which will be quite large (see more on that in Why is the sample spectrum considered inconsistent?). If you are interested in multiple frequencies, there is further analysis in 1 to show that the periodogram gives $\mathcal{O}\left(\frac{1}{N}\right)$ estimates of the NLS frequency estimates provided the finite sample periodogram has sufficient resolution such that all frequencies are spaced at least one DFT bin apart.
All of this is to show that you should absolutely expect the maximum likelihood estimate of the true frequencies, given a sinusoid in additive noise model, to be above the CRLB.
1 Stoica, P., & Moses, R. L. (2005). Spectral analysis of signals (Vol. 452, pp. 25-26). Upper Saddle River, NJ: Pearson Prentice Hall.
