Multiple sources, including the prodigal paper on the topic, "Theory of Hyperfine Structure" by Charles Schwartz, state that the hyperfine interaction can be written
$$ \hat{H}_{hfs} = \sum_k \mathbf{T}^{(k)}_n \cdot \mathbf{T}^{(k)}_e $$
where the first of these is a rank-$k$ tensor operator acting on nuclear degrees of freedom, and the second a same-rank operator acting on electronic degrees of freedom.
My feeling is that, if you're defining the hyperfine interaction to be all of the interactions between the nuclear and electronic degrees of freedom, this is the most general form of the interaction you can write and still get a scalar operator in the end. Is this correct, and is it the only justification for this form? If not, how can we derive this?
Further, we know that $k=1$ is the magnetic dipole term, $k=2$ the electric quadrupole term, and so on. But both the electric and magnetic multipole expansions in principle have all of these terms; how do we know that the electric dipole, magnetic quadrupole, etc. terms are zero? In a first-order approximation, can we treat this Hamiltonian as diagonal, take expectation values, and use parity arguments to say that this Hamiltonian is equivalent to one with only the nonzero terms? If so, what is the parity of the nuclear spin operator $\hat{I}$ and other relevant operators?