$\newcommand{\de}{\delta}
\newcommand{\De}{\Delta}
\newcommand{\ep}{\epsilon}
\newcommand{\ga}{\gamma}
\newcommand{\Ga}{\Gamma}
\newcommand{\la}{\lambda}
\newcommand{\Si}{\Sigma}
\newcommand{\thh}{\theta}
\newcommand{\R}{\mathbb{R}}
\newcommand{\F}{\mathcal{F}}
\newcommand{\E}{\operatorname{\mathsf E}}
\newcommand{\PP}{\operatorname{\mathsf P}}
\newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}}$
Your guess is correct. Indeed, it is well known (see e.g. Bertsimas--Popesku, page 781) that real numbers $m_0=1,m_1,\dots,m_{2\ell}$ are the moments of orders $0,1,\dots,2\ell$ of a real-valued random variable $X$ iff the matrix $M:=(m_{i+j})_{i,j=0}^\ell$ is nonnegative-definite, that is, iff all the principal minors of $M$ are $\ge0$; here $\ell$ is a natural number; in our case, $\ell=2$.
Also, given $\E X^4=1$, we have
\begin{equation}
\E(X-Y)^4=2-8m_3m_1+6m_2^2.
\end{equation}
Thus, the problem is a simple problem of real algebraic geometry, which can be solved algorithmically. Using the Mathematica command Maximize[], we get the result:
Added: Here is an elementary solution, without using Mathematica: Since $m_4=1$, the condition $\det M\ge0$ implies
\begin{equation}
2m_3m_1\ge m_2^2-m_4+(m_3^2+m_4m_1^4)/m_2\ge m_2^2-1,
\end{equation}
whence
\begin{equation}
\E(X-Y)^4=2-8m_3m_1+6m_2^2\le2-4(m_2^2-1)+6m_2^2=2m_2^2+6\le8,
\end{equation}
since $m_2^2\le m_4=1$. The equality in the inequality in question is attained only if $m_2^2=m_4=1$ and $2m_3m_1=m_2^2-1=0$, that is, only if $\PP(X=1)=\PP(X=-1)=1/2$.