The reach of ab initio many-body theories is rapidly extending over the nuclear chart. However, dealing fully with three-nucleon, possibly four-nucleon, interactions makes the solving of the A-body Schrödinger equation particularly cumbersome, if not impossible beyond a certain nuclear mass. Consequently, ab initio calculations of mid-mass nuclei are typically performed on the basis of the so-called normal-ordered two-body (NO2B) approximation that captures dominant effects of three-nucleon forces while effectively working with twonucleon operators. A powerful idea currently employed to extend ab initio calculations to open-shell nuclei consists of expanding the exact solution of the A-body Schrödinger equation while authorizing the approximate solution to break symmetries of the Hamiltonian. In this context, operators are normal ordered with respect to a symmetry-breaking reference state such that proceeding to a naive truncation may lead to symmetry-breaking approximate operators. The purpose of the present work is to design a normal-ordering approximation of operators that is consistent with the symmetries of the Hamiltonian while working in the context of symmetry broken (and potentially restored) methods. Focusing on many-J. Ripoche CEA, DAM,