  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\"odinger
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 normal-ordered two-body (NO2B) approximation that
captures dominant effects of three-nucleon forces while effectively working
with two-nucleon 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\"odinger 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-body formalisms in which U(1) global-gauge
symmetry associated with particle number conservation is broken (and
potentially restored), a particle-number-conserving normal-ordered k-body
(PNOkB) approximation of an arbitrary N-body operator is designed on the basis
of Bogoliubov reference states. A numerical test based on particle-number
projected Hartree-Fock-Bogoliubov calculations permits to check the
particle-number conserving/violating character of a given approximation to a
particle-number conserving operator. Using the presently proposed PNOkB
approximation, ab initio calculations based on symmetry-breaking and restored
formalisms can be safely performed.
