  Probability density functions (PDFs) of scale-dependent energy fluctuations,
$P[\delta E(\ell)]$, are studied in high-resolution direct numerical
simulations of Navier-Stokes and incompressible magnetohydrodynamic (MHD)
turbulence. MHD flows with and without a strong mean magnetic field are
considered. For all three systems it is found that the PDFs of inertial range
energy fluctuations exhibit self-similarity and monoscaling in agreement with
recent solar-wind measurements [B. Hnat et al., Geophys. Res. Lett. 29(10),
86-1 (2002)]. Furthermore, the energy PDFs exhibit similarity over all scales
of the turbulent system showing no substantial qualitative change of shape as
the scale of the fluctuations varies. This is in contrast to the well-known
behavior of PDFs of turbulent velocity fluctuations. In all three cases under
consideration the $P[\delta E(\ell)]$ resemble L\'evy-type gamma distributions
$\sim \Delta^{-1}\exp(-|\delta E|/\Delta)|\delta E|^{-\gamma}$ The observed
gamma distributions exhibit a scale-dependent width $\Delta(\ell)$ and a
system-dependent $\gamma$. The monoscaling property reflects the inertial-range
scaling of the Els\"asser-field fluctuations due to lacking Galilei invariance
of $\delta E$. The appearance of L\'evy distributions is made plausible by a
simple model of energy transfer.
