Single"Single particle spectrum can be easily deduced from the above HamiltonianHamiltonian"
This statement is a bit misleading. Your equation H$H$ is a many-body Hamiltonian involving interactions (operators quartic in $a$'s, or likewise quadratic in $n$'s) and single-particle hopping (first part that is quadratic in $a$'s). It is only the single-particle hopping that defines the single-particle spectrum in question.
The index 'i'$i$ spans over all the single-particle degrees of freedom. For example, 'i'$i$ may simply label the sites of a lattice; or if there is more than one orbital 'a'$a$ per unit cell 'x'$x$, then 'i'$i$ should be understood to represent a combined index, i.e 'i=(x,a)'$i=(x,a)$.
If the Hamiltonian h$h$ is translationally invariant -- either an infinite lattice or wrapped on a torus -- then the eigenenergies are functions of the crystal momentum 'k'$k$. As a practical matter, we can determine this single-particle spectrum by diagonalizing the matrix t_ij$t_{ij}$, which provides a unitary transformation \Omega$\Omega$
,$$a_{i}^{\dagger} = \sum_{k,\mu}b_{\mu}(k)^{\dagger}\Omega_{\mu i}(k) ,$$
which takes h$h$ into the band basis. (I've written the band basis creation/annihilation operators as b's$b$'s.) I want to stress here that \Omega$\Omega$ is a rotation of the single-particle basis, into one which diagonalizes h$h$, i.e h equals
.$$h = \sum_{k,\mu}b_\mu(k)^{\dagger}\epsilon_{\mu}(k)b_{\mu}(k) .$$
The index '\mu'$\mu$ labels the bands, of which there are N=dimension(\mu)$N=\text{dimension}(\mu)$ states at every k$k$. If there is only one orbital per unit cell, then the tight-binding model will have only 1 band, i.e N=1$N=1$. More generally however, if we remember 'i=(x,a)'$i=(x,a)$, then there are N=dimension(a)$N=\text{dimension}(a)$ bands. So that the number of bands depends on the number of orbitals in the unit cell. At this point, we have only rotated our basis, and not projected.
$$V=\sum_{ijkl}a_i^{\dagger}a_j^{\dagger}V_{ijkl}a_ka_l$$
such that our total Hamiltonian is H=h+V$H=h+V$,
.$$H=\sum_{ij}a_i^{\dagger}t_{ij}a_j+\sum_{ijkl}a_i^{\dagger}a_j^{\dagger}V_{ijkl}a_ka_l .$$
To quote from the paper you referenced (just before their Eqn 2), "The projection onto the lowest band realizes a new effective problem [...] which derives from the density-density interactions of the microscopic Hamiltonian."
"The projection onto the lowest band realizes a new effective problem [...] which derives from the density-density interactions of the microscopic Hamiltonian."
This statement might come off as confusing, since "lowest band" refers to a single-particle state, while H$H$ is a many-body Hamiltonian involving interactions. More technically, we are constraining ourselves to that section of Fock space where the number operator for any band \mu$\mu$ is zero except for the band of interest. Practically, this can be achieved by defining a new "projected-a" operator
,$$\hat{a}_i^{\dagger}\equiv\mathcal{P}_{\mu_o}a_i^{\dagger}\mathcal{P}_{\mu_o}=\sum_{k}b_{\mu_o}(k)^{\dagger}\Omega_{i,\mu_o}(k),$$
where '\mu_o'$\mu_o$ here references the band of interest. (We can also project onto a subset of bands similarly.) We can then project our Hamiltonian by replacing all 'a'$a$ operators with '\hat{a}'$\hat{a}$ operators. The effective Hamiltonian for the system projected onto '\mu_o'$\mu_o$ is
$$H_{\text{effective}(\mu_o)}=\sum_{k}b_{\mu_o}^{\dagger}(k)\epsilon_{\mu_o}(k)b_{\mu_o}(k)$$
$$+\sum_{k_1k_2k_3k_4}b_{\mu_o}(k_1)^{\dagger}b_{\mu_o}(k_2)^{\dagger}\Big(\sum_{ijkl}\Omega^*_{i\mu_o}(k_1)\Omega^*_{j\mu_o}(k_2)V_{ijkl}\Omega_{k\mu_o}(k_3)\Omega_{l\mu_o}(k_4)\Big)b_{\mu_o}(k_3)b_{\mu_o}(k_4)$$
If the band is truly flat, i.e if \epsilon_{\mu_o}=constant$\epsilon_{\mu_o}(k)=\text{constant}$, then the projected kinetic energy is simply a constant shift in the energy, and can be gauged out, leaving only the interacting part. If it has some dispersion, then one needs to ask if the interactions are much larger than its bandwidth, in which case we could solve the interacting problem (without the kinetic energy), then consider kinetic energy as a perturbation -- the strong coupling approach. I believe this is the paradigm behind the paper you shared.
This expression I've given here for projected V is pretty general, so you'll have to put in some effort in order to get its projected matrix elements into the form in the paper. Not to shamelessly advertise my own work, but a colleague and I wrote a pedagogical section in our paper in which we detail how this is done in the context of the quantum Hall 1 --(see Section 3 up through Eqn 8 for spin-half particles in continuum). This technique is common for studying quantum Hall physics, including fractional quantum Hall, and is a pretty quintessential use of projected potentials for flat bands. Since the gap between Landau levels (i.e bands) can be made much larger than the interactions, mixing between the lowest Landau level (i.e band) and the next Landau level is suppressed to the point where corrections due to inter-band mixing can be ignored. This justifies the projection.