Tensor Product Explanation I'm currently doing a research project involving 3 particle spins and have developed a simple function for the Hamiltonian:

I understand how to code my work but the physics behind it is unfamiliar territory to me. Can someone please explain what exactly the purpose of the tensor product is and how it relates to particle spins. I also am currently learning about the Bloch Sphere in relation to qubit states which has become more clear to me now but I am just unsure about the function of tensor product when performed on 2 matrices.
If anyone could point me in the direction of an explanation to tensor product i.r.t particle spin matrices I'd apreciate it
 A: To start with, I just want to clarify that np.kron implements the Kronecker product of matrices, which is basically the matrix version of the tensor product.  So you can think of it as different language for the same basic idea.
Now, the reason you take tensor products in quantum mechanics is because you are dealing with multiple particles at the same time.  In your case, your code is taking tensor products between three spin matrices, each of which represents a single particle, which tells me that you are dealing with three particles at the same time.
You might wonder why you can't just store three separate spin matrices — one for each particle.  And the reason for that is one of the big ideas in quantum mechanics: entanglement.  That means that the system cannot be described completely by describing each particle separately.  In your case, you'll see that the Hamiltonian can't be written as the sum of the energies for the separate particles; you really have to look at all three at the same time.
If you're just getting used to this idea, it might be easier to work through some examples with just two particles first, and then generalize to more particles.  The easiest entangled state for a two-particle system is $E = S_x \otimes S_x + S_y \otimes S_y$ (where I'm ignoring normalization just to make the point).  If you evaluate this with numpy you get
>>> np.kron(Sx, Sx) + np.kron(Sy, Sy)
array([[0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j]])

If you try to make up the most general possible state $A$ for the first system and $B$ for the second, and evaluate the Kronecker product, you'll quickly see that there's no possible combination such that $E = A \otimes B$.  Specifically, the upper-right and lower-left corners must be nonzero in $A \otimes B$ if our two entries of $1$ above are to appear.  Instead, you really need to take the sum of separate states "tensored" together — something more like $E = A \otimes B + C \otimes D$, which is exactly the form of our definition of $E$ above.
There's a nice discussion of all this that seems to be at the right level here in the MIT open courseware.
