You can try just working out the matrix elements between the different localised states, i.e. $a_1^{\dagger}\ldots a_m^{\dagger}|0\rangle$ etc. This may be rather tedious and require you to carefully keep track of phases due to the fermionic anti-symmetry. Probably a better way is to use the Jordan-Wigner transformation, which gives you an explicit matrix representation of the fermion operators. In 2D this representation of the Hamiltonian looks ugly, but it is still perfectly useful. You will be restricted to a fairly small lattice size, i.e. 10-20 sites depending on your computer and the efficiency of your code. It will be more efficient if you use the fermion number conservation to diagonalise within each number sector separately.
In the particularly simple case of a single particle on the lattice, the problem becomes rather easy, since exchange statistics are no longer relevant. In this case, you make the replacement
$$a_{i,j}^{\dagger}a_{m,n} \to |i,j\rangle\langle m,n|,$$
where
$$|i,j\rangle = a_{i,j}^{\dagger}|0\rangle.$$
To obtain the energy gap in the thermodynamic limit, one should take $N\rightarrow\infty, V\rightarrow\infty$ where $N$ is the number of atoms and $V$ the system size, but hold $N/V$ (i.e. density) fixed. In your case, it means simply taking $N$ to $\infty$ is enough and keep all $t, u, \alpha$ fixed for now. This tight-binding model can be solved exactly. First of all, observe that the first two terms are identity if the spin part is concerned, so one is free to choose any spin quantization axe. Let us for simplicity just rotate the spin axes such that the last term is proportional to $\sigma_z$. Assuming periodic boundary condition $a_{j+N}=a_j$, you can go to the momentum space $a_j=\frac{1}{\sqrt{N}}\sum_k a_k e^{ikj}$, the Hamiltonian becomes
$H=\sum_k (2t\cos k+u)a_{ks}^\dagger a_{ks}+2\alpha\cos k(a_{k\uparrow}^\dagger a_{k\uparrow}-a_{k\downarrow}^\dagger a_{k\downarrow})$
Basically spin up and down electrons are decoupled. Spin up/down electrons have a dispersion
$E_{\uparrow/\downarrow}(k)=(2t\pm\alpha)\cos k+u$. So there are two bands, and now you can just fill up the bands with the electrons.
Notice that $u$ is (minus) the chemical potential, which tells you basically energy levels should be filled up to $-u$. Notice that the bands for spin up and down electrons have almost the same "shape", so you can easily convince your self that as long as $-u$ crosses the dispersion at all (meaning that the bands are neither empty nor fully occupied), you will end up with a metal which has no gap in the thermodynamic limit.
Just one more comment:I'm guessing you are trying to model 1d electrons with spin-orbit coupling. Usually, in systems like semiconductor nanowires, the spin-orbit coupling (e.g. Rashba) looks like $ia_i^\dagger \sigma_y a_{i+1}+\text{h.c.}$ (notice the $i$ in front), not the one you wrote.
Best Answer
First note that while some many-body Hamiltonians can be written in that form, it is not generally the case - you also need to allow for two-body terms (and potentially higher order terms). Probably you're aware of this (if not, I recommend consulting a book like that of Altland & Simons) but it's important to remember since the spin model you wrote down consists entirely of pairwise interactions.
To obtain the Hamiltonian density you can write something like $$ H = \sum_i H_i, $$ where $H_i$ is a Hamiltonian centered around site $i$. For the Hamiltonian you mention it'd simply be $$ H_i = J\mathbf{S}_i \cdot \sum_j \mathbf{S}_j. $$ Note that this looks just like an effective field acting on $\mathbf{S}_i$, but there's no escaping that the field is due to other spins. That is, $H_i$ cannot be written as a one-body Hamiltonian unless you make a mean-field approximation. Now, regardless of whether you make that approximation or not, interpreting $H_i$ as an energy density makes most sense when the sum is restricted to some set of local terms, e.g. if $j$ and $i$ are nearest neighbors. But I wrote the sum unrestricted here, like in your post. Finally, if you want to, you can take a continuum limit (lattice spacing $a\rightarrow 0$) so that $H_i \rightarrow H(\vec{r})$.