  We present a model-independent investigation of the WMAP data with respect to
scale- dependent non-Gaussianities (NGs) by employing the method of constrained
randomization. For generating so-called surrogate maps a shuffling scheme is
applied to the Fourier phases of the original data, which allows to test for
the presence of higher order correlations (HOCs) on well-defined scales. Using
scaling indices as test statistics we find highly significant signatures for
non-Gaussianities when considering all scales. We test for NGs in the bands l =
[2,20], l = [20,60], l = [60,120] and l = [120,300]. We find highly significant
signatures for non-Gaussianities and ecliptic hemispherical asymmetries for l =
[2, 20]. We also obtain highly significant deviations from Gaussianity for the
band l = [120,300]. The result for the full l-range can be interpreted as a
superposition of the signatures found in the bands l = [2, 20] and l = [120,
300]. We find remarkably similar results when analyzing different ILC-like
maps. We perform a set of tests to investigate if the detected anomalies can be
explained by systematics. While no test can convincingly rule out the intrinsic
nature of the anomalies for the low l case, the ILC map making procedure and/or
residual noise in the maps can also lead to NGs at small scales. Our
investigations prove that there are phase correlations in the WMAP data of the
CMB. In the absence of an explanation in terms of Galactic foregrounds or known
systematic artefacts, the signatures at low l must so far be taken to be
cosmological at high significance. These findings strongly disagree with
predictions of isotropic cosmologies with single field slow roll inflation. The
task is now to elucidate the origin of the phase correlations and to understand
the physical processes leading to these scale-dependent non-Gaussianities - if
systematics as cause for them must be ruled out.
