AbstractMultiplex 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 (a multiplex network with two layers) can be described by the hyper-Laplacian in which the strength of four-body interactions among every set of four replica nodes connected in both layers can be tuned by a parameterδ11⩾ 0. The hyper-Laplacian reduces to the standard lower Laplacian, capturing pairwise interactions at the two layers, whenδ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δ11> 0, including theoretical bounds on the Fiedler and the largest eigenvalue of hyper-Laplacians and the asymptotic expansion of their spectrum forδ11≪ 1 andδ11≫ 1. 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.