Hamiltonian density of a spin Hamiltonian A many-body Hamiltonian in second quantization is written as
$$
H = \int d\vec r \Psi^\dagger_{\vec r} H_1 (r) \Psi_{\vec r}
$$
where $H_1(r)$ is one-body Hamiltonian. For example, for a non-interacting electronic system $H_1(r) = -\frac{\hbar^2\nabla^2}{2m}$. So, we get
$$
H = -\frac{\hbar^2}{2m} \int d\vec r \Psi^+_{\vec r} \nabla^2 \Psi_{\vec r} =\int d\vec r  \mathcal{E}
$$
here Hamiltonian (or energy) density $\mathcal{E}$ is taken as
$$
\mathcal{E} = -\frac{\hbar^2}{2m}\Psi^\dagger_{\vec r} \nabla^2\Psi_{\vec r} \ .
$$
Now, let we have a many-spin model, written in tight-binding notation as
$$
H = J\sum_{ij}\textbf{S}_i\cdot \textbf{S}_j \ .
$$
In this notation how do we write Hamiltonian density?
 A: 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})$.
