Context:
In non-relativistic QM and many-body theory, the second quantization formalism allows us to write a Hamiltonian for a many-body system with up to two-body terms as (up to a re-ordering definition) $$\hat H = \sum_{ij} t_{ij} a^\dagger_i a_j + h.c. + \sum_{ijkl} V_{ijkl} a^\dagger_i a^\dagger_j a_k a_l $$ where we have defined creation/annihilation operators (either bosonic or fermionic) for a single-particle basis that defines the occupation number basis $$ |n_1, n_2, n_3, \dotsc> $$
In the case of the free-particle plus Coulomb interaction, using a real-space basis would result in $$ \hat H = \sum_x a_x^\dagger t_x a_x + \text{(Hermitian conjugate of previous term)} + \sum_{xx'} V_{xx'x'x} a_x^\dagger a_{x'}^\dagger a_{x'} a_x $$ where the summations are really integrals, and $t_x$ is something like $-\nabla^2_x$, and $V_{xx'x'x} \equiv 1/|{x-x'}|$ in some units. Spin and particle species indices are omitted, but present.
This describes the Hamiltonian of a many-body Coulomb system. Second quantization is the same thing as a wavefunction method, just using a different basis.
Where is the hydrogen atom in this? If we restrict ourselves to the sector with one proton and one electron, which I presume we do by enforcing $$<N_e> = <N_p> = 1,$$ then we should be able to exactly diagonalize this sector of the Hamiltonian, right? How do we do this; what does it look like in second quantization?
How can we directly use the two-particle eigenstates of the hydrogen atom in first quantization (e.g. directly from the Schrodinger equation) to use as an intrinsically two-particle basis for the second-quantized Coulomb system? All textbook formulations of second quantization seem to just use separable one-particle states, but could we get the many-body Coulomb Hamiltonian to simplify considerably by using a two-particle basis defined by hydrogen (and the repulsive Coulomb) eigenstates?
Comments: My end goal is to see what the many-body Coulomb Hamiltonian looks like using an explicitly correlated two-particle basis composed of hydrogen basis functions (and also repulsive Coulomb eigenfunctions). I don't mean that I want to treat hydrogen as a single-particle state. I want to model bigger molecules in terms of nuclei and electrons. I don't think internal coordinates work here, so staying in Cartesian coordinates is my go-to.
For the many-particle Coulomb problem with both protons and electrons (for simplicity), I believe the Hilbert space where the two-particle basis functions would live are of the form $\mathcal{H}_e \otimes \mathcal{H}_p \oplus \mathcal{H}_e \otimes \mathcal{H}_e \oplus \mathcal{H}_p \otimes \mathcal{H}_p$. In short, we need basis functions for both hydrogen, as well as the eigenfunctions of the repulsive Coulomb interaction (available in Landau and Liftshiz, Quantum Mechanics). These details I'm happy to work out on my own, and seem feasible once I understand the two above questions.
My own attempt at finding an answer led me to finding a couple related papers on two-particle bases in second quantization using Gaussian geminals.
a) Kvasnicka, Vladimir. "Second-Quantization Formalism for Geminals." Croatica Chemica Acta 57.6 (1984): 1643-1660.
b) Sørensen, Lasse Kragh. "The minimum parameterization of the wave function for the many-body electronic Schrodinger equation. I. Theory and ansatz." arXiv preprint arXiv:1910.06633 (2019).
However, although interesting in their own right, didn't give me enough to go on to solve this problem related to hydrogen in second quantization.