Not yet reviewed by Pith; the record is open.

The probability density function (PDF) of the logarithmic density contrast, $s=\ln (\rho/\rho_0)$, with gas density $\rho$ and mean density $\rho_0$, for hydrodynamical supersonic turbulence is well-known to have significant non-Gaussian (intermittent) features that monotonically increase with the turbulent Mach number, $\mathcal{M}$. By studying the mass- and volume-weighted $s$-PDF for an ensemble of 36 sub-to-trans-Alf\'venic mean-field, supersonic, isothermal turbulence simulations with different modes of driving, relevant to molecular gas in the cool interstellar medium, we show that a more intricate picture emerges for the non-Gaussian nature of $s$. Using four independent measures of the non-Gaussian components, we find hydrodynamical-like structure in the highly magnetised plasma for $\mathcal{M} \lesssim 4$. However, for $\mathcal{M} \gtrsim 4$, the non-Gaussian signatures disappear, leaving approximately Gaussian $s$-statistics -- exactly the opposite of hydrodynamical turbulence in the high-$\mathcal{M}$ limit. We also find that the non-Gaussian components of the PDF increase monotonically with more compressive driving modes. To understand the $\mathcal{M} \lesssim 4$ non-Gaussian features we use one-dimensional (1D) pencil beams to explore the dynamics along and across the large-scale magnetic field, $\mathbf{B}_0$. We discuss kinetic, density and magnetic field fluctuations from the pencil beams, and identify physical sources of non-Gaussian components to the PDF as single, strong shocks coupled to fast magnetosonic compressions that form along $\mathbf{B}_0$. We discuss the Gaussianisation of the $\mathcal{M} \gtrsim 4$ $s$-fields through the lens of two phenomenologies: the self-similarity of the $s$-field and homogenisation of the dynamical timescales between the over- and under-dense regions in the compressible gas.