$\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=\pm1)=1/2$$\PP(X=1)=\PP(X=-1)=1/2$.