  We present a novel implementation of an extremum preserving anisotropic
diffusion solver for thermal conduction on the unstructured moving Voronoi mesh
of the AREPO code. The method relies on splitting the one-sided facet fluxes
into normal and oblique components, with the oblique fluxes being limited such
that the total flux is both locally conservative and extremum preserving. The
approach makes use of harmonic averaging points and a simple, robust
interpolation scheme that works well for strong heterogeneous and anisotropic
diffusion problems. Moreover, the required discretisation stencil is small.
Efficient fully implicit and semi-implicit time integration schemes are also
implemented. We perform several numerical tests that evaluate the stability and
accuracy of the scheme, including applications such as point explosions with
heat conduction and calculations of convective instabilities in conducting
plasmas. The new implementation is suitable for studying important
astrophysical phenomena, such as the conductive heat transport in galaxy
clusters, the evolution of supernova remnants, or the distribution of heat from
blackhole-driven jets into the intracluster medium.
