Spectral peak location estimation using complex DFT
I admit that I failed my first attempt.
The text suggests that we would take $g(m) = am^2 + b m + c$, interpolating $g(k) = X(k+l)$, for $k \in \{-1, 0, 1\}$, interpolating this would give us polynomial parameterized in $X[k], X[k-1], X[k+1]$, derive, solve and it gives the optimum value at $m = k+\delta$. Not that easy, and it gave me three roots, none of them seemed to be the desired value :(
The comment on the paper is not helpful but they give us a reference [4], from the one of the authors, and there it gives more details. There you will find the equation (3.7) and they give a concise derivation on the Appendix A.
Derivation
Here I will fill some details of the proof.
You express the interpolator polynomial centered at $k + \delta$, so delta will maximize it by construction (more clever the attempt described above).
$$g(m) = a (m - k - \delta)^2 + b$$
Then $g(m)$, determine delta so that interpolates the three samples. Since the quadratic polynomial interpolating three points is unique, we know that $\delta$ here is unique as well.
$$\begin{eqnarray}
X[k-1] &amp;=&amp; g(k-1) &amp;=&amp; b+((k-1) - k - \delta)^2 &amp;=&amp; b + a\cdot(\delta^2+2\delta + 1) \\
X[k  ] &amp;=&amp; g(k  ) &amp;=&amp; b+((k  ) - k - \delta)^2 &amp;=&amp; b + a\cdot(\delta^2) \\
X[k+1] &amp;=&amp; g(k+1) &amp;=&amp; b+((k+1) - k - \delta)^2 &amp;=&amp; b + a\cdot(\delta^2-2\delta + 1)
\end{eqnarray}$$
Start looking at the coefficient of the rightmost terms, and note that
$$2 X[k] - X[k-1] - X[k+1] = 2 a$$
Also note that
$$ X[k-1] - X[k+1] = 4 a \delta$$
$$\delta = \frac{X[k-1] - X[k+1]}{4a} = \frac{X[k-1] - X[k+1]}{2(2X[k] - X[k-1] - X[k+1])}$$
With the sentence

A simple adaptation provides an estimator suitable for estimating tone frequencies from complex DFT outputs.
They introduce the formula we see in the first paper you mentioned.

The adaptation they mention is because $\delta$ is complex, so he takes the real value. Maybe a justification for that would be that $\min_{f} |f - \delta|$ is $Re[\delta]$.