Long range spin chains

[jkudlerflam]

I was wondering if it was possible to find ground states for spin chain Hamiltonians with long-range interactions i.e. beyond nearest neighbor.

Currently, this is how I create a nearest neighbor MPO e.g. for the XX model

import quimb.tensor as qtn

builder = qtn.SpinHam(S=1/2) 
builder += -1, 'X', 'X'
builder += -1, 'Y', 'Y'
builder += -h, 'Z'
H = builder.build_mpo(L)

where L is the length of the spin chain and h is a coupling constant.

What I would like to do is create an MPO for a Hamiltonian with arbitrary long range interactions e.g. a term with three sites might look like

builder += -1, 'X', 'X', 'X'

However, when I naively try to do this, I get the following error:

NotImplementedError: 3-body+ terms are not supported yet.

Is there any way to get around this error? Thank you.

[jcmgray]

Since this is for DMRG you just need to find an MPO representation of the Hamiltonian H, i.e. entries of it tensors.

Thinking about finite state machines is one way (the following has quite a clear pic) https://arxiv.org/abs/1909.06341

Alternatively just define the single term MPOs then sum them and compress (both of which quimb can do). https://arxiv.org/abs/1611.02498

However do you just need short-range but 3-body terms like XXX or long range interactions like XIIIIX or arbitrary long range, k-body terms like XIIXIIIZ? The efficiency (bond dimension of the MPO) will depend on this quite a lot.

---

Here's an explicit example for the kind of thing you might have to do for the 3 body term you give:

import numpy as np
import quimb as qu
import quimb.tensor as qtn

L = 5
I = qu.eye(2)
X = qu.pauli('X')

# the dense hamiltonian we are aiming for
H_dense = sum(
    qu.ikron([-X, X, X], dims=[2] * L, inds=[i, i + 1, i + 2])
    for i in range(L - 2)
)

# now we define the MPO tensor
W = np.zeros([4, 4, 2, 2], dtype=complex)
W[0, 0, :, :] = I
W[0, 1, :, :] = X
W[1, 2, :, :] = X
W[2, 3, :, :] = -X
W[3, 3, :, :] = I

Those 5 lines roughly described, as we travel along the chain:

• we can either stay in the '0' (identity) state,
• or transition to state '1' with an X.
• If we transition to state '1' we can then only transition to '2' also with an X,
• and from there only to state '3' also with an X and we'll put the coefficient here as well.
• Once in state '3' we can only stay as the identity.

  
  # define the end pieces
  Wl = W[0, :, :, :]
  Wr = W[:, 3, :, :]

build the MPO

H = qtn.MatrixProductOperator([Wl] + [W] * (L - 2) + [Wr])

check alls working

np.allclose(H.to_dense(), H_dense)

True

Having this stuff automatically done in ``quimb`` would be nice! But I suspect also fiddly to handle the most general cases.

[jcmgray]

Also just to illustrate the other method, here we add our 3 body MPO to the heisenberg MPO, compress it, then find the groundstate:

H2 = qtn.MPO_ham_heis(L)
H3 = H + H2
H3.show()
# │9│9│9│9│
# ●─●─●─●─●
# │ │ │ │ │

H3.compress()
H3.show()
# │4│6│6│4│
# ●─<─<─<─<
# │ │ │ │ │

dmrg = qtn.DMRG2(H3)
dmrg.solve()
# True

dmrg.state.show()
#  2 4 4 2
# >─>─>─>─●
# │ │ │ │ │

dmrg.energy
# (-4.036124752437089+1.096793892211758e-19j)

[jkudlerflam]

Thank you very much for the explanation and references. Ideally, I would have arbitrarily long range k-body terms, but these may be truncated.

I think I understand this procedure though I need to read [1909.06341] carefully.

Say, I want to add a h * XZIX interaction. Then, I would instead create the MPO as:

# now we define the MPO tensor
W = np.zeros([5, 5, 2, 2], dtype=complex)
W[0, 0, :, :] = I
W[0, 1, :, :] = h*X
W[1, 2, :, :] = Z
W[2, 3, :, :] = I
W[3, 4, :, :] = X
W[4, 4, :, :] = I

Is this correct?

[jcmgray]

That looks right. Complications are:

1. If you need non-translational invariance you need to take care placing the operators. If it's just a matter of varying coefficients e.g. h then that's pretty simple though.
2. If you're mixing terms then making if efficient takes a bit of care. E.g. to add a j XX term to the above you could do:

   W = np.zeros([6, 6, 2, 2], dtype=complex)  # need one more bond dim
   W[0, 0, :, :] = I
   #
   W[0, 1, :, :] = h*X
   W[1, 2, :, :] = Z
   W[2, 3, :, :] = I
   W[3, 5, :, :] = X
   # a completely separate 'branch' starting from the '0' state
   W[0, 4, :, :] = j*X
   W[4, 5, :, :] = X
   #
   W[5, 5, :, :] = I

   or you could re-use the first X and branch from the '1' state instead:

   W = np.zeros([5, 5, 2, 2], dtype=complex)  # don't need one more bond dim
   W[0, 0, :, :] = I
   W[0, 1, :, :] = X
   #
   W[1, 2, :, :] = Z
   W[2, 3, :, :] = I
   W[3, 4, :, :] = h * X
   #
   W[1, 4, :, :] = j * X
   #
   W[4, 4, :, :] = I

   Alternatively, you can just define them all separately, sum them, then compress. But this might get slow for many terms and is not guaranteed to find the optimal bond dimension.

In case he has any recommendations @mattorourke17 I think knows a lot more than me about these things.

[jkudlerflam]

Thank you, that's very helpful.