  We present a novel approach to generate higher-order initial conditions (ICs)
for cosmological simulations that take into account the distinct evolution of
baryons and dark matter. We focus on the numerical implementation and the
validation of its performance, based on both collisionless N-body simulations
and full hydrodynamic Eulerian and Lagrangian simulations. We improve in
various ways over previous approaches that were limited to first-order
Lagrangian perturbation theory (LPT). Specifically, we (1) generalize nth-order
LPT to multi-fluid systems, allowing 2LPT or 3LPT ICs for two-fluid
simulations, (2) employ a novel propagator perturbation theory to set up ICs
for Eulerian codes that are fully consistent with 1LPT or 2LPT, (3) demonstrate
that our ICs resolve previous problems of two-fluid simulations by using
variations in particle masses that eliminate spurious deviations from expected
perturbative results, (4) show that the improvements achieved by going to
higher-order PT are comparable to those seen for single-fluid ICs, and (5)
demonstrate the excellent (i.e., few per cent level) agreement between Eulerian
and Lagrangian simulations, once high-quality initial conditions are used. The
rigorous development of the underlying perturbation theory is presented in a
companion paper. All presented algorithms are implemented in the Monofonic
Music-2 package that we make publicly available.
