So what you're kind of orbiting around is a central uncertainty about dealing with multiparticle systems, and that's indeed important because you do indeed have this problem of "now what is my wavefunction, if it does not have a fixed set of coordinates $\psi(\vec r_1, \vec r_2, \dots \vec r_N)$ because I don't yet know what $N$ is..." etc.
Eventually we'll define field operators which will make some sense of "position coordinates" again, but for right now you will just have to let that uncertainty go. Just be willing to do what's algebraically allowed in the bra/ket formalism.
In this case, for example, there is clearly a $|2\rangle\langle 2|$ operator which counts the number of particles in state 2. We know that this is well-defined because we see it in the Hamiltonian and you told us the Hamiltonian was well-defined, so there's no problems here. Furthermore it commutes with the
$|1\rangle\langle 1|$ operator which counts the number of particles in state 1, we know this because you are very comfortable saying "there are N particles in state 2 and 0 particles in state 1," which suggests that they are, as single-particle states, disjoint.
Now you do have to decide whether the wavefunction is exactly the same when you interchange two identical particles, or whether the wavefunction is exactly the negative of what it was, or whether the particles are not fundamentally interchangeable and the wavefunction will always know the difference. That last case seems much much harder to treat than the first two.
If the particles are bosonic and saying "all of them are in state 2" literally means "they're all in this single wavefunction state," then you get to invent bosonic annihilators $\hat a,~\hat b$ such that $[\hat a, \hat b] = 0$ but $[\hat a, \hat a^\dagger] = [\hat b, \hat b^\dagger] = 1.$ To model transport between the two you'll have terms like $\hat c~ \hat b^\dagger ~\hat a + \hat c^\dagger~\hat a^\dagger ~\hat b$ in your Hamiltonian, where usually $\hat c$ will describe some other system that the energy goes off into. And this operator we've been calling $|1\rangle\langle1|$ is actually $\hat a^\dagger ~\hat a.$ Easy peasy.
If the particles are fermionic then we generally need to say that there are secretly a collection of states that we're lumping together as "2" and a collection that we're lumping together as "1", so that gets more complicated, we have to peek at the internal states. Invent a fermionic annihilator for each state, so we don't just have $\hat a$ and $\hat b$ but $\hat a_i$ and $\hat b_i$. Now we have $\{\hat a_i, \hat b_j\} = 0$ while $\{\hat a_i, \hat a_j^\dagger\} = \{\hat b_i, \hat b_j^\dagger\} = \delta_{ij}.$ Then this operator that we've been calling $|1\rangle\langle 1|$ is actually $\sum_i \hat a_i^\dagger~\hat a_i.$ Again for transport you probably want to invent a bunch of bosonic modes with annihilator $\hat m_i$ to absorb the energy differences of jumping between the two modes, so that you have a Hamiltonian like (Einstein summation) $$\hat H = 0~ \hat a^\dagger_i ~\hat a_i + \epsilon~\hat b^\dagger_i ~\hat b_i + \epsilon~\hat m^\dagger_i ~\hat m_i + c_{ijk} ~\hat m_i^\dagger ~\hat a_j^\dagger ~\hat b_k + c_{ijk}^* ~\hat m_i ~\hat b_k^\dagger ~\hat a_j.$$So there's no deep mystery in this thing once you abandon this hope of having the explicit wavefunction over all $N$ spatial coordinates written out for you.