Turbulent motions are essential to the mixing of entrained fluids and are also capable of amplifying weak initial magnetic fields by small-scale dynamo action . Here we perform a systematic study of turbulent mixing in magnetized media , using three-dimensional magnetohydrodynamic simulations that include a scalar concentration field . We focus on how mixing depends on the magnetic Prandtl number , Pm , from 1 to 4 and the Mach number , { \mathcal { M } } , from 0.3 to 2.4 . For all subsonic flows , we find that the velocity power spectrum has a k ^ { -5 / 3 } slope in the early , kinematic phase , but steepens due to magnetic back reactions as the field saturates . The scalar power spectrum , on the other hand , flattens compared to k ^ { -5 / 3 } at late times , consistent with the Obukohov-Corrsin picture of mixing as a cascade process . At higher Mach numbers , the velocity power spectrum also steepens due to the presence of shocks , and the scalar power spectrum again flattens accordingly . Scalar structures are more intermittent than velocity structures in subsonic turbulence while for supersonic turbulence , velocity structures appear more intermittent than the scalars only in the kinematic phase . Independent of the Mach number of the flow , scalar structures are arranged in sheets in both the kinematic and saturated phases of the magnetic field evolution . For subsonic turbulence , scalar dissipation is hindered in the strong magnetic field regions , probably due to Lorentz forces suppressing the buildup of scalar gradients , while for supersonic turbulence , scalar dissipation increases monotonically with increasing magnetic field strength . At all Mach numbers , mixing is significantly slowed by the presence of dynamically-important small-scale magnetic fields , implying that mixing in the interstellar medium and in galaxy clusters is less efficient than modeled in hydrodynamic simulations .