"Single particle spectrum can be easily deduced from the above Hamiltonian"
This statement is a bit misleading. Your equation $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.
Let's consider a general (single-particle) tight-binding Hamiltonian:
$$h=\sum_{ij}a_i^{\dagger}t_{ij}a_j$$
The index $i$ spans over all the single-particle degrees of freedom. For example, $i$ may simply label the sites of a lattice; or if there is more than one orbital $a$ per unit cell $x$, then $i$ should be understood to represent a combined index, i.e $i=(x,a)$.
If the Hamiltonian $h$ is translationally invariant -- either an infinite lattice or wrapped on a torus -- then the eigenenergies are functions of the crystal momentum $k$. As a practical matter, we can determine this single-particle spectrum by diagonalizing the matrix $t_{ij}$, which provides a unitary transformation $\Omega$
$$a_{i}^{\dagger} = \sum_{k,\mu}b_{\mu}(k)^{\dagger}\Omega_{\mu i}(k) ,$$
which takes $h$ into the band basis. (I've written the band basis creation/annihilation operators as $b$'s.) I want to stress here that $\Omega$ is a rotation of the single-particle basis, into one which diagonalizes $h$, i.e
$$h = \sum_{k,\mu}b_\mu(k)^{\dagger}\epsilon_{\mu}(k)b_{\mu}(k) .$$
The index $\mu$ labels the bands, of which there are $N=\text{dimension}(\mu)$ states at every $k$. If there is only one orbital per unit cell, then the tight-binding model will have only 1 band, i.e $N=1$. More generally however, if we remember $i=(x,a)$, then there are $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.
Since we are interested in a many-body problem let's introduce some general interactions:
$$V=\sum_{ijkl}a_i^{\dagger}a_j^{\dagger}V_{ijkl}a_ka_l$$
such that our total Hamiltonian is $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."
This statement might come off as confusing, since "lowest band" refers to a single-particle state, while $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$ 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$ 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$ operators with $\hat{a}$ operators. The effective Hamiltonian for the system projected onto $\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}(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.
WARNING: One should be careful and justify the projection, seeing that it involves throwing out information about other bands in the system. There is a pretty typical hierarchy in what people refer to as "flat/narrow band" system: (1) interactions and thermal energy << band gap, for the projection to be justified; and (2) interactions >> bandwidth, for kinetic energy to be considered as a perturbation in the strong coupling problem.
If (1) is not satisfied, then interactions may substantially mix the bands, and thus the error coming from projecting out those bands will be high. This is why the Renormalization Group (RG) is so popular (2). Ultimately, at the end of the day, RG and projection have the same goal -- create an effective Hamiltonian. The difference being that RG does not throw out information about the system, as in the case of the projection. Instead, undesired modes are averaged over, which in turn leads to an effective Hamiltonian which accounts for the presence of those modes. No thermodynamic information is lost in the sense that the partition function is unchanged. This is not true for a straight projection.
Sorry for the lengthy explanation! I thought it would make sense to cover some important facts regarding the practicalities of projection, at the very least for the sake of a future googler who might come here looking for related answers. I hope this explanation kicks you in the right direction!