  The Many-Body Expansion (MBE) is a useful tool to simulate condensed phase
chemical systems, often avoiding the steep computational cost of usual
electronic structure methods. However, it often requires higher than 2-body
terms to achieve quantitative accuracy. In this work, we propose the Polarized
MBE (PolBE) method where each MBE energy contribution is treated as an
embedding problem. In each energy term, a smaller fragment is embedded into a
larger, polarized environment and only a small region is treated at the
high-level of theory using embedded mean-field theory. The role of polarized
environment was found to be crucial in providing quantitative accuracy at the
2-body level. PolBE accurately predicts non-covalent interaction energies for a
number of systems, including CO$_2$, water, and hydrated ion clusters, with a
variety of interaction mechanisms, from weak dispersion to strong
electrostatics considered in this work. We further demonstrate that the PolBE
interaction energy is predominantly pairwise unlike the usual vacuum MBE which
requires higher-order terms to achieve similar accuracy. We numerically show
that PolBE often performs better than other widely used embedded MBE methods
such as the electrostatically embedded MBE. Owing to the lack of expensive
diagonalization of Fock matrices and its embarrassingly parallel nature, PolBE
is a promising way to access condensed phase systems with hybrid density
functionals that are difficult to treat with currently available methods.
