  Multiplex networks describe systems whose interactions can be of different
nature, and are fundamental to understand complexity of networks beyond the
framework of simple graphs. Recently it has been pointed out that restricting
the attention to pairwise interactions is also a limitation, as the vast
majority of complex systems include higher-order interactions that strongly
affect their dynamics. Here, we propose hyper-diffusion on multiplex networks,
a dynamical process in which diffusion on each single layer is coupled with the
diffusion in other layers thanks to the presence of higher-order interactions
occurring when there exists link overlap. We show that hyper-diffusion on a
duplex network can be described by the Hyper-Laplacian in which the strength of
four-body interactions among every set of four replica nodes connected by a
multilink $(1,1)$ can be tuned by a parameter $\delta_{11}\ge 0$. The
Hyper-Laplacian reduces to the standard lower Laplacian, capturing pairwise
interactions at the two layers, when $\delta_{11}=0$. By combining tools of
spectral graph theory, applied topology and network science we provide a
general understanding of hyper-diffusion on duplex networks when
$\delta_{11}>0$, including theoretical bounds on the Fiedler and the largest
eigenvalue of Hyper-Laplacians and the asymptotic expansion of their spectrum
for $\delta_{11}\ll1$ and $\delta_{11}\gg1$. Although hyper-diffusion on
multiplex networks does not imply a direct ``transfer of mass" among the layers
(i.e. the average state of replica nodes in each layer is a conserved quantity
of the dynamics), we find that the dynamics of the two layers is coupled as the
relaxation to the steady state becomes synchronous when higher-order
interactions are taken into account and the Fiedler eigenvalue of the
Hyper-Laplacian is not localized in a single layer of the duplex network.
