Comparing many-body approaches against the helium atom exact solution

Answers to the Editor 2nd point:
We are all researchers, and researchers are exploring the world amazed like children, and like children we must do it by playing (a fact that research financing agencies, especially European, do not understand, unfortunately). And models are our preferred toys. We are absolutely convinced about the importance of studying models: only by playing with models we can grasp the full complexity of, e.g., many-body theory. However, also like children, at a certain point we must take awareness of the limit where the game is over and the reality starts: the confusion of the two planes is not at all good (for both researchers and children). So far as a model could be a beautiful toy to be reluctantly given up, we must be aware of its limits as portrait of the reality. I find "forced" the Editor's opinion that there is no "distinction" between models and real systems.
The Hubbard model, in its original version $\hat{H} = -t \sum_{ij \sigma} \hat{c}^\dag_{i \sigma} \hat{c}_{j \sigma} + U \sum_i \hat{n}_{i \uparrow} \hat{n}_{i \downarrow}$, is a simplification of the reality thanks to a discretization of the space $r \in \mathrm{I\!R} \to i \in \mathrm{I\!I}$ and to a reduction of the range of the interaction, from the longest-range interaction in nature, the Coulombian $1/r$, to a short-range interaction $U \hat{n}_{i \uparrow} \hat{n}_{i \downarrow}$ local to a site $i$. It is true that the site $i$ can be thought both as an effective discretized point in space $r_i$, but also as a more spread-in-space localized orbital $\phi_i(r)$. However along this way the interaction and the model become dependent on the choice.of the basis set. Suppose that we have two electrons in two different localized orbitals $\phi_i(r)$ and $\phi_j(r)$: they won't interact if we choose as basis set of the Hubbard model the $\{\phi_i\}$, but they would interact if we choose another basis set, say Wannier functions $w_k(r)$: the two wavefunctions can both have a contribution on some same Wannier function, $\langle w_k | \phi_i \rangle \ne 0 \quad \langle w_k | \phi_j \rangle \ne 0$ and so interact. In principle the physics must not depend on the chosen basis-set, one should find the same result using different basis-set in the calculation. But here it is not the case, we have a model presenting a basis-set dependent physics.
Second and main point: the raison d'être of this model is the replacement of the many-body long-distance interaction by a local interaction (or an interaction local with respect to a given basis set $\phi_i$) in the hope to have a viable simplification of the many-body problem. Two electrons interact only if they are both in the same site $r_i$ or in the same state $\phi_i$, and do not interact at all if they are in two different sites or states $\phi_i$ and $\phi_j$. Even though the two sites are close or the distributions associated to the two orbitals $\phi_i , \phi_j$ are close in space or even overlap. Which is not the case in the real world. And it cannot become the case even by playing with the range/localization of the orbitals. One can extend the range of the interaction by making the model more complex and extending the interaction to next-neighbors sites/orbitals, but the long-range nature of the Coulomb interaction is such that a finite-rank of next-neighbours cannot really capture it. Note that, along the way to be more and more realistic, a model can become more complex than the reality.
On the other hand, it is true that the Hubbard model proved to be a good representation of many-body systems with short-range interaction, like for example cold atoms systems.
The Editor proposal for a Multiband Hubbard Monomer model to simulate the real helium atom is interesting indeed!
However, as above, how physical or unphysical this model turns out to be in particular question of: which interaction are you considering among the on-site multi bands? Are the electrons on different bands all interacting by the same strength U? Or are electrons interacting by U only if they are on the very same band with different spins? I also guess that the choice of the bands is not a marginal one. Probably, to simulate the helium atom, it would be more judicious to choose not really Bloch bands, and rather hydrogenic Z=2 single-particle atomic orbitals Anyway, I encourage the Editor to develop such model really in the spirit of amazement discussed above: this is interesting per se. However at the end a comparison on the excitation spectrum between the model and the exact Hylleraas solution, is unavoidable if she/he wishes to validate the model as a faithfull representation of the reality. Notice that the model can safely be validated only against an exact solution, that is a solution like Hylleraas that coincides to any level of significant digits to the experiment. Indeed, if we only had the experiment, the model could have found an excitation spectrum with a level very close to, say, the experimental level at 19.8 eV that Hylleraas unambiguosly interpreted as $2^3\!S$ state. But the model could have interpreted that level very differently, I don't know, for example as a lower Hubbard band!