A way to characterize multipartite entanglement in pure states of a spin chain with $n$ sites and local dimension $d$ is by means of the Cayley hyperdeterminant. The latter quantity is a polynomial constructed with the components of the wave function $\psi_{i_1, \dots, i_n}$ which is invariant under local unitary transformation... For spin 1/2 chains (i.e. $d=2$) with $n=2$ and $n=3$ sites, the hyperdeterminant coincides with the concurrence and the tangle respectively. In this paper we consider spin chains with $n=4$ sites where the hyperdeterminant is a polynomial of degree 24 containing around $2.8 \times 10^6$ terms. This huge object can be written in terms of more simple polynomials $S$ and $T$ of degrees 8 and 12 respectively. In this paper we compute $S$, $T$ and the hyperdeterminant for eigenstates of the following spin chain Hamiltonians: the transverse Ising model, the XXZ Heisenberg model and the Haldane-Shastry model. Those invariants are also computed for random states, the ground states of random matrix Hamiltonians in the Wigner-Dyson Gaussian ensembles and the quadripartite entangled states defined by Verstraete et al. in 2002. Finally, we propose a generalization of the hyperdeterminant to thermal density matrices. We observe how these polynomials are able to capture the phase transitions present in the models studied as well as a subclass of quadripartite entanglement present in the eigenstates. (read more)