+1 vote
asked by (190 points)

Is there a systematic way to compute the commutator between two MPO's or an MPO and an ITensor?

For example, give two ITensors A and B the commutator is just AB -BA where * is the contraction operator in ITensor.

If say A is an ITensor and B is an MPO then we can compute AB using the apply function but cant go the other way.

I would like to calculate the Frobenius norm of the commutator between a time-evolved MPO and a single site operator.

Any feedback is much appreciated.

1 Answer

0 votes
answered by (62.7k points)

Good question. While you can use the apply function to apply a single-site operator to an MPO, there is a more "manual" way to do this also just using ITensor contractions and the fact that an MPO has an interface like an array of tensors.

A very analogous case is the manual method for applying a single-site operator to an MPS, which is given in detail here:

Following very similar steps, you can apply a single-site operator to an MPO, by paying careful attention to the prime levels of the MPO site indices and of the operator you are applying, and adjusting these prime levels appropriately.

Then by adjusting the prime levels in a different way, you can make the operator act on the "bra" indices of the MPO tensor versus the "ket" indices, which would be like doing BA instead of AB.

Finally, if you take the resulting two MPOs and subtract them, you will have the commutator.

Please try that out, first on paper then in code, and let us know if you have more questions along the way.


commented by (190 points)
edited by
Yeah I figured that was the way to go. Here is what I have been trying most of the day.

function initial(N, index)
     ampo = OpSum()
     for j=1:N
         if j<= Int(0.5*N)
           ampo += 1, "Sz", j
         elseif j > Int(0.5*N)
            ampo += -1, "Sz",j
    return MPO(ampo, index)

This defines an MPO that is a sum of Sz on each site with a weight. Now suppose I want to find the commutator of this mpo with a single-site operator.

N = 2;
s = siteinds("S=1/2", N)
rho0 = initial(N, s)
O = op("Sz", s[1])

Suppose, now we want to take find the commutator on the first site.

rho = copy(rho0)

newrho = dag(prime(rho, "Site"))
newO = dag(prime(O, "Site"))

Orho = newO * rho[1]
rho[1] = Orho

rhoO = newrho[1] * O
newrho[1] = rhoO

comm = norm(rho - newrho,2)

But this still doesn't work correctly. Hmmm
commented by (62.7k points)
Great. So this would be a case where you should draw out on paper the diagrams for each object (the MPO and the operator O) and/or print them out using println or @show and then verify that the prime levels and indices and things are matching the way you expect and accomplishing the contractions that you want.

One detail here is that you don’t need to prime the whole MPO, since you are just going to act on one of the sites. So it could be easier to just work only with that MPO tensor, only altering its prime levels (if at all) and then when it’s “done” and back into the expected form (one site with prime level zero, one site with prime level one) then you can put it back into the MPO.
commented by (190 points)
edited by
I think I found one approach but it is not the manual way which I would like to understand. But since I want to contract a single site at a time if I have my MPO called A and my single site operator B. I can do:

C = A[site]
B = op("Sz", index[site])

BA = apply([B],C)
AB = apply([C],B)

com = BA - AB

val = sqrt(tr(dag(com)*com)) #frobenius norm
commented by (190 points)
Hey Miles,

Is this the appropriate way to compute the norm of commutator at the prime level between an MPO and single-site operator?

function MPO_commutation(B::MPO, x::String, index, n::Int)
    # orthogonalize the input MPO to the site for commutation
    rho_MPO = orthogonalize(B, n)
    #define the single-site operator to perform commutation
    Ops = op(x, index, n)
    # act operator on MPO
    New_rho = dag(prime(Ops))*rho_MPO[n]
    #restore prime levels
    New_rho = swapprime(New_rho, 2=>1)
    # act MPO on operator
    New_op = dag(prime(rho_MPO[n], "Site")) * Ops
    #restore prime levels
    New_op = swapprime(New_op, 2=>1)
    return norm(New_rho - New_op)
Welcome to ITensor Support Q&A, where you can ask questions and receive answers from other members of the community.

Formatting Tips:
  • To format code, indent by four spaces
  • To format inline LaTeX, surround it by @@ on both sides
  • To format LaTeX on its own line, surround it by $$ above and below
  • For LaTeX, it may be necessary to backslash-escape underscore characters to obtain proper formatting. So for example writing \sum\_i to represent a sum over i.
If you cannot register due to firewall issues (e.g. you cannot see the capcha box) please email Miles Stoudenmire to ask for an account.

To report ITensor bugs, please use the issue tracker.