  Suitable discretizations through tensor product formulas of popular
multidimensional operators (diffusion or diffusion--advection, for instance)
lead to matrices with $d$-dimensional Kronecker sum structure. For evolutionary
Partial Differential Equations containing such operators and integrated in time
with exponential integrators, it is then of paramount importance to efficiently
approximate the actions of $\varphi$-functions of the arising matrices. In this
work, we show how to produce directional split approximations of third order
with respect to the time step size. They conveniently employ tensor-matrix
products (the so-called $\mu$-mode product and related Tucker operator,
realized in practice with high performance level 3 BLAS), and allow for the
effective usage of exponential Runge--Kutta integrators up to order three. The
technique can also be efficiently implemented on modern computer hardware such
as Graphic Processing Units. The approach has been successfully tested against
state-of-the-art techniques on two well-known physical models that lead to
Turing patterns, namely the 2D Schnakenberg and the 3D FitzHugh--Nagumo
systems, on different hardware and software architectures.
