Generating non-Gaussian maps with a given power spectrum and bispectrum

We propose two methods for generating non-Gaussian maps with fixed power spectrum and bispectrum. The first makes use of a recently proposed rigorous, non-perturbative, Bayesian framework for generating non-Gaussian distributions. The second uses a simple superposition of Gaussian distributions. The former is best suited for generating mildly non-Gaussian maps, and we discuss in detail the limitations of this method. The latter is better suited for the opposite situation, i.e. generating strongly non-Gaussian maps. The ensembles produced are isotropic and the power spectrum can be jointly fixed; however we cannot set to zero all other higher order cumulants (an unavoidable mathematical obstruction). We briefly quantify the leakage into higher order moments present in our method. We finally present an implementation of our code within the HEALPIX package