Sorry for the inconvenience,
And thanks for your response.
edit: opening the squares is not feasible I think, since there is a square term with a double summation. This introduces a square term with N^2 local spin.
The density matrix has the dimension NxN^2 and there is a certain pattern in it that will be defined by the hamiltonian. So, S is the spin at the location (i)'th row and (e+m*N)'th column. So, m, i, and e for just indexing. To be more simple I didn't go into the conceptual background of them. If you think it is necessary I can reveal more about them.
Also, I have some prework for flatten this density matrix and thread this is a 1D ground-state search. But it is not fully done yet and I understand as ITensors can work with 2D structures as well.
Other than that S is the spin itself as 0 or 1 (in this case up/down but in the computation, they should be threaded as 0/1). It has also a conceptual meaning. If you need, I can put references. I tried my best to ask the question as simply as possible and didn't go much about them. However, I myself write down the hamiltonian. So, it is not any kind of copy-paste or something. But of course, I inspired some other works.
For the third question, I want to take a linear combination of the defined 3 hamiltonians. I also need to multiply hamiltonians with a constant.
Thank you again for your response, it is a precious experience for me to discuss with you!