  We present an efficient Monte Carlo framework for perturbative calculations
of infinite nuclear matter based on chiral two-, three-, and four-nucleon
interactions. The method enables the incorporation of all many-body
contributions in a straightforward and transparent way, and makes it possible
to extract systematic uncertainty estimates by performing order-by-order
calculations in the chiral expansion as well as the many-body expansion. The
versatility of this new framework is demonstrated by applying it to chiral
low-momentum interactions, exhibiting a very good many-body convergence up to
fourth order. Following these benchmarks, we explore new chiral interactions up
to next-to-next-to-next-to-leading order (N$^3$LO). Remarkably, simultaneous
fits to the triton and to saturation properties can be achieved, while all
three-nucleon low-energy couplings remain natural. The theoretical
uncertainties of nuclear matter are significantly reduced when going from
next-to-next-to-leading order to N$^3$LO.
