  In many astrophysical environments, mixing of heavy elements occurs in the
presence of a supersonic turbulent velocity field. Here we carry out the first
systematic numerical study of such passive scalar mixing in isothermal
supersonic turbulence. Our simulations show that the ratio of the scalar mixing
timescale, $\tau_{\rm c}$, to the flow dynamical time, $\tau_{\rm dyn}$
(defined as the flow driving scale divided by the rms velocity), increases with
the Mach number, $M$, for $M \lsim3$, and becomes essentially constant for $M
\gsim3.$ This trend suggests that compressible modes are less efficient in
enhancing mixing than solenoidal modes. However, since the majority of kinetic
energy is contained in solenoidal modes at all Mach numbers, the overall change
in $\tau_{\rm c}/\tau_{\rm dyn}$ is less than 20\% over the range $1 \lsim M
\lsim 6$. At all Mach numbers, if pollutants are injected at around the flow
driving scale, $\tau_{\rm c}$ is close to $\tau_{\rm dyn}.$ This suggests that
scalar mixing is driven by a cascade process similar to that of the velocity
field. The dependence of $\tau_{\rm c}$ on the length scale at which pollutants
are injected into flow is also consistent with this cascade picture. Similar
behavior is found for the variance decay timescales for scalars without
continuing sources. Extension of the scalar cascade picture to the supersonic
regime predicts a relation between the scaling exponents of the velocity and
the scalar structure functions, with the scalar structure function becoming
flatter as the velocity scaling steepens with Mach number. Our measurements of
the volume-weighted velocity and scalar structure functions confirm this
relation for $M\lsim 2,$ but show discrepancies at $M \gsim 3$.
