Commutator between MPO

Question

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.

Answer 1

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:
https://itensor.github.io/ITensors.jl/stable/examples/MPSandMPO.html#Applying-a-Single-site-Operator-to-an-MPS

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.

Miles

Comments

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
         end
   end
    return MPO(ampo, index)
end

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]
norime!(Orho)
rho[1] = Orho

rhoO = newrho[1] * O
noprime!(rhoO)
newrho[1] = rhoO

comm = norm(rho - newrho,2)

But this still doesn't work correctly. Hmmm

---

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.

---

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

---

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)
end