Expectation value of Hamiltonian of Heisenberg spin chain

The following is Hamiltonian of Heisenberg spin half chain mapped onto Hard core bosons by Holestein-Primakoff transformation.

$$H = t\sum_{j}b_{j}^{\dagger}b_{j}+h.c + V\sum_{j}n_{j}n_{j+1}+\sum_{j}h_{j}n_{j}$$

$b^{\dagger},b$ and $n$ represent bosonic creation, annihilation and density operator. $h_{j}$ is random onsite magnetic field.Note that j varies from 1 to N which represents the number of sites in a chain.

I have some questions regarding this model:

1. On what kind of states does this Hamiltonian can operate? For instance,for single site the only valid states are |1⟩ and |0⟩ which represent spin up and spin down states. What are the valid states for N sites? Is it possible to have |2⟩,|3⟩,....,|N⟩ states or Is it N spin up and N spin down states?
2. How can I find the expectation value of this Hamiltonian numerically? The operators obey following relationships with spin operators but I couldn't come up with any way of using them to calculate expectation value of Hamiltonian. $$S_{j}^{\dagger}= b_{j}^{\dagger}\sqrt{1-n_{j}}\notag\\[1em] S_{j}= b_{j}\sqrt{1-n_{j}}\notag\\[1em] S_{j}^{z}= n_{j} + 1/2$$

Thank you.

This is a many-body hamiltonian, since the indexed operators act on each "substate" of the entire system. Imagine a total state A represented by a ket vector; then this is actually composed of $N$ states each in its corresponding Hilbert space. The total Hilbert space which is the inner product of each subspace $h_1\otimes h_2\otimes...h_n$ is totally described by A. So to sum it up, this kind of Hamiltonian will always act on states similar to A.
To find the expectation values of H, first you need to know on what basis are you working on. For example, in problems such as nuclear structure, where rotations and concept of angular moments are important, the Wigner $D_2$ base is relevant, and this you can construct a state A totally known, which will help you solve the Schrodinger equation.