[proofplan]
We first replace the real cutoff $X$ by the integer cutoff $P=\lfloor X\rfloor$. The core estimate is precisely Hua's mean value bound for kth-power Weyl sums with integer length $P$. Applying that estimate to the map $\alpha\mapsto f_k(\alpha;X)$ gives the asserted inequality in $P$, and the comparison $P\le X$ converts it to the stated form in $X$.
[/proofplan]
[step:Apply the integer-length Hua mean value bound]
Let $P:=\lfloor X\rfloor$. If $k=1$, then the only admissible integer is $j=1$. Expanding the square and using the orthogonality relation $\int_0^1 e((n-m)\alpha)\,d\mathcal L^1(\alpha)=1$ when $n=m$ and $0$ when $n\ne m$, we obtain
\begin{align*}
\int_0^1 |f_1(\alpha;X)|^2\,d\mathcal L^1(\alpha)=P.
\end{align*}
Since $P\le X\le X^{1+\varepsilon}$, the desired estimate follows for $k=1$. Hence, for the rest of the proof, assume $k\ge 2$.
If $1\le X<2$, then $P=1$ and $f_k(\alpha;X)=e(\alpha)$ for every $\alpha\in\mathbb R$, so
\begin{align*}
\int_0^1 |f_k(\alpha;X)|^{2^j}\,d\mathcal L^1(\alpha)=1\le X^{2^j-j+\varepsilon}.
\end{align*}
Assume now that $P\ge 2$. Define the integer-length Weyl sum $g_P:[0,1]\to\mathbb C$ by
\begin{align*}
g_P(\alpha):=\sum_{1\le n\le P}e(\alpha n^k).
\end{align*}
By the definition of $f_k(\cdot;X)$, one has $g_P(\alpha)=f_k(\alpha;X)$ for every $\alpha\in[0,1]$.
The [[Hua Mean Value Bound](/theorems/9075)][citetheorem:9075] applies with this integer $P$, the same integer $k$, and the chosen integer $j$ with $1\le j\le k$. Its hypotheses are satisfied because $k\ge 2$, $P\ge 1$, and $\varepsilon>0$. Therefore
\begin{align*}
\int_0^1 |f_k(\alpha;X)|^{2^j}\,d\mathcal L^1(\alpha)=\int_0^1 |g_P(\alpha)|^{2^j}\,d\mathcal L^1(\alpha)\lesssim_{j,k,\varepsilon}P^{2^j-j+\varepsilon}.
\end{align*}
[guided]
The only difference between the theorem statement and the standard integer-length form is the real cutoff $X$. We remove that difference by setting $P:=\lfloor X\rfloor$ and introducing the map $g_P:[0,1]\to\mathbb C$ by
\begin{align*}
g_P(\alpha):=\sum_{1\le n\le P}e(\alpha n^k).
\end{align*}
The definition of $f_k(\cdot;X)$ uses exactly the same integer $P$, so $g_P(\alpha)=f_k(\alpha;X)$ for every $\alpha\in[0,1]$.
We now invoke the [Hua Mean Value Bound][citetheorem:9075]. That result requires an integer $k\ge 2$, an integer length $P\ge 1$, an integer $j$ with $1\le j\le k$, and a number $\varepsilon>0$. These are precisely the data in scope after the separate $k=1$ case: $k\ge 2$, $P=\lfloor X\rfloor\ge 1$ in the nontrivial case, $j$ is fixed with $1\le j\le k$, and $\varepsilon>0$ is arbitrary. The theorem gives
\begin{align*}
\int_0^1 |g_P(\alpha)|^{2^j}\,d\mathcal L^1(\alpha)\lesssim_{j,k,\varepsilon}P^{2^j-j+\varepsilon}.
\end{align*}
Replacing $g_P$ by the equal function $f_k(\cdot;X)$ on $[0,1]$ gives
\begin{align*}
\int_0^1 |f_k(\alpha;X)|^{2^j}\,d\mathcal L^1(\alpha)\lesssim_{j,k,\varepsilon}P^{2^j-j+\varepsilon}.
\end{align*}
[/guided]
[/step]
[step:Pass from the integer cutoff to the real cutoff]
Since $P=\lfloor X\rfloor\le X$ and $X\ge 1$, the exponent $2^j-j+\varepsilon$ is positive. Hence
\begin{align*}
P^{2^j-j+\varepsilon}\le X^{2^j-j+\varepsilon}.
\end{align*}
Combining this comparison with the previous step yields
\begin{align*}
\int_0^1 |f_k(\alpha;X)|^{2^j}\,d\mathcal L^1(\alpha)\lesssim_{j,k,\varepsilon}X^{2^j-j+\varepsilon}.
\end{align*}
This proves the estimate for every integer $j$ with $1\le j\le k$. Taking $j=k$ gives the endpoint estimate.
[guided]
The preceding step gives the estimate with the integer cutoff $P=\lfloor X\rfloor$:
\begin{align*}
\int_0^1 |f_k(\alpha;X)|^{2^j}\,d\mathcal L^1(\alpha)\lesssim_{j,k,\varepsilon}P^{2^j-j+\varepsilon}.
\end{align*}
To convert this to the stated real cutoff, we use the defining property of the floor function: $P\le X$. The exponent is positive because $2^j-j\ge 1$ for every integer $j\ge 1$ and $\varepsilon>0$. Therefore monotonicity of $t\mapsto t^{2^j-j+\varepsilon}$ on $[0,\infty)$ gives
\begin{align*}
P^{2^j-j+\varepsilon}\le X^{2^j-j+\varepsilon}.
\end{align*}
Substituting this into the previous estimate yields
\begin{align*}
\int_0^1 |f_k(\alpha;X)|^{2^j}\,d\mathcal L^1(\alpha)\lesssim_{j,k,\varepsilon}X^{2^j-j+\varepsilon}.
\end{align*}
Since $j$ was arbitrary with $1\le j\le k$, this proves the full family of estimates. Choosing $j=k$ gives the endpoint estimate.
[/guided]
[/step]