$\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, if 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:

[![enter image description here][1]][1]


  [1]: https://i.sstatic.net/LuyF3.png