ABSTRACT We numerically investigate the non-Gaussianities in the late-time cosmological density field in Fourier space. We explore various statistics, including the two- and three-point probability distribution function (PDF) of phase and modulus, and their two- and three-point correlation function. Significant non-Gaussianity is observed for certain configurations. Comparing the measurement from simulation with the theoretical expansion prediction, we find that for (600 Mpc h-1)3 volume, the $\mathcal {O}(V^{-1/2})$ order term alone is sufficiently accurate to describe all the measured non-Gaussianities in not only the PDFs, but also the correlations. We also numerically find that the phase-modulus cross-correlation contributes $\sim 50{{\ \rm per\ cent}}$ to the bispectrum, further verifying the accuracy of the $\mathcal {O}(V^{-1/2})$ order prediction. This work demonstrates that the non-Gaussianity of cosmic density field is simpler in Fourier space, and may facilitate the data analysis in the era of precision cosmology.