Accurate $N$-body simulations with local Primordial non-Gaussianities: initial conditions and aliasing

New generation galaxy surveys targeting constraints on local primordial non-Gaussianity (PNG) demand $N$-body simulations that accurately reproduce its effects. In this work, we explore various prescriptions for the initial conditions of simulations with PNG, aiming to optimise accuracy and minimise numerical errors, particularly due to aliasing. We have used $186$ runs that vary the starting redshift, LPT order, and non-Gaussianities ($f^{\rm local}_{\rm NL}$ and $g^{\rm local}_{\rm NL}$). Starting with $3$LPT at a redshift as low as $z_{\rm ini}\simeq 11.5$ reproduces to $<1 \%$ the power spectrum, bispectrum and halo mass function of a high-resolution reference simulation. The aliasing induced by the PNG terms in the power spectrum produces a $\leq 3 \%$ excess small-scale power at the initial conditions, dropping below $0.1\%$ by $z=0$. State-of-the-art initial condition generators show a sub-percent agreement. We show that initial conditions for simulations with PNG should be established at a lower redshift using higher-order LPT schemes. We also show that removing the PNG aliasing signal is unnecessary for current simulations. The methodology proposed here can accelerate the generation of simulations with PNG while enhancing their accuracy.