How to store pauli's matrices inside an index of a tensor

Question

Hey all!
I am trying to write down an Hamiltonian of the form Sum(J(i,j,a,b)*S(i,a)*S(j,b)
where S is the spin operator Sx+Sy+Sz and J is just some normal rendom numbers, and i>j .
So my problem is that I can't figure out how to store an array inside an index. For example I want that whan a.b= 1 I will get Sx a,b=2 I will get Sy and for a,b=3 I will get Sz.
and to match a random number to each of them (J ).
Can anyone please help me ?
Thank you, Tomer

Comments

Hi Tomer,
Thanks for the question. I have a few to ask you before I will know how to help.

When you say "write down" a Hamiltonian, what kind of thing do you want to make? Do you want your Hamiltonian to be a single tensor? Or an MPO tensor network (which is what we use to represent Hamiltonians in many simulation methods such as DMRG)?

Are you sure you need to make a tensor S(i,a) with the property you mention? To make a spin interaction of the form you wrote above, you can write all the different terms as separate terms and use our OpSum/AutoMPO system to turn these mathematical terms into an MPO? Would that give you what you need? Pairing the right random numbers with the terms should not be too hard or we can discuss how.

Here is an example of the OpSum system making an MPO from a sum of Hamiltonian terms:
https://itensor.github.io/ITensors.jl/stable/tutorials/DMRG.html

---

What I want to do is to write an Hamiltonian (some kind of an Ising model) which sums over 4 indices. My indices are i>j, alpha, beta.
So my goal is to make a MPO which will have a indices of the site, and a index fo the physical information (pauli matrix). The problem is that I don't know how to do it.

---

Hi, thanks for explaining some more. To me, though, those sound like slightly different or incompatible goals. We can definitely help you to make a Hamiltonian MPO that represents the sum over 4 indices. The other goal about an MPO that has (perhaps?) an extra non-standard index for the kind of Pauli matrix might be a different goal and I’m not sure that the result would technically be an MPO? But perhaps I don’t understand your goal in that case.

To make an MPO which makes Heisenberg interactions connecting sites i and j with interaction coefficients J_ij you can use the OpSum and AutoMPO system to do something like this:

os = OpSum()
for i=1:N, j=1:N
  os += J[i,j], “Sz”, i, “Sz”, j
  os += J[i,j]/2, “S+”, i , “S-“, j
  os += J[i,j]/2, “S-“, i, “S+”, j
end

and then turn the OpSum into an MPO as in the ITensor documentation and code examples.

Here in the code above J[i,j] is just a regular Julia array with two dimensions that you would make earlier in the code.

You can generalize the above pattern to 4 indices by just summing over the 4 indices and including four operators in each term if that is what you want. Please let me know if you have a question about that.

If by alpha and beta, though, you just mean the indices which run over x,y,z, for the spin operators to make a Heisenberg S.S interaction, then the above code already includes these in the form of the three terms I put which are (schematically) these terms: Sz*SZ, S+*S-, S-*S+

Best regards,
Miles